# Decomposed ACE variances...now trying to run regressions

12 posts / 0 new
Offline
Joined: 02/05/2021 - 14:54
Decomposed ACE variances...now trying to run regressions

I am 100% new to twin data and OpenMx, and so far the tutorials and documentations have been superbly helpful!
I used this [genetic epi path specification annotated guide] (https://vipbg.vcu.edu/vipbg/OpenMx2/docs//OpenMx/latest/GeneticEpi_Path.html#ace-model-a-twin-analysis) to decompose the variances of ACE on a given phenotype (my DV).

Now, though, I need to run a regression and I can't figure out how to feed the ACE variances from above into the regressions. I assume I have to constrain the parameters somehow?

I have a categorical IV (an unshared environmental factor) ~ continuous Phenotype DV + covariates

And my research question is whether the E (unshared environmental factor) is significantly associated with the DV, after controlling for A and C confounders.

Is there any documentation I missed on this script?

Thank you so much for any help.

Offline
Joined: 01/24/2014 - 12:15

You should post your script as it currently stands.

Now, though, I need to run a regression and I can't figure out how to feed the ACE variances from above into the regressions. I assume I have to constrain the parameters somehow?

I have a categorical IV (an unshared environmental factor) ~ continuous Phenotype DV + covariates

And my research question is whether the E (unshared environmental factor) is significantly associated with the DV, after controlling for A and C confounders.

I don't follow this part. In particular, I think you mean "continuous Phenotype DV ~ categorical IV (an unshared environmental factor) + covariates"? Also, what makes you so sure the IV is nonshared-environmental? How many categories does it have, and are they ordered?

Offline
Joined: 02/05/2021 - 14:54
Script and clarification

Script is below, but so far it just includes the ACE variance breakdown. I do not understand how to write the regression script (which is where I am stuck). To clarify your

I don't follow this part. In particular, I think you mean "continuous Phenotype DV ~ categorical IV (an unshared environmental factor) + covariates"? Also, what makes you so sure the IV is nonshared-environmental? How many categories does it have, and are they ordered?

:
1. You are right, I meant that the continuous DV ~ categorical IV + covariates
2. The IV being non-shared is the crux of my research question. My research question asks if going outdoors into nature affects mental health. So, the behavior of going outside into nature (which may have some genetic influences, but is also the observed difference in unshared environmental, i.e. the behavior) is hypothesized to be associated with improved mental health (the phenotype DV).
3. The IV has 5 categories, ordered.

Any help feeding my ACE model summary into the regression model is needed.

# Select Variables for Analysis
selVars   <- c('cohenpss_score_A','cohenpss_score_B') #Selecting these columns for analysis
aceVars   <- c("A1","C1","E1","A2","C2","E2") #creating latent variables

# Select Data for Analysis: Creating two subsets for MZ and DZ to generate descriptive stats for each group type
mzData    <- subset(df, zyg_A==1, selVars)
dzData    <- subset(df, zyg_A==0, selVars)

# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE) #mean ~11
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete") #cov ~45
cov(dzData,use="complete")

# Path objects for Multiple Groups
manifestVars=selVars #Manifest = the Stress phenotype
latentVars=aceVars

# variances of latent variables
latVariances <- mxPath( from=aceVars, arrows=2,
free=FALSE, values=1 )

# means of latent variables... AND SPECIFYING SINGLE=HEADED ARROWS FROM TRIANGLE (WHICH IS FIXED AT 1) TO EACH OF THE LATENT VARIABLES, FIXED AT 0

latMeans     <- mxPath( from="one", to=aceVars, arrows=1,
free=FALSE, values=0 )

# means of observed variables

obsMeans     <- mxPath( from="one", to=selVars, arrows=1,
free=TRUE, values=10, labels="mean" ) #The mean cohen_pss_score is around 11

# path coefficients for twin 1
pathAceT1    <- mxPath( from=c("A1","C1","E1"), to="cohenpss_score_A", arrows=1,
free=TRUE, values=.5,  label=c("a","c","e") )
# path coefficients for twin 2
pathAceT2    <- mxPath( from=c("A2","C2","E2"), to="cohenpss_score_B", arrows=1,
free=TRUE, values=.5,  label=c("a","c","e") )

#The common environmental factors are the same for both twins (by definition), so setting the fixed correlation of 1 for C1 and C2

# covariance between C1 & C2
covC1C2      <- mxPath( from="C1", to="C2", arrows=2,
free=FALSE, values=1 )

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

# Data objects for Multiple Groups
dataMZ       <- mxData( observed=mzData, type="raw" )
dataDZ       <- mxData( observed=dzData, type="raw" )

# Combine Groups
paths        <- list( latVariances, latMeans, obsMeans,
pathAceT1, pathAceT2, covC1C2 )
modelMZ      <- mxModel(model="MZ", type="RAM", manifestVars=selVars,
latentVars=aceVars, paths, covA1A2_MZ, dataMZ )
modelDZ      <- mxModel(model="DZ", type="RAM", manifestVars=selVars,

obj          <- mxFitFunctionMultigroup(c("MZ", "DZ"))
modelACE     <- mxModel(model="ACE", modelMZ, modelDZ, obj )

# Run Model
fitACE       <- mxRun(modelACE)
sumACE       <- summary(fitACE)

Offline
Joined: 01/24/2014 - 12:15
Let me be sure I understand.

Let me be sure I understand. Your present script biometrically decomposes the variance in your continuous DV. You want to also regress the DV onto your ordinal IV and some nuisance covariates. You further want to biometrically decompose the ordinal IV. Is all that correct?

Offline
Joined: 02/05/2021 - 14:54
Correct

All correct!

Offline
Joined: 01/24/2014 - 12:15
OK, I have a quite clear idea

OK, I have a quite clear idea of what your model would look like. I'm still not sure exactly what to change about your script.

Are your covariates age and sex, perhaps?

Offline
Joined: 01/09/2020 - 14:36

How about a simple bivariate genetic model? You can set this up as a Cholesky with covariates. There are many scripts here, but if that’s all you want to do, you might want to consider umx. The function umxACE will take care of this. If you go down the openmx road, be sure to to check out joint continuous ordinal examples with definition variables.

Subsequently you can observe the e cross path and do some testing with it. I.e. does model fit decrease significantly if you fix the path to zero?

Offline
Joined: 03/01/2013 - 14:09

I would recommend umxACEv instead of umxACE. The multivariate likelihood ratio tests for umxACE (the Cholesky version) are incorrect because of the non-negative definite constraints that the Cholesky implicitly forces. See Verhulst et al 2019 for more info on this issue.

Offline
Joined: 01/24/2014 - 12:15
try this

OK, I'll assume that the two "nuisance" covariates are age and sex. In that case, since the biometrical variance decompositions of both the DV and the IV are of interest, the IV will need to be adjusted for age and sex as well.

Give the attached script a try. Obviously, I can't actually test it, since I don't have your dataset.

Note that there are alternate ways in which the location and scale of the latent liability underlying the ordinal IV could be identified.

Edit: my mistake--you'll need to use latentVars=c(aceVars,dummyVars) when you construct modelMZ and modelDZ.

File attachments:
Offline
Joined: 02/05/2021 - 14:54
covariance not positive definite

Thank you very much for your helpful script. Any suggestions on the best way to fix a covariance matrix that is not a definite positive?

 fitACE       <- mxRun(modelACE)
Running ACE with 16 parameters
constraint-adjusted standard errors could not be calculated because the coefficient matrix of the quadratic form was uninvertibleIn model 'ACE' Optimizer returned a non-zero status code 10. Starting values are not feasible. Consider mxTryHard()Error: The job for model 'ACE' exited abnormally with the error message: fit is not finite (0: The continuous part of the model implied covariance (loc3) is not positive definite in data 'DZ.data' row 89. Detail:
covariance =  matrix(c(    # 6x6
0.7575, 0.37875,   0,   0,   0,   0
, 0.37875, 0.7575,   0,   0,   0,   0
,   0,   0,   0,   0,   0,   0
,   0,   0,   0,   0,   0,   0
,   0,   0,   0,   0,   0,   0
,   0,   0,   0,   0,   0,   0), byrow=TRUE, nrow=6, ncol=6)

1: The continuous part of the model implied covariance (loc3) is not positive definite in data 'DZ.data' row 82. Detail:
covariance =  matrix(c(    # 6x6
0.7575, 0.37875,   0,   0,   0,   0
, 0.37875, 0.7575,   0,   0,   0,   0
,   0,   0,   0,   0,   0,   0
,   0,   0,   0,   0,   0,   0
,   0,   0,   0,   0,   0,   0
,   0,   0,   0,   0,   0,   0), byrow=TRUE, nrow=6, ncol=6) `
Offline
Joined: 01/24/2014 - 12:15
6x6?

How is the model-expected covariance matrix 6x6? I thought there were only two endogenous variables per twin, so it should be 4x4.

At any rate, it appears some paths that should be free are instead fixed to zero.

Offline
Joined: 02/05/2021 - 14:54
Tried collapsing categories

For the covariance not definite positive issue, I tried collapsing some of the categories but that didn't work