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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._