四分格相關係數的標準差及顯著性檢驗

    http://www.john-uebersax.com/stat/tetra.htm 提到可用兩種方法作四分格相關係數的顯著性檢驗,這裏用到的是第一種方法,也就是利用標準差的檢驗.

     <The Tetrachoric Correlation and its Asymptotic Standard Error>一文講到四分格相關係數的標準差爲s:

tetrachoric.se <- function(a, b, c, d) {
  n <- a + b + c + d
  z1 <- qnorm((a + c)/n)
  z2 <- qnorm((a + b)/n)
  r <- psych::tetrachoric(matrix(c(a, b, c, d), 2, 2))$rho
  f1 <- pnorm((z1 - r * z2)/sqrt(1 - r * r)) - 0.5
  f2 <- pnorm((z2 - r * z1)/sqrt(1 - r * r)) - 0.5
  t <- (a + d) * (b + c)/4 + 
    (a + c) * (b + d) * f2 * f2 + 
    (a + b) * (c + d) * f1 * f1 +
    2 * (a * d - b * c) * f1 * f2 +
    (c * d - b * a) * f2 +
    (b * d - a * c) * f1
  s <- n^(-3/2) * sqrt(t)/ 
    mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  r,  r, 1), 2, 2))
  c(n, z1, z2, r, f1, f2, t, s)
}

    可是給出的假定相關係數爲0時的標準差表達式有誤,具體爲下列代碼中的s1,修正爲s2:

tetrachoric.se0 <- function(a, b, c, d) {
  n <- a + b + c + d
  z1 <- qnorm((a + c)/n)
  z2 <- qnorm((a + b)/n)
  r <- 0
  f1 <- pnorm((z1 - r * z2)/sqrt(1 - r * r)) - 0.5
  f2 <- pnorm((z2 - r * z1)/sqrt(1 - r * r)) - 0.5
  t <- (a + d) * (b + c)/4 + 
    (a + c) * (b + d) * f2 * f2 + 
    (a + b) * (c + d) * f1 * f1 +
    2 * (a * d - b * c) * f1 * f2 +
    (c * d - b * a) * f2 +
    (b * d - a * c) * f1
  s <- n^(-3/2) * sqrt(t)/ 
    mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  r,  r, 1), 2, 2))
  s1 <- n^(-5/2) * sqrt((a + b) * (a + c) * (b + d) * (b + c))/ 
    mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  0,  0, 1), 2, 2))
  s2 <- n^(-5/2) * sqrt((a + b) * (a + c) * (b + d) * (d + c))/ 
    mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  0,  0, 1), 2, 2))
  c(n, z1, z2, r, f1, f2, t, s, s1, s2)
}

    其中s爲直接將相關係數設爲0以後求得的標準差, s1爲文章提供的表達式求得(估計印刷錯誤),s2爲修正以後的值,與s的值相同.

    有了相關係數爲0的標準差以後,就能夠檢驗相關係數的顯著性了,p值爲:

2 * pnorm(abs(r)/s2, lower.tail = F)

    寫一個完整的求四分格相關係數顯著性的R語言代碼:

tetrachoric.sig <- function(a, b, c, d) {
  n <- a + b + c + d
  z1 <- qnorm((a + c)/n)
  z2 <- qnorm((a + b)/n)
  r <- psych::tetrachoric(matrix(c(a, b, c, d), 2, 2))$rho
  se <- n^(-5/2) * sqrt((a + b) * (a + c) * (b + d) * (d + c))/ 
    mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  0,  0, 1), 2, 2))
  p <- 2 * pnorm(abs(r)/se, lower.tail = F)
  list(rho = r, p.value = p)
}

tetrachoric.sig <- function(m) {
  n <- sum(m)
  z1 <- qnorm(sum(m[1, ])/n)
  z2 <- qnorm(sum(m[, 1])/n)
  r <- psych::tetrachoric(m)$rho
  se <- n^(-5/2) * sqrt(prod(rowSums(m), colSums(m)))/
  mnormt::dmnorm(c(z1, z2), varcov = matrix(c(1,  0,  0, 1), 2, 2))
  p <- 2 * pnorm(abs(r)/se, lower.tail = F)
  list(rho = r, p.value = p)
}