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)
Printer-friendly version
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.
Thanks! I think that it works now...
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)
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 ofmxRun().Thank you!
I followed you advices and I think that it works now.