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

#### Documented in polyserial

```# last modified 2021-12-07 by J. Fox

polyserial <- function(x, y, ML=FALSE, control=list(), std.err=FALSE,
maxcor=.9999, bins=4, start, thresholds=FALSE){
f <- function(pars){
rho <- pars[1]
if (abs(rho) > maxcor) rho <- sign(rho)*maxcor
cts <- if (length(pars) == 1) c(-Inf, cuts, Inf)
else c(-Inf, pars[-1], Inf)
if (any(diff(cts) < 0)) return(Inf)
tau <- (matrix(cts, n, s+1, byrow=TRUE) - matrix(rho*z, n, s+1))/
sqrt(1 - rho^2)
- sum(log(dnorm(z)*(pnorm(tau[cbind(indices, y+1)]) - pnorm(tau[cbind(indices, y)]))))
}
if (!is.numeric(x)) stop("x must be numeric")
valid <- complete.cases(x, y)
x <- x[valid]
y <- y[valid]
z <- scale(x)
tab <- table(y)
n <- sum(tab)
s <- length(tab)
if (s < 2) {
warning("y has fewer than 2 levels")
return(NA)
}
if (sum(0 != tab) < 2){
warning("y has fewer than 2 levels with data")
return(NA)
}
nzeros <- sum(zeros <- 0 == tab)
if (nzeros > 0){
warning("the following ", if (nzeros == 1) " level" else " levels", " of y",
if (nzeros == 1) " has" else " have", " no cases: ", names(tab)[zeros])
}
indices <- 1:n
cuts <- qnorm(cumsum(tab)/n)[-s]
y <- as.numeric(as.factor(y))
rho <- sqrt((n - 1)/n)*sd(y)*cor(x, y)/sum(dnorm(cuts))
if (!missing(start) && (ML || std.err)) {
if (is.list(start)){
rho <- start\$rho
cuts <- start\$thresholds
} else {
rho <- start
}
if (!is.numeric(rho) || length(rho) != 1)
stop("start value for rho must be a number")
if (!is.numeric(cuts) || length(cuts) != s - 1)
stop("start values for thresholds must be ", s - 1, "numbers")
}
if (abs(rho) > maxcor) {
warning("initial correlation inadmissible, ", rho, ", set to ", sign(rho)*maxcor)
rho <- sign(rho)*maxcor
}
if (ML) {
result <- optim(c(rho, cuts), f, control=control, hessian=std.err)
if (result\$par[1] > 1){
result\$par[1] <- maxcor
}
else if (result\$par[1] < -1){
result\$par[1] <- -maxcor
warning(paste("inadmissible correlation set to -", maxcor, sep=""))
}
if (std.err){
chisq <- chisq(y, z, result\$par[1], result\$par[-1], bins=bins)
df <- s*bins - s  - 1
result <- list(type="polyserial",
rho=result\$par[1],
cuts=result\$par[-1],
var=solve(result\$hessian),
n=n,
chisq=chisq,
df=df,
ML=TRUE)
class(result) <- "polycor"
return(result)
} else if (thresholds){
result <- list(type="polyserial",
rho=result\$par[1],
cuts=result\$par[-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(rho, f, control=control, hessian=TRUE, method="BFGS")
if (result\$par > 1){
result\$par <- maxcor
}
else if (result\$par < -1){
result\$par <- -maxcor
warning(paste("inadmissible correlation set to -", maxcor, sep=""))
}
chisq <- chisq(y, z, rho, cuts, bins=bins)
df <- s*bins - s  - 1
result <- list(type="polyserial",
rho=result\$par,
cuts=cuts,
var=1/result\$hessian,
n=n,
chisq=chisq,
df=df,
ML=FALSE)
class(result) <- "polycor"
return(result)
}
else if (thresholds){
result <- list(type="polyserial",
rho=rho,
cuts=cuts,
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 Dec. 11, 2021, 3:01 a.m.