Finding likelihood-based confidence interval when a matrix contains definition variables

Attachment | Size |
---|---|
datawide.csv | 9.21 KB |
I have a problem running multilevel model using a SEM framework to get a likelihood-based CI of the effect size. I have two factors representing (random) intercept and (fixed) slope. I attach my data file here. Use the following code to get the dataset:
datawide <- read.csv("datawide.csv")
The dataset has been transformed into a wide format. Y is the response variable. X is the upper-level predictor. Z is the lower-level predictor without any centering. My target model is
L1: Y = B_0 + B_1 * Z + r
L2: B_0 = X + u
I would like to find a likelihood-based CI of B_1 / sqrt(var(r) + (B_1^2 * var(Z)))
Here is my OpenMx script:
ylab <- paste("y", 1:5, sep="")
zlab <- paste("z", 1:5, sep="")
colnames(datawide) <- c("id", "x", ylab, zlab)
latentLab <- c("intcept", "slope")
Alab <- matrix(NA, 8, 8)
Alab[7, 6] <- "groupdiff"
Alab[1:5, 8] <- paste("data.", zlab, sep="")
Aval <- matrix(0, 8, 8)
Aval[7, 6] <- 0.1 # Need parameter values
Aval[1:5, 7] <- 1
Aval[1:5, 8] <- 1
Afree <- matrix(FALSE, 8, 8)
Afree[7, 6] <- TRUE
colnames(Alab) <- colnames(Aval) <- colnames(Afree) <- rownames(Alab) <- rownames(Aval) <- rownames(Afree) <- c(ylab, "x", latentLab)
Slab <- matrix(NA, 8, 8)
diag(Slab) <- "l1error"
Slab[6, 6] <- "varX"
Slab[7, 7] <- "l2error"
Slab[8, 8] <- NA
Sval <- matrix(0, 8, 8)
diag(Sval) <- 0.8 # Need parameter values
Sval[6, 6] <- 0.25 # Need parameter values
Sval[7, 7] <- 0.2 # Need parameter values
Sval[8, 8] <- 0
Sfree <- matrix(FALSE, 8, 8)
diag(Sfree) <- TRUE
Sfree[8, 8] <- FALSE
colnames(Slab) <- colnames(Sval) <- colnames(Sfree) <- rownames(Slab) <- rownames(Sval) <- rownames(Sfree) <- c(ylab, "x", latentLab)
Fval <- cbind(diag(6), matrix(0, 6, 2))
Flab <- matrix(NA, 6, 8)
Ffree <- matrix(FALSE, 6, 8)
colnames(Flab) <- colnames(Fval) <- colnames(Ffree) <- c(ylab, "x", latentLab)
rownames(Flab) <- rownames(Fval) <- rownames(Ffree) <- c(ylab, "x")
Mlab <- matrix(c(rep(NA, 5), "meanX", "mean0", "betaW"), 1, 8)
Mval <- matrix(c(rep(0, 5), 0.5, 0.2, 0.8), 1, 8) # Need parameter values
Mfree <- matrix(c(rep(FALSE, 5), rep(TRUE, 3)), 1, 8)
colnames(Mlab) <- colnames(Mval) <- colnames(Mfree) <- c(ylab, "x", latentLab)
onecov <- mxModel("With covariate model", type="RAM",
mxData(datawide[,c(ylab, "x", zlab)], type="raw"),
mxMatrix(type="Full", nrow=8, ncol=8, values=Aval, free=Afree, labels=Alab, name="A"),
mxMatrix(type="Symm", nrow=8, ncol=8, values=Sval, free=Sfree, labels=Slab, name="S"),
mxMatrix(type="Full", nrow=6, ncol=8, values=Fval, free=Ffree, labels=Flab, name="F"),
mxMatrix(type="Full", nrow=1, ncol=8, values=Mval, free=Mfree, labels=Mlab, name="M"),
mxRAMObjective("A","S","F","M", dimnames=c(ylab, "x", latentLab)),
mxMatrix(type="Full", nrow=8, ncol=1, values=c(rep(0, 6), 1, 0), free=FALSE, name="s7"),
mxMatrix(type="Full", nrow=8, ncol=1, values=c(rep(0, 5), 1, 0, 0), free=FALSE, name="s6")
) # - end model
onecovfit <- mxRun(onecov, intervals=TRUE)
summary(onecovfit)
This code works well. Next, I added two objects into the mxModel command
mxAlgebra(expression = A[7, 6], name = "alpha"),
mxCI(c("alpha"))
When I added these lines, the confidence interval can not be computed. Here is the error message:
Error: Could not find the dataset 'data' in model 'With covariate model' used by the definition variable 'data.z1'
However, when I ran the following codes, it works fine:
mxAlgebra(expression = S[7, 7], name = "alpha"),
mxCI(c("alpha"))
Do you have any suggestions on how to find the confidence interval that involves A[7, 6] which does not really involve with the definition variable? Thank you very much in advance.
Whoops. This is a bug. It
mxCI(c("groupdiff"))
Log in or register to post comments
In reply to Whoops. This is a bug. It by mspiegel
Michael, thank you for very
Log in or register to post comments
In reply to Michael, thank you for very by psunthud
Hmm,
mxAlgebra(expression = groupdiff / sqrt(S[7,7]), name = "alpha"),
mxCI(c("alpha"))
Log in or register to post comments
In reply to Hmm, by mspiegel
It actually works! I do not
Log in or register to post comments
In reply to Whoops. This is a bug. It by mspiegel
CI's for GxE
However, I have a question about estimating CI when a definition variable has multiple values (not just one like in the above example). I have a GxE model and would like to get a CI for, let's say, genetic variance having different values for a moderator. At the moment i specified it like this:
varA = mxAlgebra(expression=(am + Mod2%x%aModm) %*% t(am + Mod2%x%aModm), name="Am")
ciVm = mxCI(c('Am'))
where
Mod2
is a definition variable. However, as a result I get only one set of values for the CI calculated as I assume for the last entry in the dataset. Is it possible to get CI's for Am for the whole range of the definition variable values?Thank you in advance!
Log in or register to post comments
In reply to CI's for GxE by Julia
Err, I assume the algebra
Mod2
is a definition variable:With regards to getting CI's across the range of a definition variable, I'll need to check with the OpenMx team.
Log in or register to post comments
In reply to Err, I assume the algebra by mspiegel
Well, yes, Mod2 is an OpenMx
Mod2
is an OpenMx label that I created before bydef2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels='data.sm', name="Mod2")
Would be great if there was a solution!
Log in or register to post comments
In reply to Well, yes, Mod2 is an OpenMx by Julia
Dummy vector?
varA <- mxAlgebra(expression=(am + data.Mod2%x%aModm) %*% t(am + data.Mod2%x%aModm), name="Am")
you could bung in a dummy vector of values and tweak the algebra a little to give you a vector of output values
myVectorOfDefVarValues <- mxMatrix("Full",nValues,1,values=c(.00, .05, .10, .15, .20),name="myVec"),
myUnitVector <- mxMatrix("Unit",nValues,1,name="myUnit"),
myVarAValues <- mxAlgebra(expression=(myUnit %x% am + myVec %x% aModm) *
(myUnit %x% am + myVec %x% aModm), name="myVarAValues"),
myCIs <- mxCI("myVarAValues")
sound reasonable? Note the replacement of %*% by just * and the lack of transpose... I am assuming that this is univariate and that am is therefore 1x1. Some tweakification required (say using vec()) otherwise.
Log in or register to post comments
In reply to Dummy vector? by neale
Thanks a lot! This is
Log in or register to post comments