You are here

ACE decomposition of growth factors code

19 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!

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

Obviously I want to estimate standardized variance components of my ACE factors for the intercept and slope. I see on the OpenMx user guide example for genetic epidemiology path specification some code to do something like that:

# Generate & Print Output
# additive genetic variance, a^2
A  <- mxEval(a*a, fitACE)
# shared environmental variance, c^2
C  <- mxEval(c*c, fitACE)
# unique environmental variance, e^2
E  <- mxEval(e*e, fitACE)
# total variance
V  <- (A+C+E)
# standardized A
a2 <- A/V
# standardized C
c2 <- C/V
# standardized E
e2 <- E/V

Is there anyway to estimate those as new parameters in my model instead, so I can get standard errors and CIs? I would love to take advantage of OpenMx's likelihood based CIs for use with standardized variance components. I already obtained CIs for my ACE factor loadings, however, it's not straightforward how I use those to obtain CIs for standardized variance components (for the upper and lower bound estimates that is). For example, if I wanted to know the 95% lower bound standardized variance component of A for the intercept, would I do something like this:

lower bounded A-intercept factor loading^2/(lower bounded A-intercept factor loading^2 + mid-bound C-intercept factor loading^2 + mid-bound E-intercept factor loading^2)

Thanks! BTW those likelihood based confidence intervals are now my new favorite feature of OpenMx!

side note- In my model above, I forgot to add ACE decomps of the residual variances of the growth factor indicators! Oops.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
some tips
Is there anyway to estimate those as new parameters in my model instead, so I can get standard errors and CIs?

Yes, there is. BTW, you're seeing firsthand the technical reason why twin analyses aren't typically taught using RAM-path specification: no matter how a twin model is parameterized, there will usually be quantities of interest that are not the free parameters themselves, but functions of the free parameters ("algebras", informally speaking). You are making the correct choice in calculating confidence intervals for variance components rather than one-headed Cholesky paths, since (at least without setting lbounds) some Cholesky paths are sign-indeterminate.

Anyhow, what you'll need to do is assign labels to the paths that are of interest to you, using the labels argument to mxPath(). You can then make MxMatrices containing elements that share labels with paths, or simply use the labels themselves (without quotes) in MxAlgebras, as though they were scalars (1x1 matrices). Then, create MxAlgebras that calculate the component covariance matrices and total phenotypic covariance matrices, and from them, the "standardized" component covariance matrices. Finally, request confidence intervals for the algebras of interest by including an mxCI() statement in the appropriate place in your MxModel, and don't forget to pass argument intervals=TRUE to mxRun().

Thanks! BTW those likelihood based confidence intervals are now my new favorite feature of OpenMx!

Yes, they are nice, aren't they!? They are recommended over standard errors for inferential purposes. But if for some reason you want standard errors, you can request standard errors on algebras with mxSE(). The smart thing to do is use imxRobustSE() (which is kind of slow) on your MxModel, with argument details=TRUE, and store its output (a list) in an object. Then, provide the 'cov' element of that list as the value for argument cov to mxSE(). That way, mxSE() will be using the robust rather than naive sampling-covariance matrix to get delta-method standard errors for the algebras of interest.

If you decide to re-parameterize your model, you may then find mxStandardizeRAMpaths() useful. It can also accept a robust sampling-covariance matrix as an argument, but it only reports results for paths as they are specified in the given MxModel, i.e. path coefficients as they are and not functions of path coefficients.

For example, if I wanted to know the 95% lower bound standardized variance component of A for the intercept, would I do something like this:

lower bounded A-intercept factor loading^2/(lower bounded A-intercept factor loading^2 + mid-bound C-intercept factor loading^2 + mid-bound E-intercept factor loading^2)

I'm not quite sure what that's supposed to mean...?

side note- In my model above, I forgot to add ACE decomps of the residual variances of the growth factor indicators! Oops.

Whoa, I'm surprised I overlooked that, too.

ethibodeau's picture
Offline
Joined: 05/16/2018 - 17:05
speaking of residual variances

I'm trying to add the ACE decomps of the residual variances right now. Is there anything special I do with the indicator means? Normally in a growth model the indicator means are set to zero. Is that still the case when I decompose their residual variances? I'm slightly confused as to why it is necessary to decompose the residual variances of the indicators. Thanks!

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
means model doesn't change
I'm trying to add the ACE decomps of the residual variances right now. Is there anything special I do with the indicator means? Normally in a growth model the indicator means are set to zero. Is that still the case when I decompose their residual variances?

Yes, the "direct" means (the one-headed paths from "one") of the indicators remain zeroes. You're just adding new random effects, which doesn't call for any change the model for the fixed effects.

I'm slightly confused as to why it is necessary to decompose the residual variances of the indicators.

I'm not sure I'd say it's "necessary," but it is often done, because there might be familial time-specific variance in the trait. You'll find out soon enough from comparing AICs whether or not the additional variance decomposition improves your model's efficiency or not.