polyserial function

Polyserial Correlation

Computes the polyserial correlation (and its standard error) between a quantitative variable and an ordinal variable, based on the assumption that the joint distribution of the quantitative variable and a latent continuous variable underlying the ordinal variable is bivariate normal. Either the maximum-likelihood estimator or a quicker ``two-step'' approximation is available. For the ML estimator the estimates of the thresholds and the covariance matrix of the estimates are also available.

polyserial(x, y, ML = FALSE, control = list(), 
  std.err = FALSE, maxcor=.9999, bins=4, start, thresholds=FALSE)

Arguments

• x: a numerical variable.
• y: an ordered categorical variable; can be numeric, logical, a factor, an ordered factor, or a character variables, but if a factor, its levels should be in proper order, and the values of a character variable are ordered alphabetically.
• ML: if TRUE, compute the maximum-likelihood estimate; if FALSE, the default, compute a quicker ``two-step'' approximation.
• control: optional arguments to be passed to the optim function.
• std.err: if TRUE, return the estimated variance of the correlation (for the two-step estimator) or the estimated covariance matrix of the correlation and thresholds (for the ML estimator); the default is FALSE.
• maxcor: maximum absolute correlation (to insure numerical stability).
• bins: the number of bins into which to dissect x for a test of bivariate normality; the default is 4.
• start: optional start value(s): if a single number, start value for the correlation; if a list with the elements rho and thresholds, start values for these parameters; start values are supplied automatically if omitted, and are only relevant when the ML estimator or standard errors are selected.
• thresholds: if TRUE (the default is FALSE) return estimated thresholds along with the estimated correlation even if standard errors aren't computed.

Returns

If std.err or thresholds is TRUE, returns an object of class "polycor" with the following components: - type: set to "polyserial".

• rho: the polyserial correlation.

• cuts: estimated thresholds for the ordinal variable (y), for the ML estimator.

• var: the estimated variance of the correlation, or, for the ML estimator, the estimated covariance matrix of the correlation and thresholds.

• n: the number of observations on which the correlation is based.

• chisq: chi-square test for bivariate normality.

• df: degrees of freedom for the test of bivariate normality.

• ML: TRUE for the ML estimate, FALSE for the two-step estimate.

Othewise, returns the polyserial correlation.

Details

The ML estimator is computed by maximizing the bivariate-normal likelihood with respect to the thresholds for $y$ ($\tau^y[j], j = 1,\ldots, c - 1$) and the population correlation ($\rho$). The likelihood is maximized numerically using the optim function, and the covariance matrix of the estimated parameters is based on the numerical Hessian computed by optim.

The two-step estimator is computed by first estimating the thresholds ($\tau^y[j], j = 1,\ldots, c - 1$) from the marginal distribution of $y$. Then if the standard error of $\rho hat$ is requested, the one-dimensional likelihood for $\rho$ is maximized numerically, using optim if standard errors are requested; the standard error computed treats the thresholds as fixed. If the standard error isn't request, $\rho hat$ is computed directly.

References

Drasgow, F. (1986) Polychoric and polyserial correlations. Pp. 68--74 in S. Kotz and N. Johnson, eds., The Encyclopedia of Statistics, Volume 7. Wiley.

Author(s)

John Fox jfox@mcmaster.ca

See Also

hetcor, polychor, print.polycor, optim

Examples

if(require(mvtnorm)){
    set.seed(12345)
    data <- rmvnorm(1000, c(0, 0), matrix(c(1, .5, .5, 1), 2, 2))
    x <- data[,1]
    y <- data[,2]
    cor(x, y)  # sample correlation
    }
if(require(mvtnorm)){
    y <- cut(y, c(-Inf, -1, .5, 1.5, Inf))
    polyserial(x, y)  # 2-step estimate
    }
if(require(mvtnorm)){
    polyserial(x, y, ML=TRUE, std.err=TRUE) # ML estimate
    }