Attachment | Size |
---|---|
datawide.csv | 9.21 KB |
Hi everyone,
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 will be fixed in OpenMx 1.2.5. As a workaround, you can use the label of the free parameter to specify the confidence interval. Use something like this:
mxCI(c("groupdiff"))
Michael, thank you for very quick response. I will wait for the next version. I need to use the "groupdiff" for further computation. That is A[7, 6] / sqrt(S[7, 7]).
Hmm, try:
It actually works! I do not know that the label can be used in the mxAlgebra. Thank you very much!
Thank you so much for fixing this bug! I was struggling with this problem for a day and then found your post here!
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:
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!
Err, I assume the algebra looks like the following, given that
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.
Well, yes,
Mod2
is an OpenMx label that I created before byWould be great if there was a solution!
I am thinking that you would like to plot say genetic variance across a range of moderator values, and also the CI of genetic variance at each point. It seems to me that you could manually evaluate CI's for a suitable set - it would be computationally intensive to do this for every definition variable value (and no point if there were duplicates). So where you have:
you could bung in a dummy vector of values and tweak the algebra a little to give you a vector of output values
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.
Thanks a lot! This is fantastic! I even think I found a bug in the way I was trying to do myself!