I'm trying to fit a very simple model using both the covariance model and the means model. I always get an error message concerning dimnames whether I specify the model via paths or using matrices. If I use matrices, it complains there are no dimnames for the F matrix - yet if I add dimnames for the F matrix it still complains.
All I want to be able to do is add a means model to the following OpenMx script that only has a covariance model:
> one <- mxModel("me",
+ type="RAM",
+ manifestVars = manifests,
+ latentVars = latents,
+ mxPath(from=latents,to=manifests,values=c(1,1,1),free=c(FALSE,TRUE,TRUE)),
+ mxPath(from=manifests,arrows=2,values=0.5,
+ free=TRUE,lbound=0),
+ mxPath(from=latents,arrows=2,values=0.5,free=TRUE,lbound=0),
+ mxAlgebra(sqrt(S),name="SD"),
+ mxData(observed=cov(OC),type="cov",numObs=42)
+ )
> summary(results <- mxRun(one))
Running me
uoaqtot uoeqtor aceqtot
Min. :1.372 Min. :1.366 Min. :1.441
1st Qu.:1.510 1st Qu.:1.369 1st Qu.:1.544
Median :1.648 Median :1.372 Median :1.648
Mean :1.619 Mean :1.393 Mean :2.146
3rd Qu.:1.742 3rd Qu.:1.407 3rd Qu.:2.498
Max. :1.836 Max. :1.441 Max. :3.349
name matrix row col Estimate Std.Error
1 A uoeqtor mu 0.8747084 0.07178385
2 A aceqtot mu 1.0500758 0.12941151
3 S uoaqtot uoaqtot 0.2671579 0.10559444
4 S uoeqtor uoeqtor 0.1660068 0.07742470
5 S aceqtot aceqtot 1.6189621 0.27089039
6 S mu mu 1.5689936 0.29003449
Observed statistics: 6
Estimated parameters: 6
Degrees of freedom: 0
-2 log likelihood: 126.5119
Saturated -2 log likelihood: 126.5119
Chi-Square: 9.023893e-12
p: 0
AIC (Mx): 9.023893e-12
BIC (Mx): 4.511946e-12
adjusted BIC:
RMSEA: Inf
Any insights would be greatly appreciated.
Could you post a copy of the script that generates the error?
From R, type
cov(OC)
Are there dimnames for the result of that calculation?
I suspect the script is doing something like:
Which don't work. But I can't know for sure until I see the script that throws the error.
I will post the code that uses the A, S, F, and M matrices but in the meantime, is there no way to add a means model to the code I showed? I want to be able to compute both intercepts (alpha's) and slopes (beta's). In addition, I would like to be able to compute the square roots of all the variance parameters and to comute the standard errors for ratios of the slopes. (I've done all this in the original Mx, I'm just having trouble converting over to OpenMx.) Also, please ignore the mxAlgebra function in the code. Thanks. I'm looking forward to coming up to speed using OpenMx - it's a great addition to R.
Oh, sorry I didn't read the question close enough. You should be able to specify a means path by adding a path where the 'to' argument is the reserved keyword "one". See http://openmx.psyc.virginia.edu/repoview/1/trunk/demo/TwoFactorModel_PathCov.R for an example of a means path.
Thanks. That helped. I show the code and some of the results below. The variance parameters (sigma1_2, etc.) and the beta's are nice, but what I really want are functions of these parameters, e.g., sqrt(sigma2_2)/beta_2. I have done this using Mx and Mx then gives confidence intervals on the results. I've tried adding an mxAlgebra, but OpenMx just seems to ignore it. Thanks. (The manual is great, but many of the examples are just too complex!)
> cov(OC)
uoaqtot uoeqtor aceqtot
uoaqtot 1.836152 1.372412 1.647562
uoeqtor 1.372412 1.366467 1.441137
aceqtot 1.647562 1.441137 3.349028
manifests <- names(OC)
latents <- "mu"
avg <- mean(OC)
grndavg <- mean(mean(OC))
one <- mxModel("me",
type="RAM",
manifestVars = manifests,
latentVars = latents,
mxPath(from=latents,to=manifests,values=c(1,1,1),free=c(FALSE,TRUE,TRUE),
labels=c("beta1","beta2","beta3")),
mxPath(from=manifests,arrows=2,values=0.5,
free=TRUE,lbound=0,labels=c("sigma1_2","sigma2_2","sigma3_2")),
mxPath(from=latents,arrows=2,values=0.5,free=TRUE,lbound=0,
labels="sigma2"),
mxPath(from="one",to=c(manifests,latents),arrows=1,free=c(T,T,T,F),
values=c(0,0,0,mean(mean(OC))),labels=c("alpha1","alpha2","alpha3","mubar")),
mxData(observed=cov(OC),type="cov",means=mean(OC),numObs=42)
)
summary(results <- mxRun(one))
1 beta2 A uoeqtor mu 0.8747087 0.07431760
2 beta3 A aceqtot mu 1.0500760 0.13397846
3 sigma1_2 S uoaqtot uoaqtot 0.2671579 0.10684218
4 sigma2_2 S uoeqtor uoeqtor 0.1660069 0.07848701
5 sigma3_2 S aceqtot aceqtot 1.6189629 0.29032748
6 sigma2 S mu mu 1.5689927 0.30605153
7 alpha1 M 1 uoaqtot 0.4636648 0.15305310
8 alpha2 M 1 uoeqtor -0.4365962 0.30146336
9 alpha3 M 1 aceqtot 0.2466312 0.52584192
I read the question with two possible interpretations, so I'm going to answer both interpretations, and I hope one of the answers is helpful.
To evaluate some expression AFTER the model has finished execution, you can use the mxEval() command. For example, you can run the following after executing the model:
mxEval(sqrt(sigma2_2)/beta2, model)
To evaluate an algebraic expression DURING the model execution, use the mxAlgebra() command. The following algebra can be added to your model:
mxAlgebra(sqrt(sigma2_2)/beta2, name = 'sqrtAlgebra')
Sadly we are not finished implementing confidence intervals. The feature is on the short list of future work. In the meantime, it is possible to implement confidence intervals in OpenMx using MxAlgebra expressions. See the following thread: http://openmx.psyc.virginia.edu/thread/153#comment-1054
Thanks but can I get the standard error? It doesn't show up in the summary.
> results$sqrtAlgebra
mxAlgebra 'sqrtAlgebra'
@formula: sqrt(sigma2_2)/beta2
@result:
[,1]
[1,] 0.4658001
dimnames: NULL
Ah sorry. Yes, standard error estimates currently show up only for the free parameters. I understand conceptually how to generate standard errors for the free parameters. The optimizer walks around the local minimum and then does some kind of math. Perhaps someone with more knowledge can comment on standard error estimates for arbitrary algebras.
There is no easy solution for arbitrary algebras because they are... arbitrary. You could put anything into the algebra including a variety of nonlinear functions.
However, likelihood based confidence intervals should be possible. Likelihood based confidence intervals are near the top of the feature queue that we have yet to implement. I expect we will have them by summer.
I'm assuming that because Mx produced confidence intervals for all the parameters and functions of parameters, then it should also be possible in OpenMx.
"Classic Mx" has likelihood based confidence intervals available. However, there is still some discussion about how these should be implemented in OpenMx so that they are not misleading. Note that just because confidence intervals are reported does not save you from thinking carefully about whether you really should be using them or not.
Perhaps Mike Neale could reply to this for a perspective on how classic Mx handles this problem?
Likelihood-based confidence intervals have some nice properties. One is that they can be obtained for arbitrarily complex functions of parameters, as was implemented in Classic Mx. A second is that inferences based on the interval of parameter a would be the same as inferences based on the interval of a^2. For example, whether zero is included in the 95% CI one would get the same answer for both a and a^2. That is not the case for inferences based on standard errors computed from the Hessian (the covariance matrix of the parameter estimates). See Neale M.C., Miller MB. (1997) The use of likelihood-based confidence intervals in genetic models. Behavior Genetics 27:113-120 (http://www.vipbg.vcu.edu/vipbg/Articles/behavgen-use-1997.pdf) for an explanation of this property.
There are of course disadvantages to the method as well. One is that two additional optimization runs are required to obtain the upper and lower CI's, which can cost CPU time. A second, corollary problem is that optimization is not an exact science, so the optimization to find intervals may not be successful. Classic Mx would deal with this by i) automatically refitting up to five times in cases where IFAIL=6 was found; and ii) reporting error codes for optimization in the output, including a special code (8) in cases where the estimate of the lower or upper CI was the same as the MLE. Similar flags, and perhaps less terse explanation of potential problems with the reported CI's, should, I think, be part of the OpenMx CI output.
Note that until there is an OpenMx implemention of an mxInterval() function, it is possible to obtain CI's "manually" by revising the fit function, as in the thread http://openmx.psyc.virginia.edu/thread/153 (see post #9). This should work for functions of parameters also, perhaps using the [] notation in order to identify a particular element of a matrix.