fid = gzfile("NHANES-600.csv.gz")
M = read.csv(fid)

PC = array(0, c(20,20))
PLOR = array(0, c(20,20))

for (j1 in 9:28) {
  for (j2 in 9:j1) {

    if (j1 == j2) { next }

    ## Pearson correlation analysis.
    c = cor(M[,j1], M[,j2], use="pairwise.complete.obs")
    fc = 0.5*log((1+c)/(1-c)) 
    n = dim(M)[1] - sum(is.na(M[,j1]) | is.na(M[,j2]))
    zc = sqrt(n-3) * fc
    pc = 2*(1 - pnorm(abs(zc)))
    PC[j1-8,j2-8] = zc

    ## Log odds ratio analysis.
    m1 = median(M[,j1], na.rm=TRUE)
    m2 = median(M[,j2], na.rm=TRUE)
    ii = which(!(is.na(M[,j1]) | is.na(M[,j2])))

    ## Get the 2x2 table of counts.
    N = array(0, c(2,2))
    N[1,1] = sum( (M[ii,j1]>=m1) & (M[ii,j2]>=m2) )
    N[1,2] = sum( (M[ii,j1]>=m1) & (M[ii,j2]<m2) )
    N[2,1] = sum( (M[ii,j1]<m1)  & (M[ii,j2]>=m2) )
    N[2,2] = sum( (M[ii,j1]<m1)  & (M[ii,j2]<m2) )

    LOR = log(N[1,1]) + log(N[2,2]) - log(N[1,2]) - log(N[2,1])
    SE = sqrt(1/N[1,1] + 1/N[1,2] + 1/N[2,1] + 1/N[2,2])
    ZLOR = LOR / SE
    plor = 2*(1-pnorm(abs(ZLOR)))

    PLOR[j1-8,j2-8] = ZLOR
  }
}