small/NaN standard errors
I'm fitting a univariate threshold model on binary data and it appears to work well (see below for the script), except that in the summary I get either very small or NaN standard errors, like this:
name matrix row col Estimate Std.Error
1 ab11 ACE.aboys 1 1 0.5690466 2.121996e-314
2 cb11 ACE.cboys 1 1 0.7157107 6.365987e-314
3 eb11 ACE.eboys 1 1 0.4049003 NaN
4 ag11 ACE.agirls 1 1 0.2620615 NaN
5 cg11 ACE.cgirls 1 1 0.8921872 NaN
6 eg11 ACE.egirls 1 1 0.3678666 1.060998e-313
7 threall ACE.Tall 1 1 -0.5358325 NaN
8 thredzm ACE.Tdzm 1 1 -0.7374305 NaN
Another thing is that when I run the submodel which equates a, c and e across sex I get the following output in OpenMx:
name matrix row col Estimate Std.Error
1 ab11 ACE.aboys 1 1 0.4833659 NaN
2 cb11 ACE.cboys 1 1 0.7886090 NaN
3 eb11 ACE.eboys 1 1 0.3800700 NaN
4 threall ACE.Tall 1 1 -0.5369859 NaN
5 thredzm ACE.Tdzm 1 1 -0.7300285 NaN
Observed statistics: 2807
Estimated parameters: 5
Degrees of freedom: 2802
-2 log likelihood: 2961.374
Whereas when I run the same model in classic Mx it says:
X 1 1 1 95.0 0.4816 0.2191 0.6461 0 1 1 1
Y 1 1 1 95.0 0.7899 0.6788 0.8761 0 1 0 1
Z 1 1 1 95.0 0.3796 0.3077 0.4583 0 1 0 0
B 1 1 1 95.0 -0.5125 -0.5789 -0.4464 0 0 0 0
R 1 1 1 95.0 -0.6879 -0.7839 -0.5934 0 0 0 1
-2 times log-likelihood of data >>> 2954.333
Degrees of freedom >>>>>>>>>>>>>>>> 2802
So there is some difference in both the parameter estimates (mostly the thresholds) and in the -2LL, should this worry me? And I'm curious about the OpenMx standard errors, why they are so small or NaN.
Thanks again!
Lot
ACE model specification (OpenMx):
univACEOrdModel <- mxModel("univACEOrd",
mxModel("ACE",
# Matrices aboys, cboys, and eboys to store a, c, and e path coefficients for boys
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="ab11", name="aboys" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="cb11", name="cboys" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="eb11", name="eboys" ),
# Matrices Aboys, Cboys, and Eboys compute variance components for boys
mxAlgebra( expression=aboys %*% t(aboys), name="Aboys" ),
mxAlgebra( expression=cboys %*% t(cboys), name="Cboys" ),
mxAlgebra( expression=eboys %*% t(eboys), name="Eboys" ),
# Matrices agirls, cgirls, and egirls to store a, c, and e path coefficients for girls
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="ag11", name="agirls" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="cg11", name="cgirls" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="eg11", name="egirls" ),
# Matrices Agirls, Cgirls, and Egirls compute variance components for girls
mxAlgebra( expression=agirls %*% t(agirls), name="Agirls" ),
mxAlgebra( expression=cgirls %*% t(cgirls), name="Cgirls" ),
mxAlgebra( expression=egirls %*% t(egirls), name="Egirls" ),
# Algebra to compute total variance for boys and girls and standard deviations (diagonal only)
mxAlgebra( expression=Aboys+Cboys+Eboys, name="Vboys" ),
mxAlgebra( expression=Agirls+Cgirls+Egirls, name="Vgirls" ),
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( expression=solve(sqrt(I*Vboys)), name="sdboys"),
mxAlgebra( expression=solve(sqrt(I*Vgirls)), name="sdgirls"),
# Constraint on variance of ordinal variables
mxConstraint( alg1="Vboys", "=", alg2="I", name="Varboys"), # varianties sommeren tot 1 (identiteitsmatrix)
mxConstraint( alg1="Vgirls", "=", alg2="I", name="Vargirls"),
# Matrix & Algebra for expected means vector
mxMatrix( type="Zero", nrow=1, ncol=nv, name="M" ),
mxAlgebra( expression= cbind(M,M), name="expMean" ),
# threshold voor de hele groep:
mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=-.6, label="threall", name="Tall" ),
mxAlgebra( expression= cbind(Tall,Tall), dimnames=list('th1',selVars), name="expThreall" ),
# threshold voor de dz jongens:
mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=-.8, label="thredzm", name="Tdzm" ),
mxAlgebra( expression= cbind(Tdzm,Tdzm), dimnames=list('th1',selVars), name="expThredzm" ),
# Algebra for expected variance/covariance matrix in mzm
mxAlgebra( expression= rbind ( cbind(Aboys+Cboys+Eboys , Aboys+Cboys),
cbind(Aboys+Cboys , Aboys+Cboys+Eboys)), name="expCovmzm" ),
# Algebra for expected variance/covariance matrix in dzm, note use of 0.5, converted to 1*1 matrix
mxAlgebra( expression= rbind ( cbind(Aboys+Cboys+Eboys , 0.5%x%Aboys+Cboys),
cbind(0.5%x%Aboys+Cboys , Aboys+Cboys+Eboys)), name="expCovdzm" ),
# Algebra for expected variance/covariance matrix in mzf
mxAlgebra( expression= rbind ( cbind(Agirls+Cgirls+Egirls , Agirls+Cgirls),
cbind(Agirls+Cgirls , Agirls+Cgirls+Egirls)), name="expCovmzf" ),
# Algebra for expected variance/covariance matrix in dzf
mxAlgebra( expression= rbind ( cbind(Agirls+Cgirls+Egirls , .5%x%Agirls+Cgirls),
cbind(.5%x%Agirls+Cgirls , Agirls+Cgirls+Egirls)), name="expCovdzf" ),
# Algebra for expected variance/covariance matrix in dosmf
mxAlgebra( expression= rbind ( cbind(Aboys+Cboys+Eboys , .5%x%(aboys%*%t(agirls))+ cboys%*%t(cgirls)),
cbind(.5%x%(agirls%*%t(aboys))+ cgirls%*%t(cboys) , Agirls+Cgirls+Egirls)), name="expCovdosmf" )
),
mxModel("mzm",
mxData( observed=mzmdata, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovmzm", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThreall" )
),
mxModel("dzm",
mxData( observed=dzmdata, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovdzm", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThredzm" )
),
mxModel("mzf",
mxData( observed=mzfdata, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovmzf", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThreall" )
),
mxModel("dzf",
mxData( observed=dzfdata, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovdzf", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThreall" )
),
mxModel("dosmf",
mxData( observed=dosmfdata, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovdosmf", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThreall" )
),
mxAlgebra( expression=mzm.objective + dzm.objective + mzf.objective + dzf.objective + dosmf.objective, name="min2sumll" ),
mxAlgebraObjective("min2sumll")
)
univACEOrdFit <- mxRun(univACEOrdModel)
univACEOrdSumm <- summary(univACEOrdFit)
univACEOrdSumm
LL_ACE <- mxEval(objective, univACEOrdFit);LL_ACE
And the submodel:
ACEequalModel = univACEOrdModel
ACEequalModel$ACE$agirls <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="ab11", name="agirls" )
ACEequalModel$ACE$cgirls <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="cb11", name="cgirls" )
ACEequalModel$ACE$egirls <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="eb11", name="egirls" )
ACEequalFit = mxRun(ACEequalModel)
ACEequalSumm = summary(ACEequalFit)
ACEequalSumm
LL_ACEeq = mxEval(objective,ACEequalFit)
LL_ACEeq
Hi, I have the same problems,
Log in or register to post comments
I suspect that the NaN
When ordinal data are being analyzed, it is possible to circumvent this problem by estimating the variances and means, while constraining the first two thresholds to equal 0 and 1. This reparameterization avoids non-linear constraints, but it does become necessary to re-standardize the estimates. Restandardization does of course mess up the standard errors.
For now, perhaps the best thing to do would be to use the hack for likelihood-based confidence intervals, described in this post:
http://openmx.psyc.virginia.edu/thread/153#comment-1054
Log in or register to post comments
In reply to I suspect that the NaN by neale
Ah thanks very much! About
About the differences in estimates en likelihood, are they maybe because of the FIML estimator? The differences seem a bit large...
Lot
Log in or register to post comments
In reply to Ah thanks very much! About by lotgeels
Hello Lot Hmmm I agree that
Hmmm I agree that the differences are a bit large. Did both scripts use the same starting values? Both classic Mx and OpenMx have the same FIML estimator, which has been shown to agree much more closely in other examples.
In a case such as this, one would want to look at the script that ends up with the larger (more positive) -2lnL, i.e. the one that has the lower likelihood. Probably its starting values are not as good as those of the classic Mx script which found the solution more directly. Possibly, there is a model specification error in one or other of the scripts. One way to test this would be to fix the starting values in the OpenMx script so that they equal the *solution* from the Classic Mx script. If that does not yield the same (or very close, say <.1 difference) -2lnL, then the model is probably not specified the same in both scripts.
Log in or register to post comments