# R/polychor.R In polycor: Polychoric and Polyserial Correlations

```# last modified 2022-01-10 by J. Fox

polychor <- function (x, y, ML=FALSE, control=list(), std.err=FALSE,
maxcor=.9999, start, thresholds=FALSE){
f <- function(pars) {
if (length(pars) == 1){
rho <- pars
if (abs(rho) > maxcor) rho <- sign(rho)*maxcor
row.cuts <- rc
col.cuts <- cc
}
else {
rho <- pars[1]
if (abs(rho) > maxcor) rho <- sign(rho)*maxcor
row.cuts <- pars[2:r]
col.cuts <- pars[(r+1):(r+c-1)]
if (any(diff(row.cuts) < 0) || any(diff(col.cuts) < 0)) return(Inf)
}
P <- binBvn(rho, row.cuts, col.cuts)
- sum(tab * log(P))
}
tab <- if (missing(y)) x else table(x, y)
zerorows <- apply(tab, 1, function(x) all(x == 0))
zerocols <- apply(tab, 2, function(x) all(x == 0))
zr <- sum(zerorows)
if (0 < zr) warning(paste(zr, " row", suffix <- if(zr == 1) "" else "s",
" with zero marginal", suffix," removed", sep=""))
zc <- sum(zerocols)
if (0 < zc) warning(paste(zc, " column", suffix <- if(zc == 1) "" else "s",
" with zero marginal", suffix, " removed", sep=""))
tab <- tab[!zerorows, ,drop=FALSE]
tab <- tab[, !zerocols, drop=FALSE]
r <- nrow(tab)
c <- ncol(tab)
if (r < 2) {
warning("the table has fewer than 2 rows")
return(NA)
}
if (c < 2) {
warning("the table has fewer than 2 columns")
return(NA)
}
n <- sum(tab)
rc <- qnorm(cumsum(rowSums(tab))/n)[-r]
cc <- qnorm(cumsum(colSums(tab))/n)[-c]

if (!missing(start) && (ML || std.err)) {
if (is.list(start)){
rho <- start\$rho
rc <- start\$row.thresholds
cc <- start\$col.thresholds
} else {
rho <- start
}
if (!is.numeric(rho) || length(rho) != 1)
stop("start value for rho must be a number")
if (!is.numeric(rc) || length(rc) != r - 1)
stop("start values for row thresholds must be ", r - 1, "numbers")
if (!is.numeric(cc) || length(cc) != c - 1)
stop("start values for column thresholds must be ", c - 1, "numbers")
}

if (ML) {
result <- optim(
c(if (missing(start)) optimise(f, interval=c(-1, 1))\$minimum else rho, rc, cc),
f, control=control, hessian=std.err
)
if (result\$par[1] > 1){
result\$par[1] <- maxcor
warning(paste("inadmissible correlation set to", maxcor))
}
else if (result\$par[1] < -1){
result\$par[1] <- -maxcor
warning(paste("inadmissible correlation set to -", maxcor, sep=""))
}
if (std.err) {
chisq <- 2*(result\$value + sum(tab * log((tab + 1e-6)/n)))
df <- length(tab) - r - c
result <- list(type="polychoric",
rho=result\$par[1],
row.cuts=result\$par[2:r],
col.cuts=result\$par[(r+1):(r+c-1)],
var=solve(result\$hessian),
n=n,
chisq=chisq,
df=df,
ML=TRUE)
class(result) <- "polycor"
return(result)
} else if (thresholds){
result <- list(type="polychoric",
rho=result\$par[1],
row.cuts=result\$par[2:r],
col.cuts=result\$par[(r+1):(r+c-1)],
var=NA,
n=n,
chisq=NA,
df=NA,
ML=TRUE)
class(result) <- "polycor"
return(result)
}
else return(as.vector(result\$par[1]))
}
else if (std.err){
result <- optim(if (missing(start)) 0 else rho,
f, control=control, hessian=TRUE, method="BFGS")
if (result\$par > 1){
result\$par <- maxcor
warning(paste("inadmissible correlation set to", maxcor))
}
else if (result\$par < -1){
result\$par <- -maxcor
warning(paste("inadmissible correlation set to -", maxcor, sep=""))
}
chisq <- 2*(result\$value + sum(tab *log((tab + 1e-6)/n)))
df <- length(tab) - r - c
result <- list(type="polychoric",
rho=result\$par,
row.cuts=rc,
col.cuts=cc,
var=1/result\$hessian,
n=n,
chisq=chisq,
df=df,
ML=FALSE)
class(result) <- "polycor"
return(result)
} else {
rho <- optimise(f, interval=c(-maxcor, maxcor))\$minimum
if (thresholds){
result <- list(type="polychoric",
rho=rho,
row.cuts=rc,
col.cuts=cc,
var=NA,
n=n,
chisq=NA,
df=NA,
ML=FALSE)
class(result) <- "polycor"
return(result)
} else {
return(rho)
}
}
}
```

## Try the polycor package in your browser

Any scripts or data that you put into this service are public.

polycor documentation built on Jan. 12, 2022, 1:08 a.m.