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

10 posts / 0 new
Offline
Joined: 02/14/2012 - 17:37
Finding likelihood-based confidence interval when a matrix contains definition variables
AttachmentSize
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:

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.

Offline
Joined: 07/31/2009 - 15:24
Whoops. This is a bug. It

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"))

Offline
Joined: 02/14/2012 - 17:37
Michael, thank you for very

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]).

Offline
Joined: 07/31/2009 - 15:24
Hmm,

Hmm, try:

mxAlgebra(expression = groupdiff / sqrt(S[7,7]), name = "alpha"),
mxCI(c("alpha"))
Offline
Joined: 02/14/2012 - 17:37
It actually works! I do not

It actually works! I do not know that the label can be used in the mxAlgebra. Thank you very much!

Offline
Joined: 03/29/2012 - 09:40
CI's for GxE

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:

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?

Offline
Joined: 07/31/2009 - 15:24
Err, I assume the algebra

Err, I assume the algebra looks like the following, given that Mod2 is a definition variable:

varA <- mxAlgebra(expression=(am + data.Mod2%x%aModm) %*% t(am + data.Mod2%x%aModm), name="Am")

With regards to getting CI's across the range of a definition variable, I'll need to check with the OpenMx team.

Offline
Joined: 03/29/2012 - 09:40
Well, yes, Mod2 is an OpenMx

Well, yes, Mod2 is an OpenMx label that I created before by

def2    <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels='data.sm', name="Mod2")

Would be great if there was a solution!

Offline
Joined: 07/31/2009 - 15:14
Dummy vector?

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:

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.

Offline
Joined: 03/29/2012 - 09:40
Thanks a lot! This is

Thanks a lot! This is fantastic! I even think I found a bug in the way I was trying to do myself!