# GxE with ordinal manifest in RAM

I'm trying to build a Purcell style GxE model in RAM. I have not found any script for this so I created it from scratch. Right now I have a two ordinal manifests (twin 1, twin 2, same phenotype, 4 cats, 3threshs). The moderator will be continuous in the end. However, right now for experimental reasons I have created a binary moderator (sex) for me to easily check if my model does produce accurate results.

I have build two models: 1) All thresholds free, two mxconstraints, means=0, variances=1 and 2) 2 thresholds fixed (0,1), third one free, mean free, variance free, no constraint.

1) is working almost perfectly, at least regarding the results. However, I don't really want to rely on this method as it is kind of messy. And I'm really unsure if I set all the ubounds, lbounds and constraints right. All I can really tell is the results are right.

2) produces sometimes similar results as 1) but does differ wildly often. Particularly the standard errors go crazy, sometimes as much as 500x estimate.

I was wondering if I could get some feedback what would have to be changed.

Here are both scripts:

1):

ordDVs = c("pop0100_w1_T1", "pop0100_w1_T2")

```
```cgr2[, ordDVs] = mxFactor(cgr2[, ordDVs], 1:4)

# subset MZ DZ data

mzData = subset(cgr2, zyg0102w1 == 1)

dzData = subset(cgr2, zyg0102w1 == 2)

latents = c("A1", "C1", "E1", "A2", "C2", "E2")

manifests = ordDVs

latVariances <- mxPath( from=latents, arrows=2,

free=F, values=1)

# means of latent variables

latMeans <- mxPath( from="one", to=latents, arrows=1,

free=F, values=0 )

# means of observed variables

obsMeans <- mxPath( from="one", to=manifests, arrows=1,

free=F, values=0, labels="mean[1,1]")

# path coefficients for twin 1

pathAceT1 <- mxPath( from=c("A1","C1","E1"), to="pop0100_w1_T1", arrows=1,

free=F, values = 1, label=c("patha[1,1]","pathc[1,1]","pathe[1,1]"))

# path coefficients for twin 2

pathAceT2 <- mxPath( from=c("A2","C2","E2"), to="pop0100_w1_T2", arrows=1,

free=F, values = 1, label=c("patha[1,1]","pathc[1,1]","pathe[1,1]"))

# covariance between C1 & C2

covC1C2 <- mxPath( from="C1", to="C2", arrows=2,

free=FALSE, values=1 )

# covariance between A1 & A2 in MZ twins

covA1A2_MZ <- mxPath( from="A1", to="A2", arrows=2,

free=FALSE, values=1 )

# covariance between A1 & A2 in DZ twins

covA1A2_DZ <- mxPath( from="A1", to="A2", arrows=2,

free=FALSE, values=.5 )

beta1 <- mxMatrix(type = "Full", 1,1, free = F, labels = "beta1", values = 0)

mean1 <- mxMatrix(type = "Full", 1,1, free = F, labels = "mean1", values = 0)

blub <- mxAlgebra(name="mean", expression = mean1+beta1*data.sex_w1_T1)

a <- mxMatrix(type = "Full", 1,1, free = T, labels = "a", lbound = 0)

beta2 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta2", ubound = 1, lbound = -1)

patha <- mxAlgebra(name="patha", expression = a+beta2*data.sex_w1_T1)

c <- mxMatrix(type = "Full", 1,1, free = T, labels = "c", lbound = 0)

beta3 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta3", ubound = 1, lbound = -1)

pathc <- mxAlgebra(name="pathc", expression = c+beta3*data.sex_w1_T1)

e <- mxMatrix(type = "Full", 1,1, free = T, labels = "e", lbound = 0)

beta4 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta4", ubound = 1, lbound = -1)

pathe <- mxAlgebra(name="pathe", expression = e+beta4*data.sex_w1_T1)

paths <- list( latVariances, latMeans, obsMeans,

pathAceT1, pathAceT2, covC1C2, blub, beta1, mean1,

a, beta2, patha, c, beta3, pathc, e, beta4, pathe)

dataMZ <- mxData(observed=mzData, type = "raw")

dataDZ <- mxData(observed=dzData, type = "raw")

threshold <- mxThreshold(vars=c("pop0100_w1_T1", "pop0100_w1_T2"), nThresh=c(3, 3), free = c(T, T, T),

labels = c("true1", "true2", "true3"), values = c(0,1,2))

A <- mxAlgebra( expression=patha*patha, name="A2" )

C <- mxAlgebra( expression=pathc*pathc, name="C2" )

E <- mxAlgebra( expression=pathe*pathe, name="E2" )

V <- mxAlgebra( expression=A2+C2+E2, name="V" )

a2 <- mxAlgebra( expression=A2/V, name="a2")

c2 <- mxAlgebra( expression=C2/V, name="c2")

e2 <- mxAlgebra( expression=E2/V, name="e2")

constraint <- mxConstraint(a*a+c*c+e*e == 1, name = "constraintordinal")

constraint2 <- mxConstraint((beta2+a)*(beta2+a)+(beta3+c)*(beta3+c)+(beta4+e)*(beta4+e) == 1, name = "constraintordinal2")

modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=manifests,

latentVars=latents, paths, covA1A2_MZ, dataMZ, threshold,

A, C, E, V, a2, c2, e2, constraint, constraint2)

modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=manifests,

latentVars=latents, paths, covA1A2_DZ, dataDZ, threshold,

A, C, E, V, a2, c2, e2)

obj <- mxFitFunctionMultigroup(c("MZ", "DZ"))

modelACE <- mxModel("ACE", modelMZ, modelDZ, obj)

fitACE <- mxRun(modelACE)

`summary(fitACE)`

script 2:

ordDVs = c("pop0100_w1_T1", "pop0100_w1_T2")

```
```cgr2[, ordDVs] = mxFactor(cgr2[, ordDVs], 1:4)

# subset MZ DZ data

mzData = subset(cgr2, zyg0102w1 == 1)

dzData = subset(cgr2, zyg0102w1 == 2)

latents = c("A1", "C1", "E1", "A2", "C2", "E2")

manifests = ordDVs

latVariances <- mxPath( from=latents, arrows=2,

free=T, values=1, labels = c("AV", "CV", "EV", "AV", "CV", "EV"))

# means of latent variables

latMeans <- mxPath( from="one", to=latents, arrows=1,

free=F, values=0 )

# means of observed variables

obsMeans <- mxPath( from="one", to=manifests, arrows=1,

free=F, values=0, labels="mean[1,1]")

# path coefficients for twin 1

pathAceT1 <- mxPath( from=c("A1","C1","E1"), to="pop0100_w1_T1", arrows=1,

free=F, values = 1, label=c("patha[1,1]","pathc[1,1]","pathe[1,1]"))

# path coefficients for twin 2

pathAceT2 <- mxPath( from=c("A2","C2","E2"), to="pop0100_w1_T2", arrows=1,

free=F, values = 1, label=c("patha[1,1]","pathc[1,1]","pathe[1,1]"))

# covariance between C1 & C2

covC1C2 <- mxPath( from="C1", to="C2", arrows=2,

free=FALSE, values=1 )

# covariance between A1 & A2 in MZ twins

covA1A2_MZ <- mxPath( from="A1", to="A2", arrows=2,

free=FALSE, values=1 )

# covariance between A1 & A2 in DZ twins

covA1A2_DZ <- mxPath( from="A1", to="A2", arrows=2,

free=FALSE, values=.5 )

beta1 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta1", values = 0)

mean1 <- mxMatrix(type = "Full", 1,1, free = T, labels = "mean1", values = 0)

blub <- mxAlgebra(name="mean", expression = mean1+beta1*data.sex_w1_T1)

a <- mxMatrix(type = "Full", 1,1, free = T, labels = "a", lbound = 0)

beta2 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta2", ubound = 1, lbound = -1)

patha <- mxAlgebra(name="patha", expression = a+beta2*data.sex_w1_T1)

c <- mxMatrix(type = "Full", 1,1, free = T, labels = "c", lbound = 0)

beta3 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta3", ubound = 1, lbound = -1)

pathc <- mxAlgebra(name="pathc", expression = c+beta3*data.sex_w1_T1)

e <- mxMatrix(type = "Full", 1,1, free = T, labels = "e", lbound = 0)

beta4 <- mxMatrix(type = "Full", 1,1, free = T, labels = "beta4", ubound = 1, lbound = -1)

pathe <- mxAlgebra(name="pathe", expression = e+beta4*data.sex_w1_T1)

paths <- list( latVariances, latMeans, obsMeans,

pathAceT1, pathAceT2, covC1C2, blub, beta1, mean1,

a, beta2, patha, c, beta3, pathc, e, beta4, pathe)

dataMZ <- mxData(observed=mzData, type = "raw")

dataDZ <- mxData(observed=dzData, type = "raw")

threshold <- mxThreshold(vars=c("pop0100_w1_T1", "pop0100_w1_T2"), nThresh=c(3, 3), free = c(F, F, T),

labels = c("true1", "true2", "true3"), values = c(0,1,2))

A <- mxAlgebra( expression=patha*patha, name="A2" )

C <- mxAlgebra( expression=pathc*pathc, name="C2" )

E <- mxAlgebra( expression=pathe*pathe, name="E2" )

V <- mxAlgebra( expression=A2+C2+E2, name="V" )

a2 <- mxAlgebra( expression=A2/V, name="a2")

c2 <- mxAlgebra( expression=C2/V, name="c2")

e2 <- mxAlgebra( expression=E2/V, name="e2")

modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=manifests,

latentVars=latents, paths, covA1A2_MZ, dataMZ, threshold,

A, C, E, V, a2, c2, e2)

modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=manifests,

latentVars=latents, paths, covA1A2_DZ, dataDZ, threshold,

A, C, E, V, a2, c2, e2)

obj <- mxFitFunctionMultigroup(c("MZ", "DZ"))

modelACE <- mxModel("ACE", modelMZ, modelDZ, obj)

fitACE <- mxRun(modelACE)

summary(fitACE)

tryhard = mxTryHardOrdinal(modelACE, exhaustive = T, extraTries = 50)

`summary(tryhard)`

any help would be much appreciated.

## model identification

latVariances <- mxPath( from=latents, arrows=2,

free=T, values=1, labels = c("AV", "CV", "EV", "AV", "CV", "EV"))

You want to fix those variance to 1.0, right? Change `free=T` to `free=F`.

Log in or register to post comments

In reply to model identification by AdminRobK

## I thought I could estimate

Log in or register to post comments

In reply to I thought I could estimate by Leo

## manifest variables

Log in or register to post comments

In reply to manifest variables by AdminRobK

## ah okay, thank you, now I

Log in or register to post comments

In reply to ah okay, thank you, now I by Leo

## See my comment here.

Log in or register to post comments

## moderator?

Log in or register to post comments

In reply to moderator? by AdminRobK

## it does differ between twins,

Log in or register to post comments