You are here

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

6 posts / 0 new
Last post
lior abramson's picture
Offline
Joined: 07/21/2017 - 13:13
Univariate ACE models for males and females separately (Heterogenous model)

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) 
neale's picture
Offline
Joined: 07/31/2009 - 15:14
Scale of input data?

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.

lior abramson's picture
Offline
Joined: 07/21/2017 - 13:13
thank you

Thanks! I think that it works now...

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
What should the path coefficients look like?

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
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
optimization issue?

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

lior abramson's picture
Offline
Joined: 07/21/2017 - 13:13
thank you

Thank you!
I followed you advices and I think that it works now.

Log in or register to post comments