You are here

Home › Forums › OpenSEM Forums › Behavioral Genetics Models

Causal Common Contingent Modeling Code Help

3 posts / 0 new
Log in or register to post comments
Last post
Mon, 10/10/2022 - 05:36
#1
szellers's picture
szellers
Offline
Joined: 04/04/2018 - 14:52
Causal Common Contingent Modeling Code Help
AttachmentSize
Binary Data Bivariate_RegularQuantity_75.R7.35 KB

Hi all,

I'm working on a model of smoking quantity given initiation of regular smoking along the lines of Bares, Kendler, and Maes 2015 among others. I'm having a bit of trouble figuring out how to move from a more straightforward bivariate model to the CCC model by implementing the estimated path from initiation to quantity. For reference, this is designated b in Figure 1 at https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4631696/

I have a few questions to help me code this properly. I've attached a script for reference that I'm using as a starting point for the CCC, its a standard bivariate model of binary regular initiation and continuous cigarette quantity, with a covariate of sex on the means. I have a few questions that I think will help me figure out the code, so I would appreciate answers to these questions and/or code alterations that would properly code the CCC model.

First question, I've regressed a covariate onto phenotypic means in fore using definition variables, is this the same process I would use to code the b path? Secondly, I've only used covariates previously where they do not vary within a twin pair (i.e., effect of sex on mean phenotype in twin models of same sex twins only). Twins can vary on smoking initiation, so how would this be coded? Lastly, there should not be cross paths between the ACE components for initiation and quantity, because the effect of initiation on quantity is through only path b, is that correct?

Thanks for your help,

Stephanie

Top
Wed, 10/12/2022 - 09:18
szellers's picture
szellers
Offline
Joined: 04/04/2018 - 14:52
Been playing around with script trying to code this

Hi again,

I've continued trying to code the CCC model and have an attempt here. It does run and produces somewhat reasonable output, but I wanted to run it past folks here to confirm that it is actually doing what I think it is doing! The estimate for path b is about 1. I would appreciate if someone could take a look and confirm if this is how the CCC should be parameterized for a binary and continuous variable. Thanks :)

File attachments: 
Binary Data CCC_RegularQuantity_75.R
Top
Mon, 11/21/2022 - 08:56
szellers's picture
szellers
Offline
Joined: 04/04/2018 - 14:52
Unsure if I have a problem in code or if this is to be expected

Hi again,

I have continued working on the CCC model code above. I noticed that if I run a univariate twin model on my stage 1 variable (here initiation of regular smoking) and compare the results (variance decomposition, threshold, age covariate effect) to that in the CCC model where the stage 2 variable is regressed on stage 1 via path B, the results change dramatically. Every parameter related to the stage 1 variable differs from the univariate results of just the stage 1 variable.

Is this something to be expected? It seems odd to me, but I cannot locate an issue in my code or interpret why this might be happening. I've tried putting lower bounds of 0 on the A, C, E matrices in the CCC model, as there is a large negative C component for the stage 1 variable when modeled without the bound. This reduced the A and C estimates for the stage 1 variable to virtually 0, resulting in an E only model. Modeling the stage 1 variable alone would not result in an E only model though. I'd appreciate any insight into this, perhaps the CCC model isn't appropriate, or I've made some mistake in my code? Or perhaps this is typical for CCC models, where the stage 1 variable estimates are different from the univariate results?

I've pasted here the output from the univariate model of stage 1 variable (top) vs. the CCC model output for reference (bottom)

Univariate Output
free parameters:
name matrix row col Estimate Std.Error A
1 threshold MZ.threG 1 RegularCigaretteF0 -2.04738347 0.63550967
2 a MZ.A 1 1 0.17848167 0.27327289
3 c MZ.C 1 1 0.35317742 0.20017502
4 e MZ.E 1 1 0.46834091 0.10745419
5 Agebeta MZ.regCoefAge 1 1 -0.02456682 0.01007581

CCC Output
free parameters:
name matrix row col Estimate Std.Error A
1 meanCig MZ.meanG 2 1 4.54581481 0.8849285663
2 AgebetaDist MZ.regCoefAge 1 1 -0.01001620 NA !
3 AgebetaCig MZ.regCoefAge 2 1 -0.01903446 0.0140772730
4 threshold MZ.threG 1 RegularCigaretteF0 -0.10372887 0.0537373169
5 a1 MZ.A 1 1 0.26592143 0.0001461964 !
6 a3 MZ.A 2 2 0.92208906 0.5904506950
7 c1 MZ.C 1 1 -1.13296070 0.0002682202 !
8 c3 MZ.C 2 2 0.18886662 0.4700930440
9 e1 MZ.E 1 1 1.86703928 NA !
10 e3 MZ.E 2 2 1.29092373 0.2065502970
11 B MZ.regCoefEver 1 2 0.81391368 0.0379732408

CCC Code

rm(list=ls())
library(OpenMx)
library(readr)
library(haven)
LongitudinalTobData <- read_csv("~/Desktop/TobaccoCohort/Data/LongitudinalTobData.csv")
#LongitudinalTobData <- read_csv("LongitudinalTobData.csv")
 
OlderCohort <- subset(LongitudinalTobData, TwoCohorts==1, select=c("person_nb", "family_nb", "birthyear", "sex", "zyg", "TenYearBands",  "TwoCohorts", "Twin", 
                                                                   "Education75", "age75", "EverCigarette75", "RegularCigarette75", "CurrentCigarette75", "FormerQuantity75", 
                                                                   "CurrentQuantity75", "EverCigarPipe75", "RegularCigarPipe75", 
                                                                   "BeerQuant75", "WineQuant75", "SpiritsQuant75", "Binge75", "Quantity75"))
 
colnames(OlderCohort) <- c("person_nb", "family_nb", "birthyear", "sex", "zyg", "TenYearBands",  "TwoCohorts", "Twin", 
                           "Education", "age", "EverCigarette", "RegularCigarette", "CurrentCigarette", "FormerQuantity", 
                           "CurrentQuantity", "EverCigarPipe", "RegularCigarPipe", 
                           "BeerQuant", "WineQuant", "SpiritsQuant", "Binge", "Quantity")
 
 
 
YoungerCohort <- subset(LongitudinalTobData, TwoCohorts==2, select=c("person_nb", "family_nb", "birthyear", "sex", "zyg", "TenYearBands",  "TwoCohorts", "Twin", 
                                                                     "Education2011", "age2011", "EverCigarette2011", "RegularCigarette2011", "CurrentCigarette2011", "FormerQuantity2011", 
                                                                     "CurrentQuantity2011", "EverCigarPipe2011", "RegularCigarPipe2011", 
                                                                     "BeerQuant2011", "WineQuant2011", "SpiritsQuant2011", "Binge2011", "Quantity2011"))
 
colnames(YoungerCohort) <- c("person_nb", "family_nb", "birthyear", "sex", "zyg", "TenYearBands",  "TwoCohorts", "Twin", 
                             "Education", "age", "EverCigarette", "RegularCigarette", "CurrentCigarette", "FormerQuantity", 
                             "CurrentQuantity", "EverCigarPipe", "RegularCigarPipe", 
                             "BeerQuant", "WineQuant", "SpiritsQuant", "Binge", "Quantity")
 
LongitudinalTobDataCohorts <- as.data.frame(rbind(OlderCohort, YoungerCohort))
 
LongitudinalTobDataCohorts$RegularCigarette[LongitudinalTobDataCohorts$EverCigarette==1] <- 0
LongitudinalTobDataCohorts$RegularCigarette[LongitudinalTobDataCohorts$RegularCigarette==1] <- 0
LongitudinalTobDataCohorts$RegularCigarette[LongitudinalTobDataCohorts$RegularCigarette==2] <- 1
 
LongitudinalTobDataCohorts$RegularCigaretteF <- as.factor(LongitudinalTobDataCohorts$RegularCigarette)
LongitudinalTobDataCohorts$RegularCigaretteF <- ordered(LongitudinalTobDataCohorts$RegularCigarette)
 
LongitudinalTobDataWide <- reshape(LongitudinalTobDataCohorts, v.names=c("person_nb", "Education", "age", "EverCigarette", "RegularCigarette", 
                                                                         "CurrentCigarette", "FormerQuantity", "CurrentQuantity", "EverCigarPipe", 
                                                                         "RegularCigarPipe", "BeerQuant", "WineQuant", "SpiritsQuant", 
                                                                         "Binge", "Quantity", "RegularCigaretteF"),
                                   timevar = 'Twin', idvar=c('family_nb'), direction="wide", sep="")
 
 
LongitudinalTobDataWide$ageavg <- (LongitudinalTobDataWide$age0 + LongitudinalTobDataWide$age1)/2
LongitudinalTobDataWide$agetwinpair <- ifelse(is.na(LongitudinalTobDataWide$age0), LongitudinalTobDataWide$age1, LongitudinalTobDataWide$ageavg) 
LongitudinalTobDataWide$agetwinpair <- ifelse(is.na(LongitudinalTobDataWide$age1), LongitudinalTobDataWide$age0, LongitudinalTobDataWide$ageavg) 
LongitudinalTobDataWide$agetwinpair[is.na(LongitudinalTobDataWide$agetwinpair)] <- mean(LongitudinalTobDataWide$agetwinpair, na.rm=T) 
 
 
LongitudinalTobDataWide <- subset(LongitudinalTobDataWide, !is.na(RegularCigarette0) & !is.na(RegularCigarette1))
 
# Select variables
varsc <- c("Quantity")
nvc <- 1
ntvc <- nvc*2
conVars <- paste(varsc, c(rep(0,nvc), rep(1,nvc)), sep="")
 
nth <- 1
varso <- c("RegularCigaretteF")
nvo <- 1
ntvo <- nvo*2
ordVars <- paste(varso, c(rep(0,nvo), rep(1,nvo)), sep="")
 
oc <- c(1,0) #1 ordinal 0 continuous
 
Vars <- c("RegularCigaretteF", "Quantity")
nv <- nvc+nvo
ntv <- nv*2
selVars <- paste(Vars, c(rep(0,nv), rep(1,nv)), sep="")
 
DefVars <- c("agetwinpair")
 
frMV <- c(F, T)
 
 
# Select Data for Analysis
MzDataOlderMale    <- subset(LongitudinalTobDataWide, zyg==1 & TwoCohorts==1 & sex==1, select=c(selVars, DefVars))
DzDataOlderMale    <- subset(LongitudinalTobDataWide, zyg==2 & TwoCohorts==1 & sex==1, select=c(selVars, DefVars))
 
MzDataOlderFemale    <- subset(LongitudinalTobDataWide, zyg==1 & TwoCohorts==1 & sex==2, select=c(selVars, DefVars))
DzDataOlderFemale    <- subset(LongitudinalTobDataWide, zyg==2 & TwoCohorts==1 & sex==2, select=c(selVars, DefVars))
 
MzDataYoungerMale    <- subset(LongitudinalTobDataWide, zyg==1 & TwoCohorts==2 & sex==1, select=c(selVars, DefVars))
DzDataYoungerMale    <- subset(LongitudinalTobDataWide, zyg==2 & TwoCohorts==2 & sex==1, select=c(selVars, DefVars))
 
MzDataYoungerFemale    <- subset(LongitudinalTobDataWide, zyg==1 & TwoCohorts==2 & sex==2, select=c(selVars, DefVars))
DzDataYoungerFemale    <- subset(LongitudinalTobDataWide, zyg==2 & TwoCohorts==2 & sex==2, select=c(selVars, DefVars))
 
# Putting MZ/DZ data into OpenMx objects
DataMZOlderMale    <- mxData( observed=MzDataOlderMale, type="raw" )
DataDZOlderMale    <- mxData( observed=DzDataOlderMale, type="raw" )
 
DataMZOlderFemale    <- mxData( observed=MzDataOlderFemale, type="raw" )
DataDZOlderFemale    <- mxData( observed=DzDataOlderFemale, type="raw" )
 
DataMZYoungerMale    <- mxData( observed=MzDataYoungerMale, type="raw" )
DataDZYoungerMale    <- mxData( observed=DzDataYoungerMale, type="raw" )
 
DataMZYoungerFemale    <- mxData( observed=MzDataYoungerFemale, type="raw" )
DataDZYoungerFemale    <- mxData( observed=DzDataYoungerFemale, type="raw" )
 
 
 
 
#######################################################################################
 
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=ntv, ncol=1, free=c(F, T, F, T), values=c(0, 4, 0, 4), labels=c("meanDist","meanCig"), name="meanG" )
 
AgeCov <- mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, labels="data.agetwinpair", name="AgeCov")
regCoefAge <- mxMatrix(type="Full", nrow=ntv, ncol=1, free=TRUE, values=.2, labels=c("AgebetaDist", "AgebetaCig", "AgebetaDist", "AgebetaCig"), name="regCoefAge")
 
EverCov <- mxMatrix(type="Full", nrow=ntv, ncol=1, free=FALSE, labels=c("data.RegularCigaretteF0", "data.RegularCigaretteF1"), name="EverCov")
regCoefEver <- mxMatrix(type="Full", nrow=1, ncol=ntv, free=c(F, T, F, T), values=c(0, .2, 0, .2), labels=c("NA","B", "NA","B"), name="regCoefEver")
 
expMean <- mxAlgebra(expression=t(meanG+(regCoefAge%*%AgeCov)+(regCoefEver%*%EverCov)), name="expMean")
#expMean <- mxAlgebra(expression=t(meanG+(regCoefAge%*%AgeCov)), name="expMean")
 
threG <- mxMatrix(type="Full", nrow=nth, ncol=ntvo, free=T, values=.01, labels=c("threshold", "threshold"), name="threG")
 
# Create Matrices for Path Coefficients
A <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=c(T, F, T), values=c(.5, 0,.5), labels=c("a1", "a2","a3"),  name="A" )
C <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=c(T, F, T), values=c(.5, 0,.5), labels=c("c1", "c2", "c3"),  name="C" )
E <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=c(T, F, T), values=c(.5, 0,.5), labels=c("e1", "e2", "e3"),  name="E" )
 
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= A+C+E, name="V" )
covMZ <- mxAlgebra( expression= A+C, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%A+ C, name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" )
 
#constrain variance of binary variables
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="matUnv" )
matOc <- mxMatrix( type="Full", nrow=1, ncol=nv, free=FALSE, values=oc, name="matOc" )
var1 <- mxConstraint( expression=diag2vec(matOc%&%V)==matUnv, name="Var1" )
 
# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars)
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars)
funML <- mxFitFunctionML()
 
# Create Model Objects for Multiple Groups
#
pars <- list( meanG, AgeCov, regCoefAge,  expMean, threG, A, C, E, covP, var1, matUnv, matOc, EverCov, regCoefEver )
modelMZ <- mxModel( pars, covMZ, expCovMZ, DataMZOlderMale, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, covDZ, expCovDZ, DataDZOlderMale, expDZ, funML, name="DZ" )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
ciACE <- mxCI( c("a1", "c1", "e1", "a3", "c3", "e3", "B") )
 
# Build Model with Confidence Intervals
modelACE <- mxModel( "oneACEc", modelMZ, modelDZ, multi, ciACE )
 
 
###################
### Test Models ###
###################
 
# Run ACE Model
#fitACE <- mxTryHardOrdinal( modelACE, intervals=F)
#fitACE <- mxOption(modelACE, 'mvnRelEps', mxOption(modelACE, 'mvnRelEps')/5)
fitACE <- mxTryHard(modelACE)
summary( fitACE)
Top

Copyright © 2007-2024 The OpenMx Project