# confidence intervals in a common pathway model (with constraint)

6 posts / 0 new
Offline
Joined: 07/31/2009 - 14:25
confidence intervals in a common pathway model (with constraint)

Reviewers have requested I give confidence intervals on estimates from a common pathway model (I am using the umx::umxCP function to build my common pathway model).

When I try and use Hunter's nice mxSE(top.es_std[1,1], model = cp3) function on to request the SE on the specific-e for variable 1), I get the following error:

Model does not have a reasonable Hessian or standard errors.
Model has at least one mxConstraint. This prevented standard error computation.
Try mxCI()


Is there any work-around for this (to use mxSE with a constraint in the model)?

Computing profile-based CIs will take a considerable time…

Can the CP model be implemented without constraints? the constraint used is

    mxConstraint(diagL == nFac_Unit   , name = "fix_CP_variances_to_1"),

Many thanks!

PS: Why do constraints invalidate SEs?

Offline
Joined: 01/24/2014 - 12:15
constraints
Can the CP model be implemented without constraints?

I think so. You could instead just fix the nonshared-environmental variance of the CP to some value (1.0, say). Or, you could fix one of the loadings onto the CP to 1.0, and freely estimate the CP's A, C, and E variances.

Is there any work-around for this (to use mxSE with a constraint in the model)?

PS: Why do constraints invalidate SEs?

As you know, standard errors are calculated as the square roots of the diagonal elements of the inverted Fisher information matrix (which in OpenMx is typically twice the inverted Hessian). The sampling covariance matrix you get from inverting the information matrix simply doesn't reflect the influence of the constraints on the repeated-sampling distribution of the parameter estimates.

You must have missed my recent email to the OpenMx developer list about valid standard errors in the presence of constraints. It turns out it's possible to calculate an asymptotically valid sampling covariance matrix from the usual Hessian, and a matrix providing a basis for the nullspace of the constraint Jacobian matrix. The current plan is to implement this approach in the backend, as part of the MxComputeNumericDeriv and MxComputeStandardError compute steps. Here's a demonstration script:

library(OpenMx)
mxOption(NULL,"Default optimizer","SLSQP")
library(MASS)

#Custom compute plan, to calculate Hessian despite constraints:
plan <- omxDefaultComputePlan()
plan$steps <- list(plan$steps$GD, plan$steps$ND, plan$steps$RD, plan$steps$RE) #MxModel to estimate multinomial proportions from frequencies (something that can be done by hand): m1 <- mxModel( "MultinomialWithLinearConstraints", plan, mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=0.25,labels="pred",name="Pred",lbound=0,ubound=1), mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=0.25,labels="pyellow",name="Pyellow",lbound=0,ubound=1), mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=0.25,labels="pgreen",name="Pgreen",lbound=0,ubound=1), mxMatrix(type="Full",nrow=1,ncol=1,free=T,values=0.25,labels="pblue",name="Pblue",lbound=0,ubound=1), mxAlgebra( -2*(43*log(Pred) + 22*log(Pyellow) + 20*log(Pgreen) + 15*log(Pblue)), name="fitfunc"), mxAlgebra( cbind(-2*43/Pred,-2*22/Pyellow,-2*20/Pgreen,-2*15/Pblue), name="objgrad", dimnames=list(NULL,c("pred","pyellow","pgreen","pblue"))), mxFitFunctionAlgebra(algebra="fitfunc",gradient="objgrad",numObs=100), mxConstraint(Pred + Pyellow + Pgreen + Pblue - 1 == 0,name="indentifying") ) m1run <- mxRun(m1) summary(m1run) #Estimated proportions: p <- c(0.43,0.22,0.2,0.15) #Basis for the nullspace of the constraint Jacobian: U <- Null(t(m1run$output$constraintJacobian)) #Repeated-sampling covariance matrix: U%*%solve(t(U)%*%(m1run$output\$hessian/2)%*%U)%*%t(U)
#Analytic value of covariance matrix, for comparison:
(diag(p)-outer(p,p))/100
Offline
Joined: 07/31/2009 - 14:25
Did the Hessian %&% nullspace method get implemented?

Hi Rob,
Did computing the sampling covariance matrix from the usual Hessian and the nullspace of the constraint Jacobian matrix get implemented in the backend?

If so, can mxSE now be allowed to work with constraints? Would be very useful.

Offline
Joined: 01/24/2014 - 12:15
yes

Yes to both questions.

Offline
Joined: 07/31/2009 - 14:25
Great that mxSE works with constraints!

Fantastic!

cp1 = umxCP(selDVs = c('anx', 'dep', 'ocd'), sep = "", mzData = mz, dzData = dz)
mxSE(top.as_std, cp1)

#          [,1]      [,2]       [,3]
# [1,] 0.1458348 0.0000000 0.00000000
# [2,] 0.0000000 0.3616857 0.00000000
# [3,] 0.0000000 0.0000000 0.09989182
Offline
Joined: 01/24/2014 - 12:15
thanks

Thanks. The whole standard-errors-with-constraints feature was a few weeks' worth of effort on my part.