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.

Attachment | Size |
---|---|

datawide.csv | 9.21 KB |

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!