# inst/models/passing/multinomSE.R In OpenMx/OpenMx: Extended Structural Equation Modelling

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

# #Custom compute plan, to calculate Hessian despite constraints:
# plan <- omxDefaultComputePlan()
# plan\$steps <- list(GD=plan\$steps\$GD)#, ND=plan\$steps\$ND, SE=plan\$steps\$SE, RD=plan\$steps\$RD, RE=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"),
dimnames=list(NULL,c("pred","pyellow","pgreen","pblue"))),
mxConstraint(Pred + Pyellow + Pgreen + Pblue - 1 == 0,name="indentifying")
)
m1run <- mxRun(m1)
summary(m1run)
omxCheckCloseEnough(coef(m1run),c(0.43,0.22,0.2,0.15),1e-7)
#Estimated proportions:
p <- round(coef(m1run),2)

#Basis for the nullspace of the constraint Jacobian:
U <- Null(t(m1run\$output\$constraintJacobian))
#Repeated-sampling covariance matrix:
( cov1 <- U%*%solve(t(U)%*%(m1run\$output\$hessian/2)%*%U)%*%t(U) )
omxCheckCloseEnough(cov1, m1run\$output\$vcov)
#Analytic value of covariance matrix, for comparison:
( cov2 <- (diag(p)-outer(p,p))/100 )
omxCheckCloseEnough(cov2, m1run\$output\$vcov, 1e-8)

#Compare to the "covariance matrix" you'd get directly from the Hessian (which is wrong):
solve(m1run\$output\$hessian/2)
#Compare to the "covariance matrix" you get from the optimizer's BFGS approximation to the Hessian (which is even wronger):
if(mxOption(NULL,"Default optimizer")=="NPSOL"){
print(solve( t(m1run\$output\$hessianCholesky) %*% m1run\$output\$hessianCholesky/2 ))
} else{
print(solve(m1run\$output\$LagrHessian/2))
}
```
OpenMx/OpenMx documentation built on Dec. 5, 2019, 4:22 a.m.