You are here

ACE decomposition of growth factors code

15 posts / 0 new
Last post
ethibodeau's picture
Offline
Joined: 05/16/2018 - 17:05
ACE decomposition of growth factors code

Does anyone have example path specification code to run a basic linear growth model where the slope and intercept factors are decomposed into ACE based on twin data? Thanks!

ethibodeau's picture
Offline
Joined: 05/16/2018 - 17:05
cholesky approach

I should specify, using the cholesky bivariate approach for the slope and intercept association modeling

AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
No, but...

I have emailed people who work in this area to find out. It seems to me that what is needed would be a combination of the umxRAM (which is used by umx to specify the non-twin LGC models) and umxACE which is matrix-based. It seems to me that a genetic LGC function ought to allow NGrowthFactors (to set up level, slope and quadratic at least), Nvar variables and Nocc occasions. None of this requires more than a few hours of scripting. However, to build a function for submission to umx (which could prove really useful to the scientific community) would take a bit more effort, but hopefully not too much.

Obviously, allowing raw data input would be much better as longitudinal studies that claim not to have any missing data should be viewed with deep suspicion. Better still is to avoid the biases inherent in declaring people to be of the same exact age at testing when such precise logistics are not really achievable with human research subjects. Doing that takes us to using observed ages to replace the factor loadings for each case individually - easily done with definition variables, see, e.g., Schmitt et al 2018, which features the OP's desired model, but is specified with matrices instead of paths.

HTH
Mike

ethibodeau's picture
Offline
Joined: 05/16/2018 - 17:05
Maybe I should clarify..

Thanks for the response. I have one continuous variable measured across 4 waves. I'm working with the raw data, if I understand you correctly. Wide format data. The variable measured at each wave, age at each wave, zygosity, etc. There is missing data.

I want to run a simple linear growth model, indeed using the age definition variable approach. That is, as you indicate, fixing the slope factor loadings to the age at each wave for each person. So this would be a linear LGM with individually-varying times of observation. Generally I freely estimate the indicator residual variances across time, but constrained to equality between twins. I also constrain the growth factor means to equality between twin pairs.

Then I want to decompose the slope and intercept factors into ACE factors (freely estimating factor loadings, setting factor variances to 1, means to 0, setting A factor covariances between twins to 1 for MZ and .5 for DZ, setting C factor covariances between twins to 1, setting residual variance of the slope and intercept to 0 implying that the ACE factors explain all the variance, etc) and using the Cholesky approach to model the bivariate nature of the slope and intercept. That is, regressing the slope onto the intercept ACE factors. I am not decomposing the indicator residual variances into ACE factors. I have done this in Mplus and I want to replicate in OpenMx for piece of mind and to use some of OpenMx's nice features. :)

I was trying to meld the Openmx example code of univariate ACE decomp with openmx's example code of linear LGM, but I got myself confused.

any help would be appreciated!

ethibodeau's picture
Offline
Joined: 05/16/2018 - 17:05
Also..

Im happy to do this in matrix specification if it's easier. The only reason I mentioned path was to make it easier to compare to Mplus, which is only path specification as you are well aware.

ethibodeau's picture
Offline
Joined: 05/16/2018 - 17:05
What do people think of this code?

The output essentially matches that of MPlus, however, the only difference is the factor loading for the 'E' factor loading onto the intercept is the same number (magnitude), but the sign is negative in OpenMx and positive in Mplus.

### OpenMx script for linear latent growth model with ACE decomposition and Cholesky approach 
 
 
## Load data
 
mnmj_30d <- read.csv(file ="/Desktop/mnmj_30d_twin.csv", header=TRUE, sep=",")
 
 
## Rename data to get rid of dots in variable names
 
names(mnmj_30d)<-c("X","famid","famid_id","famid_school","psuscid",      
                   "region","region_na","gswgt1","gswgt1_na","cohort","cohort_na","mny", 
                   "mno", "mny_na","mno_na","zyg","mj1use","sex","data","mj11", "mj21",    
                   "mj31","mj41","age11","age21","age31","age41","mj12","mj22",     
                   "mj32","mj42","age12","age22","age32","age42" )
 
 
## subset to relevant data
 
mj_twins<-mnmj_30d[c("mj11", "mj21","mj31","mj41",
                     "age11","age21","age31","age41",
                     "mj12","mj22","mj32","mj42",
                     "age12","age22","age32","age42", "zyg")]
 
 
## Load various packages 
 
library(OpenMx)
 
## center age at 16.5yrs (mean age of first use) for twin 1 and twin 2 ages
 
mj_twins$age11<-mj_twins$age11-16.5
mj_twins$age21<-mj_twins$age21-16.5
mj_twins$age31<-mj_twins$age31-16.5
mj_twins$age41<-mj_twins$age41-16.5
 
mj_twins$age12<-mj_twins$age12-16.5
mj_twins$age22<-mj_twins$age22-16.5
mj_twins$age32<-mj_twins$age32-16.5
mj_twins$age42<-mj_twins$age42-16.5
 
# subset data by zygosity
mj_twins_mz<-mj_twins[which(mj_twins$zyg==1),]
mj_twins_dz<-mj_twins[which(mj_twins$zyg==2),]
 
cholesky_ace_bivariate_multigroup_lgm <- mxModel('Multiple group, Path Specification',
 
                                                 ### MZ Growth Model Specification
 
 
                                                 mxModel('mzmodel',type='RAM', mxData(observed=mj_twins_mz, type='raw'),
                                                         manifestVars=c("mj11", "mj21",    
                                                                        "mj31","mj41",
                                                                        "mj12", "mj22",    
                                                                        "mj32","mj42"),
                                                         latentVars=c(
                                                           # growth factors for twin 1 and twin 2
                                                           'intercept1','slope1','intercept2','slope2',
                                                           # Intercept ACE factors for twin 1 and twin 2
                                                           "AI1","CI1","EI1","AI2","CI2","EI2",
                                                           # Slope ACE factors for twin 1 and twin 2
                                                           "AS1","CS1","ES1","AS2","CS2","ES2"),
 
 
                                                         # twin1 residual variance paths, freely estimated across time, constrained to equality between twin 1 and twin 2 models
                                                         mxPath(from=c("mj11", "mj21",    
                                                                       "mj31","mj41"),
                                                                arrows=2, free=TRUE, values=10, labels=c('rv1','rv2','rv3','rv4')),
 
                                                         # twin2 residual variance paths, freely estimated across time,constrained to equality between twin 1 and twin 2 models
                                                         mxPath(from=c("mj12", "mj22",    
                                                                       "mj32", "mj42"),
                                                                arrows=2, free=TRUE, values=10,labels=c('rv1','rv2','rv3','rv4')),
 
                                                        # twin1 intercept residual variance set to 0 (all consumed by ACE factors)
                                                         mxPath(from=c('intercept1'),arrows=2, connect='unique.pairs',free=F, labels=c('intercept_variance'),values=0),
 
                                                         # twin2 intercept residual variance set to 0 (all consumed by ACE factors)
                                                         mxPath(from=c('intercept2'),arrows=2, connect='unique.pairs',free=F, labels=c('intercept_variance'),values=0),
 
                                                         # twin1 slope residual variance set to 0 (all consumed by ACE factors)
                                                         mxPath(from=c('slope1'),arrows=2, connect='unique.pairs',free=F, labels=c('slope_variance'),values=0),
 
                                                         # twin2 slope residual variance set to 0 (all consumed by ACE factors)
                                                         mxPath(from=c('slope2'),arrows=2, connect='unique.pairs',free=F, labels=c('slope_variance'),values=0),
 
 
                                                         # twin1 set intercept factor loadings set to 1
                                                         mxPath(from='intercept1', to=c("mj11", "mj21",    
                                                                                        "mj31","mj41"),
                                                                arrows=1, free=FALSE, values=1),
 
                                                         # twin2 set intercept factor loadings set to 1
                                                         mxPath(from='intercept2', to=c("mj12", "mj22",    
                                                                                        "mj32", "mj42"),
                                                                arrows=1, free=FALSE, values=1),
 
                                                         # twin1 set slope factor loadings as age definition variables
                                                         mxPath(from='slope1', to=c("mj11", "mj21",    
                                                                                    "mj31","mj41"),
                                                                arrows=1, free=FALSE,
                                                                labels=c('data.age11','data.age21','data.age31','data.age41')),
 
                                                         # twin2 set slope factor loadings as age definition variables
                                                         mxPath(from='slope2', to=c("mj12", "mj22",    
                                                                                    "mj32", "mj42"),
                                                                arrows=1, free=FALSE,
                                                                labels=c('data.age12','data.age22','data.age32','data.age42')),
 
 
                                                         # twin1 estimate growth factor means, constraining to equality between twins
                                                         mxPath(from='one', to=c('intercept1','slope1'),
                                                                arrows=1, free=TRUE, values=c(1, .2), labels=c('intercept_mean','slope_mean')),
 
                                                         # twin1 estimate growth factor means, constraining to equality between twins
                                                         mxPath(from='one', to=c('intercept2','slope2'),
                                                                arrows=1, free=TRUE, values=c(1, .2), labels=c('intercept_mean','slope_mean')),
 
 
 
                                                         ### MZ Intercept ACE Specification
 
                                                         # variances of intercept ACE factors set to 1
                                                         IlatVariances <- mxPath( from=c("AI1","CI1","EI1","AI2","CI2","EI2"), arrows=2,
                                                                                  free=FALSE, values=1 ),
 
                                                        # means of intercept ACE factors set to 0
                                                         IlatMeans     <- mxPath( from="one", to=c("AI1","CI1","EI1","AI2","CI2","EI2"), arrows=1,
                                                                                  free=FALSE, values=0 ),
 
                                                         # A path coefficients for the Intercept (with slope cholesky cross paths)
                                                         IpathAT1    <- mxPath( from=c("AI1"), to=c("intercept1","slope1"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ai','sona') ),
 
                                                         # A path coefficients for the Intercept (with slope cholesky cross paths)
                                                         IpathAT2    <- mxPath( from=c("AI2"), to=c("intercept2","slope2"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ai','sona') ),
 
                                                         # C path coefficients for the Intercept (with slope cholesky cross paths)
                                                         IpathCT1    <- mxPath( from=c("CI1"), to=c("intercept1","slope1"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ci','sonc') ),
 
                                                         # C path coefficients for the Intercept (with slope cholesky cross paths)
                                                         IpathCT2    <- mxPath( from=c("CI2"), to=c("intercept2","slope2"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ci','sonc') ),
 
                                                         # E path coefficients for the Intercept (with slope cholesky cross paths)
                                                         IpathET1    <- mxPath( from=c("EI1"), to=c("intercept1","slope1"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ei','sone') ),
 
                                                         # E path coefficients for the Intercept (with slope cholesky cross paths)
                                                         IpathET2    <- mxPath( from=c("EI2"), to=c("intercept2","slope2"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ei','sone') ),
 
                                                         # covariance between C1 & C2
                                                         IcovC1C2      <- mxPath( from="CI1", to="CI2", arrows=2,
                                                                                  free=FALSE, values=1 ),
 
                                                         # covariance between A1 & A2 in MZ twins set to 1
                                                         IcovA1A2_MZ   <- mxPath( from="AI1", to="AI2", arrows=2,
                                                                                  free=FALSE, values=1 ),
 
                                                         ### MZ Slope ACE Specification
 
                                                         # variances of intercept ACE factors set to 1
                                                         SlatVariances <- mxPath( from=c("AS1","CS1","ES1","AS2","CS2","ES2"), arrows=2,
                                                                                  free=FALSE, values=1 ),
                                                         # means of intercept ACE factors set to 0
                                                         SlatMeans     <- mxPath( from="one", to=c("AS1","CS1","ES1","AS2","CS2","ES2"), arrows=1,
                                                                                  free=FALSE, values=0 ),
                                                         # path coefficients for intercept ACE factor of twin 1
                                                         SpathAceT1    <- mxPath( from=c("AS1","CS1","ES1"), to="slope1", arrows=1,
                                                                                  free=TRUE, values=.5,  label=c("as","cs","es") ),
                                                         # path coefficients for intercept ACE factor of twin 2
                                                         SpathAceT2    <- mxPath( from=c("AS2","CS2","ES2"), to="slope2", arrows=1,
                                                                                  free=TRUE, values=.5,  label=c("as","cs","es") ),
                                                         # covariance between C1 & C2
                                                         ScovC1C2      <- mxPath( from="CS1", to="CS2", arrows=2,
                                                                                  free=FALSE, values=1 ),
                                                         # covariance between A1 & A2 in MZ twins
                                                         ScovA1A2_MZ   <- mxPath( from="AS1", to="AS2", arrows=2,
                                                                                  free=FALSE, values=1 )
 
 
 
                                                 ), # close MZ Model
 
                                                 ### MZ Growth Model Specification
 
                                                 mxModel('dzmodel',type='RAM', mxData(observed=mj_twins_dz, type='raw'),
                                                         manifestVars=c("mj11", "mj21",    
                                                                        "mj31","mj41",
                                                                        "mj12", "mj22",    
                                                                        "mj32","mj42"),
                                                         latentVars=c(
                                                           # growth factors
                                                           'intercept1','slope1','intercept2','slope2',
                                                           # Intercept ACE factors for twin 1 and twin 2
                                                           "AI1","CI1","EI1","AI2","CI2","EI2",
                                                           # Slope ACE factors for twin 1 and twin 2
                                                           "AS1","CS1","ES1","AS2","CS2","ES2"),
 
 
                                                         # twin1 residual variance paths, freely estimated across time but constrained to equality between twins
                                                         mxPath(from=c("mj11", "mj21",    
                                                                       "mj31","mj41"),
                                                                arrows=2, free=TRUE, values=10, labels=c('rv1','rv2','rv3','rv4')),
 
                                                         # twin2 residual variance paths, freely estimated across time but constrained to equality between twins
                                                         mxPath(from=c("mj12", "mj22",    
                                                                       "mj32", "mj42"),
                                                                arrows=2, free=TRUE, values=10,labels=c('rv1','rv2','rv3','rv4')),
 
 
                                                         # twin1 intercept residual variance set to 0 (all consumed by ACE factors)
                                                         mxPath(from=c('intercept1'),arrows=2, connect='unique.pairs',free=F, labels=c('intercept_variance'),values=0),
 
                                                         # twin2 intercept residual variance set to 0 (all consumed by ACE factors)
                                                         mxPath(from=c('intercept2'),arrows=2, connect='unique.pairs',free=F, labels=c('intercept_variance'),values=0),
 
                                                         # twin1 slope residual variance set to 0 (all consumed by ACE factors)
                                                         mxPath(from=c('slope1'),arrows=2, connect='unique.pairs',free=F, labels=c('slope_variance'),values=0),
 
                                                         # twin2 slope residual variance set to 0 (all consumed by ACE factors)
                                                         mxPath(from=c('slope2'),arrows=2, connect='unique.pairs',free=F, labels=c('slope_variance'),values=0),
 
                                                         # twin1 set intercept factor loadings to 1
                                                         mxPath(from='intercept1', to=c("mj11", "mj21",    
                                                                                        "mj31","mj41"),
                                                                arrows=1, free=FALSE, values=1),
 
                                                         # twin2 set intercept factor loadings to 1
                                                         mxPath(from='intercept2', to=c("mj12", "mj22",    
                                                                                        "mj32", "mj42"),
                                                                arrows=1, free=FALSE, values=1),
 
                                                         # twin1 set slope factor loadings as age definition variables
                                                         mxPath(from='slope1', to=c("mj11", "mj21",    
                                                                                    "mj31","mj41"),
                                                                arrows=1, free=FALSE,
                                                                labels=c('data.age11','data.age21','data.age31','data.age41')),
 
                                                         # twin2 set slope factor loadings as age definition variables
                                                         mxPath(from='slope2', to=c("mj12", "mj22",    
                                                                                    "mj32", "mj42"),
                                                                arrows=1, free=FALSE,
                                                                labels=c('data.age12','data.age22','data.age32','data.age42')),
 
 
                                                         # twin1 estimate growth factor means constraining to equality between twins
                                                         mxPath(from='one', to=c('intercept1','slope1'),
                                                                arrows=1, free=TRUE, values=c(1, .2), labels=c('intercept_mean','slope_mean')),
 
                                                         # twin1 estimate growth factor means constraining to equality between twins
                                                         mxPath(from='one', to=c('intercept2','slope2'),
                                                                arrows=1, free=TRUE, values=c(1, .2), labels=c('intercept_mean','slope_mean')),
 
                                                         ### DZ Intercept ACE Specification
 
                                                         # variances of intercept ACE factors set to 1
                                                         IlatVariances <- mxPath( from=c("AI1","CI1","EI1","AI2","CI2","EI2"), arrows=2,
                                                                                  free=FALSE, values=1 ),
                                                         # means of intercept ACE factors set to 0
                                                         IlatMeans     <- mxPath( from="one", to=c("AI1","CI1","EI1","AI2","CI2","EI2"), arrows=1,
                                                                                  free=FALSE, values=0 ),
 
                                                         # A path coefficients for the Intercept (with slope Cholesky cross paths)
                                                         IpathAT1    <- mxPath( from=c("AI1"), to=c("intercept1","slope1"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ai','sona') ),
 
                                                         # A path coefficients for the Intercept (with slope Cholesky cross paths)
                                                         IpathAT2    <- mxPath( from=c("AI2"), to=c("intercept2","slope2"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ai','sona') ),
 
                                                         # C path coefficients for the Intercept (with slope Cholesky cross paths)
                                                         IpathCT1    <- mxPath( from=c("CI1"), to=c("intercept1","slope1"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ci','sonc') ),
 
                                                         # C path coefficients for the Intercept (with slope Cholesky cross paths)
                                                         IpathCT2    <- mxPath( from=c("CI2"), to=c("intercept2","slope2"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ci','sonc') ),
 
                                                         # E path coefficients for the Intercept (with slope Cholesky cross paths)
                                                         IpathET1    <- mxPath( from=c("EI1"), to=c("intercept1","slope1"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ei','sone') ),
 
                                                         # E path coefficients for the Intercept (with slope Cholesky cross paths)
                                                         IpathET2    <- mxPath( from=c("EI2"), to=c("intercept2","slope2"), arrows=1,
                                                                                free=TRUE, values=.5,  label=c('ei','sone') ),
 
                                                         # covariance between C1 & C2
                                                         IcovC1C2      <- mxPath( from="CI1", to="CI2", arrows=2,
                                                                                  free=FALSE, values=1 ),
                                                         # covariance between A1 & A2 in DZ twins
                                                         IcovA1A2_DZ   <- mxPath( from="AI1", to="AI2", arrows=2,
                                                                                  free=FALSE, values=.5 ),
 
                                                         ### DZ Slope ACE Specification
 
                                                         # variances of intercept ACE factors set to 1
                                                         SlatVariances <- mxPath( from=c("AS1","CS1","ES1","AS2","CS2","ES2"), arrows=2,
                                                                                  free=FALSE, values=1 ),
                                                         # means of intercept ACE factors set to 0
                                                         SlatMeans     <- mxPath( from="one", to=c("AS1","CS1","ES1","AS2","CS2","ES2"), arrows=1,
                                                                                  free=FALSE, values=0 ),
                                                         # path coefficients for intercept ACE factor of twin 1
                                                         SpathAceT1    <- mxPath( from=c("AS1","CS1","ES1"), to="slope1", arrows=1,
                                                                                  free=TRUE, values=.5,  label=c("as","cs","es") ),
                                                         # path coefficients for intercept ACE factor of twin 2
                                                         SpathAceT2    <- mxPath( from=c("AS2","CS2","ES2"), to="slope2", arrows=1,
                                                                                  free=TRUE, values=.5,  label=c("as","cs","es") ),
                                                         # covariance between C1 & C2
                                                         ScovC1C2      <- mxPath( from="CS1", to="CS2", arrows=2,
                                                                                  free=FALSE, values=1 ),
                                                         # covariance between A1 & A2 in MZ twins
                                                         ScovA1A2_DZ   <- mxPath( from="AS1", to="AS2", arrows=2,
                                                                                  free=FALSE, values=.5 )
 
 
 
 
 
                                                 ), # Close DZ  Model
 
 
 
                                                 mxAlgebra(mzmodel.objective + dzmodel.objective, name='mg_objective'), mxFitFunctionAlgebra('mg_objective')
 
 
 
 
 
) # Close Model
 
 
cholesky_ace_bivariate_multigroup_lgm.fit <- mxRun(cholesky_ace_bivariate_multigroup_lgm)
summary(cholesky_ace_bivariate_multigroup_lgm.fit)
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
As for the script itself, it

As for the script itself, it looks pretty good. My only concern is that it defines two latent C factors (1 per twin) each for the latent intercept and slope, and in the MZ model, two A factors each for the latent intercept and slope. Offhand, I am not sure how well OpenMx handles perfectly correlated latent variables in type="RAM" models. In the few times I've ever written a twin-analysis OpenMx script using RAM-path specification, I've simply replaced each pair of perfectly correlated latent variables with a single latent variable.

One suggestion: replace this,

mxAlgebra(mzmodel.objective + dzmodel.objective, name='mg_objective'), mxFitFunctionAlgebra('mg_objective')

, with this,

mxFitFunctionMultigroup(c("mzmodel","dzmodel"))

. It should make no difference for simply running the model, but there are a few useful functions in OpenMx that aren't compatible with MxModels using algebra fitfunctions.

ethibodeau's picture
Offline
Joined: 05/16/2018 - 17:05
great

Yes, I was wondering if OpenMx would bark at me about the perfectly correlated factors. Mplus does, it still runs the model. OpenMx, surprisingly did not bark. The -2(LL) between Mplus and OpenMx matched perfectly as well. Thanks for looking this over.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
Oh, nice! Thanks for giving

Oh, nice! Thanks for giving it a try.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
Age & sex?

BTW, you will probably also want to adjust for the main effects of age and sex.

ethibodeau's picture
Offline
Joined: 05/16/2018 - 17:05
sex yes, but age?

yes I wanted to start simple, then ill add sex as a covariate. However, why would I use age as a covariate when age is my time metric? I'm using age definition variables for the slope factor loadings.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
my mistake
However, why would I use age as a covariate when age is my time metric? I'm using age definition variables for the slope factor loadings.

Ah, good point. You're already adjusting for age.

ethibodeau's picture
Offline
Joined: 05/16/2018 - 17:05
ignore the weird code display with respect to the subsetting

Not sure why that looks weird, I think because of the use of $ signs. Just assume the data gets read in correctly :)

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
HTML

I just edited your post with the code to use the code HTML tags. code tags make the literal text display in a fixed-width font, without line-wrapping or table-breaking text, and with syntax highlighting if the page supports it. The fixed-width font and the R syntax highlighting both go a long way toward making the code readable...especially because the default format for these forums is full HTML with markdown support. Without the code tag, some of the symbols in your code--such as the #, $, and % get interpreted as markdown, which in this case is a hindrance to readability.

ProTip: use code tags, or at the very least, use the dropdown menu when writing a post to change its text format to filtered HTML.

ethibodeau's picture
Offline
Joined: 05/16/2018 - 17:05
thanks!

Good to know! That will be helpful!

Log in or register to post comments