Confidence intervals (granted)
Posted on

Forums
Confidence intervals are basic to most results sections now, so implementing Mx's 'Intervals' likelihood-based confidence intervals (Neale & Miller, 1997) (plus standard error-based CIs for speed) would be high on the wish-list
maybe something like
mxInterval(..., percent=95, type="ML") where ... consists of matrix addresses
mxInterval(z.A) would mean all cells of A
mxInterval(z.A[1:2,], type="SE") first two rows of A, standard-error based method
Agreed, and I like the syntax
Log in or register to post comments
From the standpoint of the
We also would need to provide pathic versions of the interface.
Log in or register to post comments
In reply to From the standpoint of the by Steve
I guess the CI matrices would
In that case, instead of a type slot, might be best allow both sets of CI matrices to be added: that way people could say
mxCI(z.A, "SE");
mxCI(z.A, "ML");
and get both kinds to compare etc.
Log in or register to post comments
If I'm reading Tim's
Log in or register to post comments
In reply to If I'm reading Tim's by Ryne
Hi Ryne, I thought CIs will
I thought CIs will have to be part of mxModel, as asking for confidence intervals causes the optimiser to fit the model, then move each (requested) parameter estimate down (and up), holding the others constant, until fit declines by some number of units related to the confidence requested. But perhaps it could be handled in a separate mxEval-like mxCI() function, which would take the model and current parameters as input?
Log in or register to post comments
In reply to Hi Ryne, I thought CIs will by tbates
Right you are. It's going to
Log in or register to post comments
In reply to Right you are. It's going to by Steve
I would say that if a model
The code for choosing new values and homing in on the point where the fit degrades below the desired CI must be in old mx for something to work off.
Log in or register to post comments
In reply to I would say that if a model by tbates
Hey guys, I just discovered
But is it a goal in OpenMx to keep names that probably come from MX syntax like mxInterval?
If no, you could maybe consider R generic functions intuitively known to users such as confint()?
Log in or register to post comments
In reply to I would say that if a model by tbates
Yes, the Mx1 code is there.
Possibly, these revised fit functions won't be too bad to implement with an Algebra type fit function.
For example, using the UnivariateTwinAnalysis_MatrixRaw.R script from the demo directory, we can obtain the likelihood-based CI's as follows:
CIaupper<-mxModel(twinACEFit,mxMatrix("Full",values=mxEval(objective,twinACEFit),name="oldfit"),mxAlgebra(((oldfit+3.84)-(MZ.objective+DZ.objective))^2-A,name="upperCIa"),mxAlgebraObjective("upperCIa"))
CIalower<-mxModel(twinACEFit,mxMatrix("Full",values=mxEval(objective,twinACEFit),name="oldfit"),mxAlgebra(((oldfit+3.84)-(MZ.objective+DZ.objective))^2+A,name="lowerCIa"),mxAlgebraObjective("lowerCIa"))
then running
runCIalower<-mxRun(CIalower)
runCIalower@algebras$A
runCIaupper<-mxRun(CIaupper)
runCIaupper@algebras$A
which are the lower and upper CI's on the computed matrix element A[1,1]
Obviously, making this an automated function would be much nicer for the user! But it's pretty neat that it can be done without any programing on behalf of the OpenMx team.
Log in or register to post comments
In reply to Yes, the Mx1 code is there. by neale
Thanks mike: those are
I think it would still be very well worth while to support this core functionality inside an openMx-package supported mxInterval() function, able to be added to an mxModel, calling itself recursively, or as a parameter of mxRun
t
Log in or register to post comments
In reply to Thanks mike: those are by tbates
Just a thought too: The
Would be good to have the ability to specify a range of cells to get CIs on.
and to be able to specify CIs on RAM mxPaths by c(name,)
Log in or register to post comments
In reply to Just a thought too: The by tbates
Agreed - syntax of the type
Second, on the issue of computation time, you are right, likelihood-based CI's are expensive. They can be done in parallel though, so perhaps this is another area where Drs. Wilde & Kenny may be able to help. With a 48 core laptop on the horizon, it's tempting indeed to make this example run locally in parallel.
Log in or register to post comments
In reply to Agreed - syntax of the type by neale
With the zero influence I
-Scott
Log in or register to post comments
In reply to With the zero influence I by svrieze
Ah, that would be useful. In
model <- mxModel(....)
modelOut <- mxRun(model)
library(snowfall)
sfInit(parallel = TRUE, cpus = ???)
sfLibrary(OpenMx)
intervals1 <- mxModel(modelOut, mxCI('a'), name = "intervals1", independent = TRUE)
intervals2 <- mxModel(modelOut, mxCI('b'), name = "intervals2", independent = TRUE)
intervals3 <- mxModel(modelOut, mxCI('c'), name = "intervals3", independent = TRUE)
intervals4 <- mxModel(modelOut, mxCI('d'), name = "intervals4", independent = TRUE)
intervals5 <- mxModel(modelOut, mxCI('e'), name = "intervals5", independent = TRUE)
topModel <- mxModel('topModel', intervals1, intervals2, intervals3, intervals4, intervals5)
results <- mxRun(topModel, intervals = TRUE)
sfStop()
UPDATE: I forgot to give the submodels new names. The name collision would have resulted in only one submodel stored as a child of the topModel.
Log in or register to post comments
In reply to Ah, that would be useful. In by mspiegel
That's a great idea, thanks!
The renaming of submodels is a problem, to be sure. My ACE common pathway model has several submodels ('ACE', 'DZ','MZ'), and renaming them for each confidence interval is time-consuming. It seems some kind of built-in parallelizing confidence interval function would save a lot of people a lot of time. As others are saying, reporting confidence intervals is (and should) becoming routine. Each submodel could be automatically renamed to make it easier on the user (e.g., for the first confidence interval requested it could have 'ACE1', 'DZ1', 'MZ1', and so on).
Just thinking out loud. THanks for an extremely useful product thusfar!
-Scott
Log in or register to post comments
In reply to With the zero influence I by svrieze
Ah yes the model naming
So when I tested it using
foo1 <- modelNameAddPrefix(twinACEModel, "foo1_")
the resulting model is named 'foo1_twinACE' and has three submodels 'foo1_common', 'foo1_MZ', 'foo1_DZ'. It's not exactly user-friendly at the moment, but it will help in running the confidence intervals.Oh, plus I found a bug in the function mxRename(). It will currently not rename the appropriate confidence intervals, because confidence intervals had not been implemented when mxRename() was written. But you should be specifying the confidence intervals manually after creating all the submodels for this particular workflow.
UPDATE: Oops, I forgot that child models could refer to entities in the parent models. I had to change the recursion step to apply to the entire nested model. The code has been updated.
Log in or register to post comments
In reply to Just a thought too: The by tbates
ML CI's are inherently
require(boot)
mles<-function(dataset,wt){
manifests <- names(dataset)
latents <- c("G")
covwt<-cov.wt(dataset,wt)
mlevals<-mxRun(mxModel("One Factor", type="RAM",
manifestVars = manifests,
latentVars = latents,
mxPath(from=latents, to=manifests),
mxPath(from=manifests, arrows=2),
mxPath(from=latents, arrows=2, free=F, values=1.0),
mxData(covwt$cov, type="cov",
numObs=500)))
return(as.vector(mlevals@output$estimate))}
boot(demoOneFactor,mles,R=1000)
Possibly a wrapper function would facilitate joint use of the boot and OpenMx packages. Note that the mlevals in the above example might be augmented to have the results of algebras included in the return statement.
Log in or register to post comments
In reply to ML CI's are inherently by neale
could someone please post a
thanks
Log in or register to post comments
In reply to could someone please post a by smedland
Same here. Having svn upped
> require(OpenMx)
> data(demoOneFactor)
> factorModel <- mxModel("One Factor",
+ mxMatrix("Full", 5, 1, values=0.2,free=T, name="A"),
+ mxMatrix("Symm", 1, 1, values=1,free=F, name="L"),
+ mxMatrix("Diag", 5, 5, values=1,free=T, name="U"),
+ mxAlgebra(A %*% L %*% t(A) + U, name="R"),
+ mxCI(c('A[1,1]','R[1,1]')),
+ mxMLObjective("R", dimnames = names(demoOneFactor)),
+ mxData(cov(demoOneFactor), type="cov", numObs=500))
> factorModelOuput <- (mxRun(factorModel,intervals=T))
Running One Factor
Warning messages:
1: In model 'One Factor' while computing the lower bound for 'One Factor.A[1,1]' NPSOL returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)
2: In model 'One Factor' while computing the lower bound for 'One Factor.R[1,1]' NPSOL returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)
3: In model 'One Factor' while computing the upper bound for 'One Factor.A[1,1]' NPSOL returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)
4: In model 'One Factor' while computing the upper bound for 'One Factor.R[1,1]' NPSOL returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)
> factorModelSummary <- summary(factorModelOutput)
> factorModelOutput@output$confidenceIntervals
[,1] [,2]
>
And that's it. Summary is no help either.
Log in or register to post comments
In reply to Same here. Having svn upped by neale
Umm, I'm getting confidence
R 2.11.0 and OpenMx subversion revision 1297. I'll logon to some other machines and try running it on older versions of R. Works on R 2.9.2 and 2.10.1 also.
Log in or register to post comments
In reply to Umm, I'm getting confidence by mspiegel
I don't think it's an R
?OpenMx indicates Package OpenMx version 0.3.3-1264
but I suspect that the version number is only updated with binary releases. It coincides with what I get during the build:
tar: trunk/build/OpenMx_0.3.3-1264.tar: Can't add archive to itself
and I'm pretty sure the building is going according to plan:
Michael-Neales-Macbook-Pro-4% svn up
At revision 1297.
Michael-Neales-Macbook-Pro-4% cd trunk
Michael-Neales-Macbook-Pro-4% sudo make install
I retried this procedure after quitting R, then fired up R again and still no CI's in the summary:
> require(OpenMx)
> data(demoOneFactor)
> factorModel <- mxModel("One Factor",
+ mxMatrix("Full", 5, 1, values=0.2,free=T, name="A"),
+ mxMatrix("Symm", 1, 1, values=1,free=F, name="L"),
+ mxMatrix("Diag", 5, 5, values=1,free=T, name="U"),
+ mxAlgebra(A %*% L %*% t(A) + U, name="R"),
+ mxCI(c('A[1,1]','R[1,1]')),
+ mxMLObjective("R", dimnames = names(demoOneFactor)),
+ mxData(cov(demoOneFactor), type="cov", numObs=500))
> factorModelOuput <- (mxRun(factorModel,intervals=T))
> factorModelSummary <- summary(factorModelOutput)
> factorModelSummary
data:
$`One Factor.data`
$`One Factor.data`$cov
x1 x2 x3 x4 x5
x1 0.1985443 0.1999953 0.2311884 0.2783865 0.3155943
x2 0.1999953 0.2916950 0.2924566 0.3515298 0.4019234
x3 0.2311884 0.2924566 0.3740354 0.4061291 0.4573587
x4 0.2783865 0.3515298 0.4061291 0.5332788 0.5610769
x5 0.3155943 0.4019234 0.4573587 0.5610769 0.6703023
free parameters: A x1 G 0.39715212 0.015549769 A x2 G 0.50366111 0.018232514 A x3 G 0.57724141 0.020448402 A x4 G 0.70277369 0.024011418 A x5 G 0.79624998 0.026669452 S x1 x1 0.04081419 0.002812717 S x2 x2 0.03801999 0.002805794 S x3 x3 0.04082718 0.003152308 S x4 x4 0.03938706 0.003408875 S x5 x5 0.03628712 0.003678561
name matrix row col Estimate Std.Error
1
2
3
4
5
6
7
8
9
10
observed statistics: 15
estimated parameters: 10
degrees of freedom: 5
-2 log likelihood: -3648.281
saturated -2 log likelihood: -3655.665
number of observations: 500
chi-square: 7.384002
p: 0.1936117
AIC (Mx): -2.615998
BIC (Mx): -11.84452
adjusted BIC:
RMSEA: 0.03088043
timestamp: 2010-07-01 10:24:30
frontend time: 0.08384323 secs
backend time: 0.01688194 secs
independent submodels time: 5.698204e-05 secs
wall clock time: 0.1007822 secs
cpu time: 0.1007822 secs
openmx version number: 0.3.3-1264
Weird, huh?
Log in or register to post comments
In reply to I don't think it's an R by neale
Same errors as Mike. I'm on
svn info tells me: Revision: 1298
Log in or register to post comments
In reply to Same errors as Mike. I'm on by tbates
Er, those seem like different
rm(list = ls())
DOES help.
Now I get the errors reported by Tim! But I am able to see the CI's if I don't go via the summary object. So basically the CI's are causing the creation of a summary object to fail, and my earlier confusion was caused by my inspection of a historical summary object of the same name (created by an earlier version of OpenMx which didn't have this bug).
Log in or register to post comments
In reply to Er, those seem like different by neale
Yes, having something fail
Log in or register to post comments
In reply to Yes, having something fail by tbates
Copy/pasting Mike's input
factorModelOuput <- (mxRun(factorModel,intervals=T))
. The "t" is missing in Output.Log in or register to post comments
In reply to Copy/pasting Mike's input by mspiegel
Oops! Sorry. Much ado about
Log in or register to post comments
In reply to ML CI's are inherently by neale
CIs by bootstrapping
Log in or register to post comments
How did you end up solving
Log in or register to post comments
Can I check that I understand
Based on the code Mike wrote above, If I have 4 variables and a best fit AE model, then by specifying A for each variable (a11*a11, a22*a22, etc), I should be able to obtain corresponding CI's, as in;
A <- mxEval(a11*a11, multiCholAEFit)
CIaupper<-mxModel(multiCholAEFit,mxMatrix("Full",values=mxEval(objective, multiCholAEFit),name="oldfit"),mxAlgebra(((oldfit+3.84)-(MZ.objective+DZ.objective+Sib.objective))^2-A,name="upperCIa"),mxAlgebraObjective("upperCIa"))
CIalower<-mxModel(multiCholAEFit,mxMatrix("Full",values=mxEval(objective, multiCholAEFit),name="oldfit"),mxAlgebra(((oldfit+3.84)-(MZ.objective+DZ.objective+Sib.objective))^2+A,name="lowerCIa"),mxAlgebraObjective("lowerCIa"))
runCIalower<-mxRun(CIalower)
runCIalower@algebras$A
runCIaupper<-mxRun(CIaupper)
runCIaupper@algebras$A
If this is right, like Cathy and Marcel's recent post, I keep running into code reds (even when re-running an optimised fit).
Log in or register to post comments
In reply to Can I check that I understand by pgseye
I wouldn't use mxEval() for
Also, please be aware that a better interface for requesting & estimating confidence intervals is due to be released shortly.
Log in or register to post comments
In reply to I wouldn't use mxEval() for by neale
That's correct. mxEval()
And yes, a confidence interval interface will be released in OpenMx 1.0.
Log in or register to post comments