Login

News

Univariate ACE models for males and females separately (Heterogenous model)

Posted on
No user picture. lior abramson Joined: 07/21/2017
Forums
Behavioral Genetics Models
Dear openMX users and forum administrators,

I tried to write a script using paths specifications for multigroup ACE univariate analysis with two covariates (age and place of experiment), for boys and girls separately (no opposite sex twins in the analysis). I tried to adapt the script in the user guide for heterogeneous models, which is a general model that does not refer specifically to twins analyses.
In the heterogenous model, I specified different labels for boys and girls for the means of the observed variables as well as the effects of A C and E on the observed variables. I managed to run the script, but the coefficient estimates that I got were so small… The only explanation for such results, I think, is that something is wrong with the syntax, especially when the homogeneity model results turned out fine. The problem is, I cannot figure out what went wrong…
I attach the script and the output. Any help would be very much appreciated.
Thank you!

SCRIPT:

#Heterogenity model
#telling R that place of experiment is a nominal variable
D.impWide $ LabHome6 <- mxFactor(D.impWide $LabHome6, levels=c(0,1), label="LabHome6")
threshold <-mxThreshold(vars="LabHome6", nThresh=1, free=T, values=c(0))

#Finding reasonable start values
meanFAC1 <-mean(na.omit(D.impWide$ OBSERVED_1))
meanFAC2 <-mean(na.omit(D.impWide$ OBSERVED_2))
StartMean <-mean(c(meanFAC1,meanFAC2))
FAC <-c(D.impWide$OBSERVED_1,D.impWide$ OBSERVED_2)
StartVar.ACE <-sqrt(var(na.omit(FAC))/3)

require("OpenMx")

# Select Variables for Analysis
selVars <- c('OBSERVED_1','OBSERVED_2', 'Age6','LabHome6')
aceVars <- c("A1","C1","E1","A2","C2","E2")

# Select Data for Analysis
mzM <- subset(D.impWide, zygosity==1 & sex_a=='sex1',selVars)
mzF <- subset(D.impWide, zygosity==1 & sex_a=='sex2',selVars)
dzM <- subset(D.impWide, zygosity==2 & sex_a=='sex1',selVars)
dzF <- subset(D.impWide, zygosity==2 & sex_a=='sex2',selVars)

# Path objects for Multiple Groups
manifestVars=selVars
latentVars=aceVars

# specify paths for the variances and means of the latent variables
latVariances <- mxPath( from=aceVars, arrows=2, free=FALSE, values=1 )

varCovar <- mxPath( from= c('Age6','LabHome6'), arrows=2, free=(c(TRUE,FALSE)), values=c(1,1), label=c('agevar' ,'placeVar') )

# means of latent variables
latMeans <- mxPath( from="one", to=aceVars, arrows=1, free=FALSE, values=0 )

#specify paths for the means of the observed variables. Observed variables are separate for boys and girls, while boys and girls' covariates (age and place of experiment) receive the same label
obsMeansM <- mxPath( from="one", to=selVars, arrows=1, free=c(T,T,T,F), values=c(StartMean, StartMean, StartMean.Age, 0), labels=c("meanFACm","meanFACm","meanAge" ,"meanPlace") )
obsMeansF <- mxPath( from="one", to=selVars, arrows=1, free=c(T,T,T,F), values=c(StartMean, StartMean, StartMean.Age,0), labels=c("meanFACf","meanFACf","meanAge" ,"meanPlace") )

# path coefficients for twin 1. Separate for boys and girls
pathAceTM1 <- mxPath( from=c("A1","C1","E1", "Age6","LabHome6"), to="OBSERVED_1", arrows=1, free=TRUE, values=StartVar.ACE, label=c("aM","cM","eM", "age", "place") )
pathAceTF1 <- mxPath( from=c("A1","C1","E1", "Age6","LabHome6"), to="OBSERVED_1", arrows=1, free=TRUE, values=StartVar.ACE, label=c("aF","cF","eF", "age","place") )

# path coefficients for twin 2. Separate for boys and girls
pathAceTM2 <- mxPath( from=c("A2","C2","E2", "Age6","LabHome6"), to="OBSERVED_2", arrows=1,free=TRUE, values=StartVar.ACE, label=c("aM","cM","eM", "age","place") )
pathAceTF2 <- mxPath( from=c("A2","C2","E2", "Age6","LabHome6"), to="OBSERVED_2", arrows=1,free=TRUE, values=StartVar.ACE, label=c("aF","cF","eF", "age","place") )

# covariance between C1 & C2
covC1C2 <- mxPath( from="C1", to="C2", arrows=2,free=FALSE, values=1 )
# covariance between A1 & A2 in MZ's
covA1A2_MZ <- mxPath( from="A1", to="A2", arrows=2, free=FALSE, values=1 )
# covariance between A1 & A2 in DZ's
covA1A2_DZ <- mxPath( from="A1", to="A2", arrows=2, free=FALSE, values=.5 )

#call the data.frame with the MZ raw data, mzData, and the DZ raw data
dataMZM <- mxData( observed=mzM, type="raw" )
dataMZF <- mxData( observed=mzF, type="raw" )
dataDZM <- mxData( observed=dzM, type="raw" )
dataDZF <- mxData( observed=dzF, type="raw" )

# Combine Groups
#arrange together all the paths that are the same for MZ and DZ
pathsM <- list(latVariances, latMeans, obsMeansM, pathAceTM1, pathAceTM2, covC1C2, varCovar )
pathsF <- list( latVariances, latMeans, obsMeansF, pathAceTF1, pathAceTF2, covC1C2, varCovar )

#Build a model seperately for MZ and DZ and for each sex
modelMZM <- mxModel(model="MZM", type="RAM", manifestVars=selVars,
latentVars=aceVars, pathsM, threshold, covA1A2_MZ, dataMZM )
modelDZM <- mxModel(model="DZM", type="RAM", manifestVars=selVars,
latentVars=aceVars, pathsM, threshold, covA1A2_DZ, dataDZM )
modelMZF <- mxModel(model="MZF", type="RAM", manifestVars=selVars,
latentVars=aceVars, pathsF, threshold, covA1A2_MZ, dataMZF )
modelDZF <- mxModel(model="DZF", type="RAM", manifestVars=selVars,
latentVars=aceVars, pathsF, threshold, covA1A2_DZ, dataDZF )

minus2ll <- mxAlgebra(expression=MZM.fitfunction + DZM.fitfunction +MZF.fitfunction + DZF.fitfunction, name="minus2loglikelihood" )
obj <- mxFitFunctionAlgebra( "minus2loglikelihood" )
modelACE.Hetero <- mxModel(model="ACE_Hetero", modelMZM, modelDZM,modelMZF, modelDZF, minus2ll, obj )
fitACE.Hetero <- mxTryHardOrdinal(modelACE.Hetero, intervals=TRUE)
sumACE.Hetero <- summary(fitACE.Hetero)

OUTPUT:

Summary of ACE_Hetero

free parameters:
name matrix row col Estimate Std.Er
1 age MZM.A OBSERVED_1 Age6 6.143232e-03 0.012
2 place MZM.A OBSERVED_1 LabHome6 2.172180e-01 0.051
3 aM MZM.A OBSERVED_1 A1 4.955210e-01 0.136
4 cM MZM.A OBSERVED_1 C1 -4.193252e-01 0.139
5 eM MZM.A OBSERVED_1 E1 -5.520680e-01 0.041
6 agevar MZM.S Age6 Age6 7.773803e+00 0.49818
7 meanFACm MZM.M 1 OBSERVED_1 2.937491e+00 1.008
8 meanAge MZM.M 1 Age6 7.931822e+01 0.12634
9 MZM.Thresholds[1,1] MZM.Thresholds 1 LabHome6 -8.392083e-01 0.15320
10 DZM.Thresholds[1,1] DZM.Thresholds 1 LabHome6 -1.068696e+00 0.12395
11 aF MZF.A OBSERVED_1 A1 -6.033133e-07 0.555
12 cF MZF.A OBSERVED_1 C1 5.712114e-01 0.054
13 eF MZF.A OBSERVED_1 E1 7.030873e-01 0.033
14 meanFACf MZF.M 1 OBSERVED_1 2.903178e+00 1.012
15 MZF.Thresholds[1,1] MZF.Thresholds 1 LabHome6 -8.781690e-01 0.16578
16 DZF.Thresholds[1,1] DZF.Thresholds 1 LabHome6 -1.038198e+00 0.11548

Model Statistics:
| Parameters | Degrees of Freedom | Fit (-2lnL units)
Model: 16 1877 5108.699
Saturated: NA NA NA
Independence: NA NA NA
Number of observations/statistics: 487/1893

Information Criteria:
| df Penalty | Parameters Penalty | Sample-Size Adjusted
AIC: 1354.699 5140.699 NA
BIC: -6506.673 5207.711 5156.928
CFI: NA
TLI: 1 (also known as NNFI)
RMSEA: 0 [95% CI (NA, NA)]
Prob(RMSEA <= 0.05): NA
To get additional fit indices, see help(mxRefModels)
timestamp: 2017-08-12 22:26:49
Wall clock time (HH:MM:SS.hh): 00:00:00.33
optimizer: CSOLNP
OpenMx version number: 2.7.15
Need help? See help(mxSummary)

Replied on Mon, 08/21/2017 - 11:50
Picture of user. neale Joined: 07/31/2009

Hi - is it possible that your phenotypes have very small variances? This would give small variance component estimates.

PS Thank you for your patience - summer can slow things down a bit.

Log in or register to post comments
Replied on Mon, 08/21/2017 - 13:31
Picture of user. AdminRobK Joined: 01/24/2014

I don't see anything wrong with your script. Approximately what sort of values are you expecting the path coefficients to be? I'm curious to see the output of the homogeneity model (which apparently looked OK to you), and comparisons of (for example)

cov(mzM, use="pairwise")
modelACE.Hetero$MZM$expectation$output
Log in or register to post comments
Replied on Mon, 08/21/2017 - 13:49
Picture of user. AdminRobK Joined: 01/24/2014

Come to think of it, I guess it's possible that the optimizer is getting stuck, probably at a local minimum. If that's the case, there are a few things you could try:
1. Mean age is about 79 years? Try dividing the age variable by 100, i.e. convert its measuring units from years to centuries. This is a trick that can improve the numerical stability of the optimization.
2. Switch to a different default optimizer with mxOption(), e.g. mxOption(NULL,"Default optimizer","SLSQP") will switch to SLSQP.
3. Use mxTryHard() or one of its wrapper functions in place of mxRun().
Log in or register to post comments
  Login

News