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

9 replies [Last post]
psunthud's picture
Joined: 02/14/2012

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)

This code works well. Next, I added two objects into the mxModel command

mxAlgebra(expression = A[7, 6], name = "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"),

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.

datawide.csv9.21 KB
mspiegel's picture
Joined: 07/31/2009
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"))

Julia's picture
Joined: 03/29/2012
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?

Thank you in advance!

mspiegel's picture
Joined: 07/31/2009
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.

Julia's picture
Joined: 03/29/2012
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='', name="Mod2") 

Would be great if there was a solution!

neale's picture
Joined: 07/31/2009
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.

Julia's picture
Joined: 03/29/2012
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!

psunthud's picture
Joined: 02/14/2012
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]).

mspiegel's picture
Joined: 07/31/2009

Hmm, try:

mxAlgebra(expression = groupdiff / sqrt(S[7,7]), name = "alpha"),

psunthud's picture
Joined: 02/14/2012
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!