# Author: Frühling Rijsdijk # MODEL: MV quantitative & qualitative sex-dif Sex-Limitation script (ACE Correlated Factors model) # Correlation Approach to ensure that the order of the variables does NOT affect the ability # of the model to account for the DZOS data. # Ref: Neale et al., Multivariate genetic analysis of sex-lim and G x E interaction, Twin Research & Human Genetics, 2006 # Restrictions: Assumes means and variances can be equated across birth order within zygosity groups # ------------------------------------------------------------------------------------------------------------------------------------------------------ #setwd("C:/Users/Karen Burton/Dropbox/Analyses/DASS") require(psych) require(OpenMx) #source("GenEpiHelperFunctions.R") #source("RFunctions.R") DASS <- read.csv("~/mx/openmx/problems/DASS_Data_Dummy_1.csv", header=TRUE) # Read in data file # ----------------------------------------------------------------------- nvar <- 3 # number of variables tnvar <- 2*nvar # number of variables*max family size nlower <-tnvar*(tnvar+1)/2 # number of free elements in a lower matrix tnvar*tnvar mzData <- as.data.frame(subset(DASS, Zygosity==0 , names(DASS))) names(mzData) dzData <- as.data.frame(subset(DASS, Zygosity==1 , names(DASS))) mzmData <- as.data.frame (subset (mzData, Sex ==0, names(DASS))) mzfData <- as.data.frame (subset (mzData, Sex ==1, names(DASS))) dzmData <- as.data.frame (subset (dzData, Sex ==0, names(DASS))) dzfData <- as.data.frame (subset (dzData, Sex ==1, names(DASS))) selVars<-c("DepressALn","AnxietyALn","StressALn","DepressBLn","AnxietyBLn","StressBLn") lMeans<-c("Depress","Anxiety","Stress", "Depress","Anxiety","Stress") stM<-c(mean(mzData$DepressALn, na.rm=T)+ sd(mzData$DepressALn, na.rm=T), mean(mzData$AnxietyALn, na.rm=T)+ sd(mzData$AnxietyALn, na.rm=T), mean(mzData$StressALn, na.rm=T)+ sd(mzData$StressALn, na.rm=T)) stMeans<-cbind(stM,stM) stMDepress <-c(mean(mzData$DepressALn, na.rm=T)+ sd(mzData$DepressALn, na.rm=T)) stMAnxiety <-c( mean(mzData$AnxietyALn, na.rm=T)+ sd(mzData$AnxietyALn, na.rm=T)) stMStress <-c( mean(mzData$StressALn, na.rm=T)+ sd(mzData$StressALn, na.rm=T)) stVar<- c(var(mzData$DepressALn, na.rm=T)+ sd(mzData$DepressALn, na.rm=T), var(mzData$AnxietyALn, na.rm=T)+ sd(mzData$AnxietyALn, na.rm=T), var(mzData$StressALn, na.rm=T)+ sd(mzData$StressALn, na.rm=T)) stA<-sqrt(stVar*.6) stC<-sqrt(stVar*.2) stE<-sqrt(stVar*.5) nf <-1 time<-10000 # 1 Nonscalar Sex Limitation ##################################################################################### # Quantitative sex Diff & Qualitative sex Diff for A # In this first model there are male and female paths, plus male and female Ra, Rc and Re between variables # Male-Female correlations in DZO group between A factors FREE # ---------------------------------------------------------------------------------------------------------------- ACENonScSLCoraModel <- mxModel("ACENonScSLCora", mxModel("ACE", # Matrices a, c, and e to store a, c, and e path coefficients mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=stA, label=c("am11","am22", "am33"), name="am" ), mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=FALSE, values=0, label=c("cm11","cm22", "cm33"), name="cm" ), mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=stE, label=c("em11","em22", "em33"), name="em" ), mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=stA, label=c("af11","af22", "af33"), name="af" ), mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=FALSE, values=0, label=c("cf11","cf22", "cf33"), name="cf" ), mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=stE, label=c("ef11","ef22", "ef33"), name="ef" ), mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("ram21", "ram31", "ram32"), lbound=-.9999, ubound=.9999, name="Ram" ), mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=FALSE, values=0, label=c("rcm21", "rcm31", "rcm32"), lbound=-.9999, ubound=.9999, name="Rcm" ), mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("rem21", "rem31", "rem32"), lbound=-.9999, ubound=.9999, name="Rem" ), mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("raf21", "raf31", "raf32"), lbound=-.9999, ubound=.9999, name="Raf" ), mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=FALSE, values=0, label=c("rcf21", "rcf31", "rcf32"), lbound=-.9999, ubound=.9999, name="Rcf" ), mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("ref21", "ref31", "ref32"), lbound=-.9999, ubound=.9999, name="Ref" ), # Matrices A, C, and E compute variance components mxAlgebra( expression=am %*% (Ram) %*% t(am), name="Am" ), mxAlgebra( expression=cm %*% (Rcm) %*% t(cm), name="Cm" ), mxAlgebra( expression=em %*% (Rem) %*% t(em), name="Em" ), mxAlgebra( expression=af %*% (Raf) %*% t(af), name="Af" ), mxAlgebra( expression=cf %*% (Rcf) %*% t(cf), name="Cf" ), mxAlgebra( expression=ef %*% (Ref) %*% t(ef), name="Ef" ), # Algebra to compute total variances and standard deviations (diagonal only) mxAlgebra( expression=Am+Cm+Em, name="Vm" ), mxAlgebra( expression=Af+Cf+Ef, name="Vf" ), mxMatrix( type="Iden", nrow=nvar, ncol=nvar, name="I"), mxAlgebra( expression=solve(sqrt(I*Vm)), name="iSDm"), mxAlgebra( expression=solve(sqrt(I*Vf)), name="iSDf"), # Matrix & Algebra for expected means vector mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= stM, name="Mm" ), mxAlgebra( cbind(Mm,Mm), name="expMeanm"), mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= stM, name="Mf" ), mxAlgebra( cbind(Mf,Mf), name="expMeanf"), #matrix for the regresion coefficent -education mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values=-0.030, label=c("bedu_Depress","bedu_Anxiety","bedu_Stress"), name="b_edu" ), #matrix for the regresion coefficent - Age mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values=-0.015, label=c("bage_Depress","bage_Anxiety","bage_Stress"), name="b_age" ), # Algebra for expected variance/covariance matrix in MZ mxAlgebra( rbind ( cbind(Am+Cm+Em , Am+Cm), cbind(Am+Cm , Am+Cm+Em)), name="expCovMZm" ), mxAlgebra( rbind ( cbind(Af+Cf+Ef , Af+Cf), cbind(Af+Cf , Af+Cf+Ef)), name="expCovMZf" ), # Algebra for expected variance/covariance matrix in DZ mxAlgebra( rbind ( cbind(Am+Cm+Em , 0.5%x%Am+Cm), cbind(0.5%x%Am+Cm, Am+Cm+Em)), name="expCovDZm" ), mxAlgebra( rbind ( cbind(Af+Cf+Ef , 0.5%x%(Af) + Cf ), cbind(0.5%x%(Af) + Cf , Af+Cf+Ef)), name="expCovDZf" ), # Compute Variance Components #mxAlgebra( expression=rbind(cbind(Am/Vm,Cm/Vm,Em/Vm),cbind(Af/Vf,Cf/Vf,Ef/Vf)), name="VarComp"), mxModel("MZm", #matrix for the observed age mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.AgeA", "data.AgeB"), name="obs_age" ), #matrix for the observed Education mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.EducationA", "data.EducationB"), name="obs_edu" ), #Algebra to add covariates to mean # mxAlgebra( expression= ACE.expMeanm + ACE.b_age %x% obs_age + ACE.b_edu %x% obs_edu , name="expMeanmMZ" ), mxAlgebra( expression= ACE.expMeanm , name="expMeanmMZ" ), mxData( observed=mzmData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovMZm", means="expMeanmMZ", dimnames=selVars ) ), mxModel("DZm", #matrix for the observed age mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.AgeA", "data.AgeB"), name="obs_age" ), #matrix for the observed Education mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.EducationA", "data.EducationB"), name="obs_edu" ), #Algebra to add covariates to mean # mxAlgebra( expression= ACE.expMeanm + ACE.b_age %x% obs_age + ACE.b_edu %x% obs_edu , name="expMeanmDZ" ), mxAlgebra( expression= ACE.expMeanm , name="expMeanmDZ" ), mxData( observed=dzmData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovDZm", means="expMeanmDZ", dimnames=selVars ) ), mxModel("MZf", #matrix for the observed age mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.AgeA", "data.AgeB"), name="obs_age" ), #matrix for the observed Education mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.EducationA", "data.EducationB"), name="obs_edu" ), #Algebra to add covariates to mean # mxAlgebra( expression= ACE.expMeanf + ACE.b_age %x% obs_age + ACE.b_edu %x% obs_edu , name="expMeanfMZ" ), mxAlgebra( expression= ACE.expMeanf , name="expMeanfMZ" ), mxData( observed=mzfData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovMZf", means="expMeanfMZ", dimnames=selVars ) ), mxModel("DZf", #matrix for the observed age mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.AgeA", "data.AgeB"), name="obs_age" ), #matrix for the observed Education mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.EducationA", "data.EducationB"), name="obs_edu" ), #Algebra to add covariates to mean # mxAlgebra( expression= ACE.expMeanf + ACE.b_age %x% obs_age + ACE.b_edu %x% obs_edu , name="expMeanfDZ" ), mxAlgebra( expression= ACE.expMeanf , name="expMeanfDZ" ), mxData( observed=dzfData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovDZf", means="expMeanfDZ", dimnames=selVars ) ) ), mxAlgebra( expression=MZm.objective + DZm.objective + MZf.objective + DZf.objective , name="modelfit" ), mxAlgebraObjective("modelfit") #, mxCI (c('ACE.VarComp')) ) #ACENonScSLCoraFit <- mxRun(ACENonScSLCoraModel, intervals=F) ACENonScSLCoraFit <- mxRun(ACENonScSLCoraModel) ACENonScSLCoraSumm <- summary(ACENonScSLCoraFit) ACENonScSLCoraSumm