Login

News

Moderation Model

Posted on
No user picture. JuanJMV Joined: 07/20/2016
Forums
Behavioral Genetics Models
Hi everyone,

I am trying to fit a moderation model. I have modified one script but in this script the moderation variable is the same for twin1 and twin 2 (age). I want to fit a moderation model with 1 variable with different values for twin1 and twin2 (like IQ). I have done all the changes but I am not sure if it is correct. Specially this line:

meanG <- mxAlgebra( expression= cbind((int + DR1),(int + DR2)), name="expMeanG" )

Also, I don’t know if the plot is correct since I don’t know how to add the two variables in this line (I think the plot doesn’t change if I use SL1 or SL2 but I am not sure if it is the correct way to do the plot):

SL <- sort(unique(c(mzData$SL2)))

Thank you so much in advance I am stuck with that a little bit of time. Any suggestion and guidance are very welcome

Here the full script:

# Matrices to store a, c, and e Path Coefficients
pathA <- mxMatrix( type="Full", nrow=nv, ncol=na, free=TRUE, values=.6,
label="a11", name="a" )
pathC <- mxMatrix( type="Full", nrow=nv, ncol=nc, free=TRUE, values=.6,
label="c11", name="c" )
pathE <- mxMatrix( type="Full", nrow=nv, ncol=ne, free=TRUE, values=.6,
label="e11", name="e" )

# Matrices to store the moderating a, c, and e Path Coefficients
modPathA <- mxMatrix( type='Full', nrow=nv, ncol=na, free=TRUE, values=.1,
label="aM11", name="aM" )
modPathC <- mxMatrix( type='Full', nrow=nv, ncol=nc, free=TRUE, values=.1,
label="cM11", name="cM" )
modPathE <- mxMatrix( type='Full', nrow=nv, ncol=ne, free=TRUE, values=.1,
label="eM11", name="eM" )

# Matrices for the expected means vector
intercept <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1,
label="interc", name="int" )
PathM <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=0,
label=c("m11"), name="m" )

# Matrices to compute the non-moderated A, C, and E Variance Components
covA <- mxAlgebra( expression=a %*% t(a), name="nA" )
covC <- mxAlgebra( expression=c %*% t(c), name="nC" )
covE <- mxAlgebra( expression=e %*% t(e), name="nE" )

# Algebra to compute the total variance per twin
covP <- mxAlgebra( expression=nA+nC+nE, name="nV" )

# Matrix for the moderator variable
mod1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.SL1"), name="Mod1")
mod2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.SL2"), name="Mod2")

# Matrices to compute the moderated A, C, and E variance components

covAmod1 <- mxAlgebra( expression=(a+ Mod1%*%aM) %*% t(a+ Mod1%*%aM), name="A" )
covCmod1 <- mxAlgebra( expression=(c+ Mod1%*%cM) %*% t(c+ Mod1%*%cM), name="C" )
covEmod1 <- mxAlgebra( expression=(e+ Mod1%*%eM) %*% t(e+ Mod1%*%eM), name="E" )
covAmod2 <- mxAlgebra( expression=(a+ Mod2%*%aM) %*% t(a+ Mod2%*%aM), name="AI2" )
covCmod2 <- mxAlgebra( expression=(c+ Mod2%*%cM) %*% t(c+ Mod2%*%cM), name="CI2" )
covEmod2 <- mxAlgebra( expression=(e+ Mod2%*%eM) %*% t(e+ Mod2%*%eM), name="EI2" )
covAmod12 <- mxAlgebra( expression=(a+ Mod1%*%aM) %*% t(a+ Mod2%*%aM), name="AI12" )
covCmod12 <- mxAlgebra( expression=(c+ Mod1%*%cM) %*% t(c+ Mod2%*%cM), name="CI12" )
covAmod21 <- mxAlgebra( expression=(a+ Mod2%*%aM) %*% t(a+ Mod1%*%aM), name="AI21" )
covCmod21 <- mxAlgebra( expression=(c+ Mod2%*%cM) %*% t(c+ Mod1%*%cM), name="CI21" )

# Algebra to compute the total variance per twin
covPmod <- mxAlgebra( expression=A+C+E, name="V" )

# Algebra for the expected mean vector and expected variance/covariance matrices
covMZ <- mxAlgebra( expression= rbind ( cbind(A+C+E , A+C),
cbind(A+C, A+C+E)), name="expCovMZ" )
covDZ <- mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C, A+C+E)), name="expCovDZ" )
wmod1 <- mxAlgebra( expression= m %*% Mod1, name="DR1" )
wmod2 <- mxAlgebra( expression= m %*% Mod2, name="DR2" )

meanG <- mxAlgebra( expression= cbind((int + DR1),(int + DR2)), name="expMeanG" )
# Data objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups
objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMeanG",
dimnames=selVars )
objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanG",
dimnames=selVars )
fitMZ <- mxFitFunctionML()
fitDZ <- mxFitFunctionML()

# Combine Groups
parsZ <- list( pathA, pathC, pathE, modPathA, modPathC, modPathE, PathM,
intercept)
modelMZ <- mxModel( parsZ, meanG, covAmod1, covCmod1, covEmod1,covAmod2,covCmod2,covEmod2,covAmod12,covCmod12,covAmod21,covCmod21, covPmod, mod1, mod2,
wmod1,wmod2, covMZ, dataMZ, objMZ, fitMZ, name="MZ" )
modelDZ <- mxModel( parsZ, meanG, covAmod1, covCmod1, covEmod1,covAmod2,covCmod2,covEmod2,covAmod12,covCmod12,covAmod21,covCmod21, covPmod, mod1,mod2,
wmod1,wmod2, covDZ, dataDZ, objDZ, fitDZ, name="DZ" )
minus2ll <- mxAlgebra( expression=MZ.objective+ DZ.objective, name="m2LL" )
obj <- mxFitFunctionAlgebra( "m2LL" )
ModAceModel <- mxModel( "ModACE", parsZ, modelMZ, modelDZ, minus2ll, obj )

#=======================================================================#
# RUN MODEL #
#=======================================================================#

# Run Moderation ACE model
ModAceFit <- mxRun(ModAceModel)
ModAceSumm <- summary(ModAceFit)
ModAceSumm
round(ModAceFit@output$estimate,4)

# Generate Output using Helper Functions
tableFitStatistics(ModAceFit)

# Generate Table of Parameter Estimates using mxEval
pathEstimatesACE <- round(mxEval(cbind(a,c,e,aM,cM,eM), ModAceFit),4)
rownames(pathEstimatesACE) <- 'pathEstimates'
colnames(pathEstimatesACE) <- c('a','c','e','aM','cM','eM')
pathEstimatesACE
varComponentsACE <- round(mxEval(cbind(A/V,C/V,E/V), ModAceFit@submodels$MZ),4)
#varComponentsACE <- round(mxEval(cbind(A,C,E), ModAceFit@submodels$MZ),4)
rownames(varComponentsACE) <- 'varComponents'
colnames(varComponentsACE) <- c('%a^2','%c^2','%e^2')
varComponentsACE
fitGofs(ModAceFit)
fitEstS(ModAceFit)
#=======================================================================#
# PLOT MODERATION #
#=======================================================================#

a <- pathEstimatesACE[,'a']
c <- pathEstimatesACE[,'c']
e <- pathEstimatesACE[,'e']
aM <- pathEstimatesACE[,'aM']
cM <- pathEstimatesACE[,'cM']
eM <- pathEstimatesACE[,'eM']

SL <- sort(unique(c(mzData$SL2)))

Va <- (a + aM*SL)^2
Vc <- (c + cM*SL)^2
Ve <- (e + eM*SL)^2
Vt <- Va + Vc + Ve
SVa <- Va/Vt
SVc <- Vc/Vt
SVe <- Ve/Vt
out <- as.matrix(cbind(Va,Vc,Ve,Vt,SVa,SVc,SVe))
rownames(out)=paste('SL=',1:2,sep='')
out

# par(mfrow=c(2,1))
matplot(SL, out[,1:4], type="l", lty=4:1, col=4:1, xlab="SL",
ylab="Variance Component", main="SL Moderation - unstandardized
variance components")
legend(2.7, 1.2, c("Va","Vc","Ve","Vt"), lty=4:1, col=4:1)

Replied on Sat, 12/23/2017 - 05:15
No user picture. JuanJMV Joined: 07/20/2016

Thank you so much Rob! The document and the script are extremely helpful. Now I am a little bit less confused.
I have been working on it and I have fitted the full bivariate model with 17parameters. But I still have some doubts.

1)I don’t know if I have to tried to fit the extended univariate moderation model. I think in my data set M and T are correlated but I don’t know how I can check it.
In your script you checked it with the true values for "a0m","c0m","e0m…." but how can I get the true values with my data?

2) Also, I don’t know how can I get the results from my model, I mean I want to know a c e for SL=1 and for SL=2. Have I to create something like in that? But what values should I use?:

a <- pathEstimatesACE[,'a']
c <- pathEstimatesACE[,'c']
e <- pathEstimatesACE[,'e']
aM <- pathEstimatesACE[,'aM']
cM <- pathEstimatesACE[,'cM']
eM <- pathEstimatesACE[,'eM']

SL <- sort(unique(c(mzData$SL)))

Va <- (a + aM*SL)^2
Vc <- (c + cM*SL)^2
Ve <- (e + eM*SL)^2
Vt <- Va + Vc + Ve
SVa <- Va/Vt
SVc <- Vc/Vt
SVe <- Ve/Vt
out <- as.matrix(cbind(Va,Vc,Ve,Vt,SVa,SVc,SVe))
rownames(out)=paste('SL=',1:2,sep='')

3)I am trying to learn it fitting several models but now I am focused in one of them where the moderation variable is dichotomous. This model (the script below) is ok for that? If the trait is dichotomous or ordinal I have to use a threshold model?

4) the order of the variables in the data set should be moderation variable 1, trait 1, moderation variable 2, trait 2, right?

Thank you so much Rob and sorry for the inconveniences I am not very good with that but I need it and I like it!
Here the full script adapted for my data:

# Select trait Variables
varsc <- c('RU') # list of continuous variables names
nvc <- 1 # number of continuous variables
ntvc <- nvc*2 # number of total continuous variables
conVars <- paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="")

# Select mod Variables
nth <- 1 # number of thresholds
varso <- c('SL') # list of ordinal variables names
nvo <- 1 # number of ordinal variables
ntvo <- nvo*2 # number of total ordinal variables
ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="")

# Select Variables for Analysis
vars <- c("SL","RU") # list of variables names
nv <- nvc+nvo # number of variables
ntv <- nv*2 # number of total variables
oc <- c(0,1) # 0 for continuous, 1 for ordinal variables
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")

datmz1 <- data[data$Zyg==1, c(selVars)]
datdz1 <- data[data$Zyg==2, c(selVars)]

datmz <- datmz1[!is.na(datmz1$SL1),]
datdz <- datdz1[!is.na(datdz1$SL2),]
datmz <- datmz1[!is.na(datmz1$SL2),]
datdz <- datdz1[!is.na(datdz1$SL1),]

describe(data)
dim(data)
head(data)#

# big model MZ and DZ
#
modmodelMZ=mxModel("modmodMZ",
# def .. Definition variables which we need to create the moderated paths (from M->T)
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE,
labels=c("data.SL1"), name="m1" ) ,
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE,
labels=c("data.SL2"), name="m2" ) ,
# ------------------------------------------------------------ the means .. Not starting values are good.
mxMatrix( type="Full", nrow=1, ncol=4, free=TRUE, values=0,
labels=c("meanm","meant","meanm","meant"), name="expMean" ),
# ------------------------------------------------------------------------
# ACE model moderator - cholesky elements. now defined as scalars (1x1), called lower for no good reason.
# moderator parameters
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="a0m", name="A0M" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="c0m", name="C0M" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="e0m", name="E0M" ),
# ACE model trait 0 trait parameters: main effects
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="a0t", name="A0T" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="c0t", name="C0T" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="e0t", name="E0T" ),
# ACE model trait 1 moderator effects
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="a1t", name="A1T" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="c1t", name="C1T" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="e1t", name="E1T" ),
#
# ACE model moderator-trait 0 moderated corvariances. First main effects
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="a0mt", name="A0MT" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="c0mt", name="C0MT" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="e0mt", name="E0MT" ),
# ACE model moderator-trait 1 moderated covariances.. moderation
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="a1mt", name="A1MT" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="c1mt", name="C1MT" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="e1mt", name="E1MT" ),

# assumble parameters in to a cholesky decomp for twin 1 and twin 2
# do this in three parts : A, C, and E.
#
mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values=
c(1,0,1,0,1,0,1,1,0,1),name="PsAmz"),
mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values=
c(1,0,1,0,1,0,1,1,0,1),name="PsC"),
#
# Matrices generated to hold A and E computed Variance Components
# A This is a Cholesky decomposition of A for twin 1 and twin 2 (mz and dz)
# not that m1 appear in the top left part (tinw 1) and m2 in the bottom right part (twin 2)
# A..
mxAlgebra( rbind( cbind(A0M,0,0,0),
cbind(A0MT+m1%x%A1MT,A0T+m1%x%A1T,0,0),
cbind(0,0,A0M,0),
cbind(0,0,A0MT+m2%x%A1MT,A0T+m2%x%A1T)
), name="chA"),
# C chol
mxAlgebra( rbind(
cbind(C0M,0,0,0),
cbind(C0MT+m1%x%C1MT,C0T+m1%x%C1T,0,0),
cbind(0,0,C0M,0),
cbind(0,0,C0MT+m2%x%C1MT,C0T+m2%x%C1T)
), name="chC"),
# E chol
mxAlgebra( rbind(
cbind(E0M,0,0,0),
cbind(E0MT+m1%x%E1MT,E0T+m1%x%E1T,0,0),
cbind(0,0,E0M,0),
cbind(0,0,E0MT+m2%x%E1MT,E0T+m2%x%E1T)
), name="chE"),
#
# get the expected covariance matrices
mxAlgebra(chA%*%PsAmz%*%t(chA),name="Amz" ), # variance component A
mxAlgebra(chC%*%PsC%*%t(chC),name="C" ), # variance component A
mxAlgebra(chE%*%t(chE),name="E" ), # variance component A
#
# and assemble in to phenotypic 4x4 covariance matrix of MZ and DZ
mxAlgebra(Amz+C+E, name="expCovMZ" ),
mxData( observed=datmz, type="raw" ),
mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars)
)
# Same as above
modmodelDZ=mxModel("modmodDZ",
# def
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE,
labels=c("data.SL1"), name="m1" ) ,
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE,
labels=c("data.SL2"), name="m2" ) ,
# ------------------------------------------------------------ means
mxMatrix( type="Full", nrow=1, ncol=4, free=TRUE, values=0,
labels=c("meanm","meant","meanm","meant"), name="expMean" ),
# ------------------------------------------------------------------------
# ACE model moderator
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="a0m", name="A0M" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="c0m", name="C0M" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="e0m", name="E0M" ),
# ACE model trait 0
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="a0t", name="A0T" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="c0t", name="C0T" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="e0t", name="E0T" ),
# ACE model trait 1
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="a1t", name="A1T" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="c1t", name="C1T" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="e1t", name="E1T" ),

# ACE model moderator-trait 0
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="a0mt", name="A0MT" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="c0mt", name="C0MT" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.6,
label="e0mt", name="E0MT" ),
# ACE model moderator-trait 1
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="a1mt", name="A1MT" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="c1mt", name="C1MT" ) , # path coefficient A
mxMatrix( type="Lower", nrow=1, ncol=1, free=TRUE, values=.06,
label="e1mt", name="E1MT" ),

# ,
mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values=
c(1,0,.5,0,1,0,.5,1,0,1),name="PsAdz"),
mxMatrix(type="Symm", nrow=4, ncol=4, free=F, values=
c(1,0,1,0,1,0,1,1,0,1),name="PsC"),

# Matrices generated to hold A and E computed Variance Components
# A
mxAlgebra( rbind( cbind(A0M,0,0,0),
cbind(A0MT+m1%x%A1MT,A0T+m1%x%A1T,0,0),
cbind(0,0,A0M,0),
cbind(0,0,A0MT+m2%x%A1MT,A0T+m2%x%A1T)
), name ="chA"),
# C chol
mxAlgebra( rbind(
cbind(C0M,0,0,0),
cbind(C0MT+m1%x%C1MT,C0T+m1%x%C1T,0,0),
cbind(0,0,C0M,0),
cbind(0,0,C0MT+m2%x%C1MT,C0T+m2%x%C1T)
), name ="chC"),
# E chol
mxAlgebra( rbind(
cbind(E0M,0,0,0),
cbind(E0MT+m1%x%E1MT,E0T+m1%x%E1T,0,0),
cbind(0,0,E0M,0),
cbind(0,0,E0MT+m2%x%E1MT,E0T+m2%x%E1T)
), name ="chE"),
#
#
mxAlgebra(chA%*%PsAdz%*%t(chA),name="Adz" ), # variance component A
mxAlgebra(chC%*%PsC%*%t(chC),name="C" ), # variance component A
mxAlgebra(chE%*%t(chE),name="E" ), # variance component A
#
# Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
# this should look familiar! expected MZ covariance matrix as function of A,E
mxAlgebra(Adz+C+E, name="expCovDZ" ),
mxData( observed=datdz, type="raw" ),
mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars)
)

modmodelA=mxModel("modmodel",modmodelMZ, modmodelDZ,
mxAlgebra( modmodMZ.objective + modmodDZ.objective, name="m2LL" ),
mxAlgebraObjective( "m2LL" ))

modmodelAfit=mxTryHardOrdinal(modmodelA)
tableFitStatistics(modmodelAfit)

Log in or register to post comments
Replied on Fri, 01/05/2018 - 13:38
Picture of user. AdminRobK Joined: 01/24/2014

1)I don’t know if I have to tried to fit the extended univariate moderation model. I think in my data set M and T are correlated but I don’t know how I can check it.
In your script you checked it with the true values for "a0m","c0m","e0m…." but how can I get the true values with my data?

First off, I can't take credit (or blame) for the example script, since I didn't write it--Conor V. Dolan did. Secondly, you may find Conor's slides from Boulder 2016 useful. An important take-home message is that you should always use the "full bivariate" model (in which the moderator is an endogenous random regressor), unless (1) the moderator correlates 1.0 between twins, (2) the moderator correlates 0 between twins, or (3) the covariance between the moderator and trait is equally due to A, C, and E. In practice, cases (2) and (3) are exceedingly rare, and case (3) can really only be checked by model-fitting anyhow.

2) Also, I don’t know how can I get the results from my model, I mean I want to know a c e for SL=1 and for SL=2. Have I to create something like in that? But what values should I use?:

Once you have parameter estimates from your selected moderation model, you can plot the values of standardized and raw A, C, and E variance components at different values of the moderator. Take a look at Figure 2 in my relevant publication. The exact code used to create those plots is:

tiff(filename="C:/Grad school/Biometric moderation project/manuscript/second revision/Figures/Fig2a.tif",
width=6, height=6, units="in",
type="cairo",compression="lzw",
pointsize=10,
res=300
)
x <- (0:100)/100
ya <- (cbind(1,x) %*% run_15@output$estimate[c("gammaA0","gammaA2")])^2
yc <- rep(run_15@output$estimate[c("gammaC0")]^2,101)
ye <- (cbind(1,x) %*% run_15@output$estimate[c("gammaE0","gammaE2")])^2
plot(ya ~ x, type="l", lwd=2, ylab="Variance", xlab="SES", ylim=c(0,135), col="red", main="A")
points(yc ~ x, type="l", col="forestgreen", lwd=2, lty="dashed")
points(ye ~ x, type="l", col="blue", lwd=2, lty="dotted")
legend(x=0, y=135, legend=c("A","C","E"), lty=c("solid","dashed","dotted"),
col=c("red", "forestgreen", "blue"))
dev.off()
tiff(filename="C:/Grad school/Biometric moderation project/manuscript/second revision/Figures/Fig2b.tif",
width=6, height=6, units="in",
type="cairo",compression="lzw",
pointsize=10,
res=300
)
yt <- ya+yc+ye
yap <- ya/yt
ycp <- yc/yt
yep <- ye/yt
plot(yap ~ x, type="l", lwd=2, ylab="Variance", xlab="SES", ylim=c(0,1), col="red", main="B")
points(ycp ~ x, type="l", col="forestgreen", lwd=2, lty="dashed")
points(yep ~ x, type="l", col="blue", lwd=2, lty="dotted")
legend(x=0, y=1, legend=c("A","C","E"), lty=c("solid","dashed","dotted"),
col=c("red", "forestgreen", "blue"))
dev.off()

The scale of your moderator and trait will be different, the symbol for your MxModel will be different, the labels for your free parameters will be different, and the function of the free parameters that maps x to y will be different, but the essential idea is there.

3)I am trying to learn it fitting several models but now I am focused in one of them where the moderation variable is dichotomous. This model (the script below) is ok for that? If the trait is dichotomous or ordinal I have to use a threshold model?

I haven't looked over your script yet., but if you're fitting a moderation model in which the moderator itself is biometrically decomposed, then treating an ordinal moderator as such (i.e., using thresholds) is a good idea. However, you might also be able to side-step all this Purcell-style moderation modeling with savvy use of multiple groups. You'd have the following groups: MZ twins, concordant for score of 1 on the moderator; MZ twins, concordant 0; MZ twins, discordant; DZ twins, concordant 1; DZ twins, concordant 0; DZ twins, discordant. I'm not entirely sure of the details offhand, but in principle this should work.

4) the order of the variables in the data set should be moderation variable 1, trait 1, moderation variable 2, trait 2, right?

Yes, that works.

Log in or register to post comments
Replied on Wed, 01/24/2018 - 10:33
Picture of user. AdminRobK Joined: 01/24/2014

This reminds me...two miscellaneous points.

First, replace

mxAlgebra( modmodMZ.objective + modmodDZ.objective, name="m2LL" ),
mxAlgebraObjective( "m2LL" )

with

mxFitFunctionMultigroup(c("modmodMZ","modmodDZ"))

. They are nearly equivalent, but the important thing is that an MxFitFunctionMultigroup is compatible with certain utility functions in OpenMx, whereas an algebra fitfunction isn't.

Second, does your script actually run? I had thought that it would be necessary to enter the moderator, 'SL', into the dataset twice, with one copy used as an endogenous variable, and the other as a "definition" variable. But I could be wrong about that. It looks like you're using the same column as an endogenous variable AND as a definition variable, right?

Log in or register to post comments
  Login

News