[R] tetrachoric correlations

[Timothy R. Johnson]

I have a function that I wrote for R to compute tetrachoric/polychoric
correlations via maximum likelihood. It is given below. I have not tested it
extensively but it seems to work. It could probably be cleaned-up a bit. It
needs the mvtnorm library which is available at CRAN.

Tim Johnson
--
Assistant Professor of Statistics
University of Idaho
Moscow, Idaho 83844-1104

http://www.uidaho.edu/~trjohns

polychor <- function(x, trace = F, se = F) {

 # Purpose: Provides maxiumum likelihood estimates of tetrachoric/
 #  polychoric correlation and thresholds from a contingency table.
 #  Standard errors from the numerical hessian may also be obtained.
 #  The mvtnorm library is necessary for the pmvnorm() function.
 # Note: I have not tested this function extensively. You should
        #  be wary of potential numerical problems in either pmvnorm()
        #  optim().
 # Author: Tim R. Johnson (trjohns at uidaho.edu)
        # Last Revised: Feb 7, 2001

 n <- nrow(x)
 m <- ncol(x)
 tri.row <- matrix(0, n-1, n-1)
 tri.col <- matrix(0, m-1, m-1)
 tri.row[row(tri.row) >= col(tri.row)] <- 1
 tri.col[row(tri.col) >= col(tri.col)] <- 1

 negloglik <- function(theta) {
  row.cut <- c(+Inf, rev(tri.row %*% theta[1:(n-1)]), -Inf)
  col.cut <- c(-Inf, tri.col %*% theta[n:(n+m-2)],    +Inf)
  prb <- matrix(NA, n, m)
  rho <- theta[n + m - 1]
  sig <- matrix(rho, 2, 2)
  diag(sig) <- 1
  for(i in 1:n) {
   for(j in 1:m) {
    prb[i,j] <- pmvnorm(c(0,0), sig,
     c(row.cut[i+1], col.cut[j]),
     c(row.cut[i], col.cut[j+1]))$value
   }
  }
  sum(x * log(prb)) * (-1)
 }

 theta.start <- c(seq(-2, 2, length = n-1), seq(-2, 2, length = m-1),
  runif(1, -0.8, 0.8))

 fit <- optim(theta.start, negloglik, method = "L-BFGS-B",
  lower = c(-Inf, rep(0, n-2), -Inf, rep(0, m-2), -1),
  hessian = se, control = list(trace = trace, REPORT = 1))
 est <- c(tri.row %*% fit$par[1:(n-1)],
   tri.col %*% fit$par[n:(n+m-2)], fit$par[n+m-1])
 ste <- rep(NA, length(est))
 if(se) {
  covmat <- solve(fit$hessian)
  ste[1:(n-1)] <- sqrt(diag(tri.row %*% covmat[1:(n-1),1:(n-1)] %*%
t(tri.row)))
  ste[n:(n+m-2)] <- sqrt(diag(tri.col %*% covmat[n:(n+m-2),n:(n+m-2)] %*%
t(tri.col)))
  ste[n + m - 1] <- sqrt(covmat[n + m - 1, n + m - 1])
  }
 out <- cbind(est, ste)
 dimnames(out) <- list(c(paste("rowcut",1:(n-1),sep=""),
  paste("colcut",1:(m-1),sep=""),"rho"), c("estimate","se"))
 list(loglik = -fit$value, fit = out)
}

----- Original Message -----
From: "Michael Campbell" <M.Campbell at dial.pipex.com>
To: <r-help at stat.math.ethz.ch>
Sent: Thursday, March 15, 2001 4.54 AM
Subject: [R] tetrachoric correlations


> Does anyone have code (preferably S or R) to calculate tetrachoric
correlations?
>
> Thanks
> Michael Campbell
>
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
-.-.-
> r-help mailing list -- Read
http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
>
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._
>

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._