# Constrain Total Variance instead of Residual Variance in Ordinal Threshold Path Mdel

5 posts / 0 new
Offline
Joined: 01/14/2010 - 16:47
Constrain Total Variance instead of Residual Variance in Ordinal Threshold Path Mdel

The path specification for an ordinal model in the documentation constrains the residual variances for the ordinal variables to 1. I would like to instead constrain the total variance for the ordinal variables to 1. But I'm not sure how to accomplish this. When the total variance is constrained to 1, there residual variance is constrained to 1 minus the factor loading squared (the factor loadings are constrained to be between -1 and 1) as I understand it. I was able to accomplish this in the matrix model but would like to be able to do it in the path type model.

Offline
Joined: 07/31/2009 - 15:14

Hi Richard

You could extract the diagonal of the expected covariance matrix and use mxConstraint to equate this to a unit vector. You lose SE's when using mxConstraints, however. Here's an example:

#
#  OpenMx Ordinal Data Example
#  Revision history:
#     Michael Neale 14 Aug 2010
#

require(OpenMx)

#
# Step 2: set up simulation parameters
# Note: nVariables>=3, nThresholds>=1, nSubjects>=nVariables*nThresholds (maybe more)
# and model should be identified
#
nVariables<-3
nFactors<-1
nThresholds<-3
nSubjects<-500
isIdentified<-function(nVariables,nFactors) as.logical(1+sign((nVariables*(nVariables-1)/2) -  nVariables*nFactors + nFactors*(nFactors-1)/2))
# if this function returns FALSE then model is not identified, otherwise it is.
isIdentified(nVariables,nFactors)

mu <- matrix(0,nrow=nVariables,ncol=1)
# Step 3: simulate multivariate normal data
set.seed(1234)
continuousData <- mvtnorm::rmvnorm(n=nSubjects,mu,sigma)

# Step 4: chop continuous variables into ordinal data
# with nThresholds+1 approximately equal categories, based on 1st variable
quants<-quantile(continuousData[,1],  probs = c((1:nThresholds)/(nThresholds+1)))
ordinalData<-matrix(0,nrow=nSubjects,ncol=nVariables)
for(i in 1:nVariables)
{
ordinalData[,i] <- cut(as.vector(continuousData[,i]),c(-Inf,quants,Inf))
}

# Step 5: make the ordinal variables into R factors
ordinalData <- mxFactor(as.data.frame(ordinalData),levels=c(1:(nThresholds+1)))

# Step 6: name the variables
fruitynames<-paste("banana",1:nVariables,sep="")
names(ordinalData)<-fruitynames

thresholdModel <- mxModel("thresholdModel",
mxMatrix("Full", nVariables, nFactors, values=0.2, free=TRUE, lbound=-.99, ubound=.99, name="L"),
mxMatrix("Unit", nVariables, 1, name="vectorofOnes"),
mxMatrix("Zero", 1, nVariables, name="M"),
mxMatrix("Diag", free=T, values=.5, nVariables, nVariables, name="E"),
mxAlgebra(L %*% t(L) + (E), name="impliedCovs"),
mxConstraint(vectorofOnes == diag2vec(impliedCovs)),
mxMatrix("Full",
name="thresholdDeviations", nrow=nThresholds, ncol=nVariables,
values=.2,
free = TRUE,
lbound = rep( c(-Inf,rep(.01,(nThresholds-1))) , nVariables),
dimnames = list(c(), fruitynames)),
mxMatrix("Lower",nThresholds,nThresholds,values=1,free=F,name="unitLower"),
mxAlgebra(unitLower %*% thresholdDeviations, name="thresholdMatrix"),
mxFitFunctionML(),mxExpectationNormal(covariance="impliedCovs", means="M", dimnames = fruitynames, thresholds="thresholdMatrix"),
mxData(observed=ordinalData, type='raw')
)

summary(thresholdModelrun <- mxRun(thresholdModel))
Offline
Joined: 01/24/2014 - 12:15
Paths; constraint not necessary

Richard wants to know how to do it with path specification, though, and not matrix specification. I'm not sure, though, that there's a good way to do it with path specification that doesn't involve creating MxAlgebras and MxConstraints.

Also, with matrix specification, it's possible to do it without using MxConstraints,

thresholdModel2 <- mxModel(
"thresholdModel",
mxMatrix("Full", nVariables, nFactors, values=0.2, free=TRUE, lbound=-.99, ubound=.99, name="L"),
mxMatrix("Unit", nVariables, 1, name="vectorofOnes"),
mxMatrix("Zero", 1, nVariables, name="M"),
mxAlgebra(L %*% t(L), name="Communality"),
mxAlgebra(vec2diag(vectorofOnes - diag2vec(Communality)), name="Uniqueness"),
mxAlgebra(Communality + Uniqueness, name="impliedCovs"),
mxMatrix(
"Full",
name="thresholdDeviations", nrow=nThresholds, ncol=nVariables,
values=.2,
free = TRUE,
lbound = rep( c(-Inf,rep(.01,(nThresholds-1))) , nVariables),
dimnames = list(c(), fruitynames)),
mxMatrix("Lower",nThresholds,nThresholds,values=1,free=F,name="unitLower"),
mxAlgebra(unitLower %*% thresholdDeviations, name="thresholdMatrix"),
mxFitFunctionML(),mxExpectationNormal(covariance="impliedCovs", means="M", dimnames = fruitynames, thresholds="thresholdMatrix"),
mxData(observed=ordinalData, type='raw')
)

though this is slightly more difficult to optimize, since the constraint is only implicit.

Offline
Joined: 07/31/2009 - 15:14
Fix all residuals to 1

Another equivalent approach (same model fit) is simply to fix the residuals to a positive constant like 1. This has the effect of returning unstandardized thresholds and parameters etc. but standardized parameters can be returned for RAM models. That doesn't solve the threshold ordering problem, for which I don't see an alternative except the algebra. Might be good for threshold RAM models to do the ordering internally and automagically.

Offline
Joined: 01/24/2014 - 12:15
one possible way

In the single-factor case, something you could do is as follows. Give parameter labels to the paths that represent your factor loadings and your unique ("residual") variances. Create an MxMatrix that also contains the factor loadings, by assigning its elements the same parameter labels that you gave the loading paths. Also create a diagonal MxMatrix that contains the unique variances, again by use of labels. Create an MxAlgebra that multiplies the matrix of loadings by its transpose, and adds that product to the diagonal matrix of uniquenesses. Finally, using diag2vec(), create an MxConstraint that sets the diagonal elements of the algebra's result to unity.

You can probably see that the sort of constraint you're asking about isn't something that's easy to do with MxPath specification, which is probably why the example fixes the unique variances instead.