# # # #----------------------------------------------------------------------------------------------------------------- #R SCRIPTS USED FOR THE EMPIRICAL DATA APPLICATION IN #"Environment-by-PGS Interaction: Combining the classical twin design and Polygenic Scores # to test for Genotype x Environment Interaction" #BY BRUINS, HOTTENGA, NEALE, POOL, BOOMSMA, & DOLAN (2022) # # #OPENMX SCRIPTS FOR MODEL BY C.V. DOLAN #EDITED BY S.BRUINS # #----------------------------------------------------------------------------------------------------------------- # # # #----------------------------------------------------------------------------------------------------------------- #------------------------------------------------------PREPARATION------------------------------------------------ #----------------------------------------------------------------------------------------------------------------- require(OpenMx) mxOption(NULL, "Default optimizer","NPSOL") #DATA FORMAT: #colnames(mzdat) = varnamesmz = c('P1','T1','T2') #columns: PGS, Phenotype T1, Phenotype T2 #colnames(dzdat) = varnamesdz = c('P1','T1','P2','T2') #columns: PGS T1, Phenotype T1, PGS T2, Phenotype T2 #TV: for MZs---P1=AD PGS, T1=cog twin 1, T2=cog twin 2. for DZs---PGS T1= AD PGS twin 1, T1= cog twin 1, #TV: PGS T2= AD PGS twin 2, T2=cog twin 2. #data files with PGS and WL only library(haven) #data files with PGS, age, sleep, WL----------------------------------------------- #PGS= z_ad_pgs for twin 1 and twin 2 #Age= age_X for twin 1 and twin 2 #Sleep= slp_durX for twin 1 and twin 2 mzdat_wl_sleep <- read_sav("//data/BruinsMZfile_wordlist_sleep_age_pgs_17dec23.sav") View(mzdat_wl_sleep) newmzdat_wl_sleep<-as.data.frame(mzdat_wl_sleep) class(newmzdat_wl_sleep) dzdat_wl_sleep <- read_sav("//data/BruinsDZfile_wordlist_sleep_age_pgs_17dec23.sav") View(dzdat_wl_sleep) newdzdat_wl_sleep<-as.data.frame(dzdat_wl_sleep) class(newdzdat_wl_sleep) #new comment: reorganized the datafile to be in the exact same order as bruins' code # added age and sleep to these lines of code. Bruins had it in order of twin 1 vars # followed by twin 2 vars. followed the same format when adding in age and sleep colnames(newmzdat_wl_sleep)= varnamesmz = c('z_ad_pgs_1', 'age_X_1', 'slp_durX_1', 'T_wl_1', 'age_X_2','slp_durX_2', 'T_wl_2') #columns: PGS, age1, sleep1, phen t1, age2, sleep2, phen t2 colnames(newdzdat_wl_sleep)= varnamesdz = c('z_ad_pgs_1', 'age_X_1', 'slp_durX_1', 'T_wl_1', 'z_ad_pgs_2', 'age_X_2', 'slp_durX_2', 'T_wl_2') #columns: PGS1 age1, sleep1 phen t1, PGS2, age2, sleep2, phen t2 #attempt to bring in covariates #based on hermine's script (https://ibg.colorado.edu/cdrom2022/day3/practicalDay3.R) #need to create covVars for age and sleep. Bruins original script does not have covariates or additional moderators #Hermine's script only has one selVars even though she has data MZ and data DZ selVarsMZ <- c("z_ad_pgs_1", "T_wl_1", "T_wl_2") selVarsDZ <- c('z_ad_pgs_1', 'T_wl_1', 'z_ad_pgs_2', 'T_wl_2') covVars <- c('age_X_1', 'slp_durX_1', 'age_X_2', 'slp_durX_2') defLage <- mxMatrix( type="Full", nrow=1, ncol=3, free=FALSE, labels=c("data.ageX_1", "data.ageX_1", "data.ageX_2"), name="defage" ) #defL from hermine's pathbage <- mxMatrix( type="Full", nrow=1, ncol=3, free=c(FALSE,TRUE,TRUE), values=c(0, 0.01, 0.01) , label=c("bageA", "bageB", "bageB") , name="betaage" ) #A=pgs, B=Phen #created to match van der Sluis mplus approach with the treatment of sleep twin and co-twin #Clarification question: Why is it nrow=1 and ncol=3 for MZ but ncol=4 for DZ? Is this due to the # PGS appearing for only once for twin 1 and 2 for MZ group? yes. defLsleepMZ <- mxMatrix( type="Full", nrow=1, ncol=3, free=FALSE, labels=c("data.slp_durX_1", "data.slp_durX_1", "data.slp_durX_2"), name="defsleepMZ" ) #defL from hermine's pathbsleepMZ <- mxMatrix( type="Full", nrow=1, ncol=3, free=c(FALSE, TRUE, TRUE), values=c(0, 0.1, 0.01), label=c("bsleepAMZ", "bsleepBMZ", "bsleepBMZ"), name="betasleepMZ" ) defLsleep2MZ <- mxMatrix( type="Full", nrow=1, ncol=3, free=FALSE, labels=c("data.slp_durX_1", "data.slp_durX_2", "data.slp_durX_1"), name="defsleepMZ2" ) #defL from hermine's pathbsleep2MZ <- mxMatrix( type="Full", nrow=1, ncol=3, free=c(FALSE, TRUE, TRUE), values=c(0, 0.1, 0.01), label=c("bsleepAMZ", "bsleepBMZ2", "bsleepBMZ2"), name="betasleepMZ2" ) defLsleepDZ <- mxMatrix( type="Full", nrow=1, ncol=4, free=FALSE, labels=c("data.slp_durX_1", "data.slp_durX_1", "data.slp_durX_2", "data.slp_durX_2"), name="defsleepDZ" ) #defL from hermine's pathbsleepDZ <- mxMatrix( type="Full", nrow=1, ncol=4, free=c(FALSE, TRUE, FALSE, TRUE), values=c(0, 0.1, 0, 0.01), label=c("bsleepADZ", "bsleepBDZ","bsleepADZ", "bsleepBDZ"), name="betasleepDZ" ) defLsleep2DZ <- mxMatrix( type="Full", nrow=1, ncol=4, free=FALSE, labels=c("data.slp_durX_1", "data.slp_durX_2", "data.slp_durX_1", "data.slp_durX_2"), name="defsleepDZ2" ) #defL from hermine's pathbsleep2DZ <- mxMatrix( type="Full", nrow=1, ncol=4, free=c(FALSE, TRUE, FALSE, TRUE), values=c(0, 0.1, 0, 0.01), label=c("bsleepADZ", "bsleepBDZ2", "bsleepADZ", "bsleepBDZ2"), name="betasleepDZ2" ) #1x3 for MZ 1x4 for DZ bc of the PGS meanMZ <- mxMatrix( "Full", 1, 3, TRUE, values=c(0, 50, 50), labels=c("mPGS", "mMZ1","mMZ2"), name="meanMZ" ) meanDZ <- mxMatrix( "Full", 1, 4, TRUE, values=c(0, 50, 0, 50), labels=c("mPGS1", "mDZ1", "mPGS2", "mDZ2"), name="meanDZ" ) expMeanMZ <- mxAlgebra( expression= meanMZ + betaage%*%defage + betasleepMZ%*%defsleepMZ + betasleepMZ2%*%defsleepMZ2 , name="expMeanMZ" ) expMeanDZ <- mxAlgebra( expression= meanDZ + betaage%*%defage + betasleepDZ%*%defsleepDZ + betasleepDZ2%*%defsleepDZ2 , name="expMeanDZ" ) # we have to create a second dataset to bring the varnamesmz and covVars into MZ and DZ datasets dataMZ <- mxData( observed=newmzdat_wl_sleep[,c(selVarsMZ, covVars)], "raw" ) dataDZ <- mxData( observed=newdzdat_wl_sleep[,c(selVarsDZ, covVars)], "raw" ) #not sure if needed...also not fully correct. would need to modify if needed #modelMZ <-mxModel( defLage, pathbage, defLsleep, pathbsleep, meanMZ, expMeanMZ, covMZ, dataMZ, expMZ, funML, name="MZ" ) #modelDZ <- mxModel( defLage, pathbage, defLsleep, pathbsleep, meanDZ, expMeanDZ, covDZ, dataDZ, expDZ, funML, name="DZ" ) # #Note from Bruin's code: the version of this script as provided on OSF is set to fit an AE model with E-by-PGS-interaction. #To include C in the model, change starting values and free parameters accordingly. # #starting values: c(From latent PGS to observed PGS, from latent PGS to observed phenotype, from latent A/C/E variance to observed phenotype) #constants a_st<-c(1,.01,.7) c_st<-c(0,0,.0) e_st<-c(0,0,.7) #moderation betas b1a_st<-c(0,0,0)#pgs b1c_st<-c(0,0,.0) b1e_st<-c(0,0,.01) #add b's for other covariates b2a_st<-c(0,0,0) #age b2c_st<-c(0,0,.0) b2e_st<-c(0,0,.01) b3a_st<-c(0,0,0) #sleep b3c_st<-c(0,0,.0) b3e_st<-c(0,0,.01) # nphen1<-2 nphen2<-nphen1*2 # # #----------------------------------------------------------------------------------------------------------------- #-----------------------------------------------------OPENMX MODEL------------------------------------------------ #----------------------------------------------------------------------------------------------------------------- # ACEsumstatMZ <- mxModel("ACEMZ", # # Matrix expMean for expected mean vector for MZ and DZ twins mxMatrix(type="Full", nrow=1, ncol=3, free=c(TRUE,TRUE,TRUE),values=c(0,0,0),label=c("mprs","mPh","mPh"), name="expMeanMZ"), #bruin's script has "expMean" for both MZ & DZ models. # # Matrices a, c, and e to store the a, c, and e path coefficients mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,TRUE,TRUE), values=a_st, label=c("a11","a21","a22"),name="a"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=c_st, #FALSE,FALSE,TRUE label=c("c11","c21","c22"),name="c"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,TRUE), values=e_st, label=c("e11","e21","e22"),name="e"), mxMatrix(type="Full", nrow=3, ncol=4, free=FALSE, values= matrix(c(1,0,0,0, 0,1,0,0, 0,0,0,1),3,4,byrow=TRUE),name="F"), #filter matrix only for MZ for one PGS #moderation coefficients mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=b1a_st, label=c("b1a11","b1a21","b1a22"),name="b1a"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=b1c_st, #FALSE,FALSE,TRUE label=c("b1c11","b1c21","b1c22"),name="b1c"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,TRUE), values=b1e_st, label=c("b1e11","b1e21","b1e22"),name="b1e"), #AGE mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,TRUE), values=b2a_st, label=c("b2a11","b2a21","b2a22"),name="b2a"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=b2c_st, #FALSE,FALSE,TRUE label=c("b2c11","b2c21","b2c22"),name="b2c"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,TRUE), values=b2e_st, label=c("b2e11","b2e21","b2e22"),name="b2e"), #SLEEP mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,TRUE), values=b3a_st, label=c("b3a11","b3a21","b3a22"),name="b3a"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=b3c_st, #FALSE,FALSE,TRUE label=c("b3c11","b3c21","b3c22"),name="b3c"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,TRUE), values=b3e_st, label=c("b3e11","b3e21","b3e22"),name="b3e"), #moderator mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.z_ad_pgs_1","data.z_ad_pgs_1"),name="Modp1"), mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.z_ad_pgs_1","data.z_ad_pgs_1"),name="Modp2"), #P1 for both twins because this is MZ mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.age_X_1","data.age_X_1"),name="Moda1"), mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.age_X_2","data.age_X_2"),name="Moda2"), mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.slp_durX_1","data.slp_durX_1"),name="Mods1"), mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.slp_durX_2","data.slp_durX_2"),name="Mods2"), # # Matrices A, C, and E to compute A, C, and E variance components #leaving correction moderation b for PGS -stays the same # mxAlgebra( expression=( b1e[2,2]*( (a[2,1]^2+a[2,2]^2)/a[2,1]^2 )^.5 ), name="w1be"), #corrected moderation beta # mxAlgebra( expression=( b1c[2,2]*( (a[2,1]^2+a[2,2]^2)/a[2,1]^2 )^.5 ), name="w1bc"), #corrected moderation beta #basically moderation effect of PGS for the a component, even if it isn't being #estimated, will still do for age and sleep. #exactly the same on both sides (the transposed side) because within twin #twin 1 mxAlgebra( expression=(a+b1a%*%Modp1 + b2a%*%Moda1 + b3a%*%Mods1) %*% t((a+b1a%*%Modp1 + b2a%*%Moda1 + b3a%*%Mods1)), name="A11"), # a^2 #exactly the same on both sides because within twin mxAlgebra( expression=(c+b1c%*%Modp1 + b2c%*%Moda1 + b3c%*%Mods1) %*% t((c+b1c%*%Modp1 + b2c%*%Moda1 + b3c%*%Mods1)), name="C11"), # c^2 mxAlgebra( expression=(e+b1e%*%Modp1 + b2e%*%Moda1 + b3e%*%Mods1) %*% t((e+b1e%*%Modp1 + b2e%*%Moda1 + b3e%*%Mods1)), name="E11"), # e^2 #twin 2. mxAlgebra( expression=(a+b1a%*%Modp2 + b2a%*%Moda2 + b3a%*%Mods2) %*% t((a+b1a%*%Modp2 +b2a%*%Moda2 + b3a%*%Mods2)), name="A22"), # a^2 mxAlgebra( expression=(c+b1c%*%Modp2 + b2c%*%Moda2 + b3c%*%Mods2) %*% t((c+b1c%*%Modp2 +b2c%*%Moda2 + b3c%*%Mods2)), name="C22"), # c^2 mxAlgebra( expression=(e+b1e%*%Modp2 + b2e%*%Moda2 + b3e%*%Mods2) %*% t((e+b1e%*%Modp2 +b2e%*%Moda2 + b3e%*%Mods2)), name="E22"), # e^2 # twin 1 on the left, twin 2 on the right mxAlgebra( expression=(a+b1a%*%Modp1 + b2a%*%Moda1 +b3a%*%Mods1) %*% t((a+b1a%*%Modp2 +b2a%*%Moda2 +b3a%*%Mods2 )), name="A12"), # a^2 mxAlgebra( expression=(c+b1c%*%Modp1 + b2c%*%Moda1 +b3c%*%Mods1) %*% t((c+b1c%*%Modp2 +b2c%*%Moda2 +b3c%*%Mods2 )), name="C12"), # c^2 mxAlgebra( expression=(e+b1e%*%Modp1 + b2e%*%Moda1 +b3e%*%Mods1) %*% t((e+b1e%*%Modp2 +b2e%*%Moda2 +b3e%*%Mods2 )), name="E12"), # e^2 # Matrix expCovMZ for expected covariance matrix for MZ twins # mxAlgebra( expression= F%*% rbind( cbind(A11+C11+E11, A12+C12), cbind(t(A12+C12), A22+C22+E22))%*%t(F),name="expCovMZ"), mxData( observed=dataMZ, type="raw"), # the data, in bruin's adapted it was newmzdat_wl mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", #in bruins it was expMean dimnames=selVarsMZ), # the fit function #in bruins it was varnamesMZ mxFitFunctionML() ) ACEsumstatDZ <- mxModel("ACEDZ", # # Matrix expMean for expected mean vector for MZ and DZ twins mxMatrix(type="Full", nrow=1, ncol=4, free=c(TRUE,TRUE,TRUE,TRUE),values=c(0,0,0,0),label=c("mprs","mPh",'mprs',"mPh"), name="expMeanDZ"), #bruins code has expMean should we change to expMeanDZ?, # Matrices a, c, and e to store the a, c, and e path coefficients #constants mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,TRUE,TRUE), values=a_st, label=c("a11","a21","a22"),name="a"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=c_st, #FALSE,FALSE,TRUE label=c("c11","c21","c22"),name="c"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,TRUE), values=e_st, label=c("e11","e21","e22"),name="e"), # #moderation betas mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=b1a_st, label=c("b1a11","b1a21","b1a22"),name="b1a"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=b1c_st, label=c("b1c11","b1c21","b1c22"),name="b1c"), #FALSE,FALSE,TRUE mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,TRUE), values=b1e_st, label=c("b1e11","b1e21","b1e22"),name="b1e"), #AGE mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=b2a_st, label=c("b2a11","b2a21","b2a22"),name="b2a"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=b2c_st, label=c("b2c11","b2c21","b2c22"),name="b2c"), #FALSE,FALSE,TRUE mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,TRUE), values=b2e_st, label=c("b2e11","b2e21","b2e22"),name="b2e"), #SLEEP mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=b3a_st, label=c("b3a11","b3a21","b3a22"),name="b3a"), mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,FALSE), values=b3c_st, label=c("b3c11","b3c21","b3c22"),name="b3c"), #FALSE,FALSE,TRUE mxMatrix(type="Lower", nrow=nphen1, ncol=nphen1, free=c(FALSE,FALSE,TRUE), values=b3e_st, label=c("b3e11","b3e21","b3e22"),name="b3e"), # #moderator mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.z_ad_pgs_1","data.z_ad_pgs_1"),name="Modp1"), mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.z_ad_pgs_2","data.z_ad_pgs_2"),name="Modp2"), mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.age_X_1","data.age_X_1"),name="Moda1"), mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.age_X_2","data.age_X_1"),name="Moda2"), mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.slp_durX_1","data.slp_durX_1"),name="Mods1"), mxMatrix(type="Diag", nrow=nphen1, ncol=nphen1, free=FALSE, label=c("data.slp_durX_2","data.slp_durX_2"),name="Mods2"), # # Matrices A, C, and E to compute A, C, and E variance components # leaving correction moderation b for PGS- stays the same #mxAlgebra( expression=( b1e[2,2]*( (a[2,1]^2+a[2,2]^2)/a[2,1]^2 )^.5 ), name="w1be"), #corrected moderation beta #mxAlgebra( expression=( b1c[2,2]*( (a[2,1]^2+a[2,2]^2)/a[2,1]^2 )^.5 ), name="w1bc"), #corrected moderation beta #expanding to include age and sleep #Twin 1 mxAlgebra( expression=(a+b1a%*%Modp1 + b2a%*%Moda1 + b3a%*%Mods1) %*% t((a+b1a%*%Modp1 + b2a%*%Moda1 + b3a%*%Mods1)), name="A11"), # a^2 mxAlgebra( expression=(c+b1c%*%Modp1 + b2c%*%Moda1 + b3c%*%Mods1) %*% t((c+b1c%*%Modp1 + b2c%*%Moda1 + b3c%*%Mods1)), name="C11"), # c^2 mxAlgebra( expression=(e+b1e%*%Modp1 + b2e%*%Moda1 + b3e%*%Mods1) %*% t((e+b1e%*%Modp1 + b2e%*%Moda1 + b3e%*%Mods1)), name="E11"), # e^2 # TWin 2 mxAlgebra( expression=(a+b1a%*%Modp2 +b2a%*%Moda2 + b3a%*%Mods2) %*% t((a+b1a%*%Modp2 +b2a%*%Moda2 + b3a%*%Mods2)), name="A22"), # a^2 mxAlgebra( expression=(c+b1c%*%Modp2 +b2c%*%Moda2 + b3c%*%Mods2) %*% t((c+b1c%*%Modp2 +b2c%*%Moda2 + b3c%*%Mods2)), name="C22"), # c^2 mxAlgebra( expression=(e+b1e%*%Modp2 +b2e%*%Moda2 + b3e%*%Mods2) %*% t((e+b1e%*%Modp2 +b2e%*%Moda2 + b3e%*%Mods2)), name="E22"), # e^2 # Twin 1 on the left, twin 2 on the right mxAlgebra( expression=(a+b1a%*%Modp1 +b2a%*%Moda1 +b3a%*%Mods1) %*% t((a+b1a%*%Modp2 +b2a%*%Moda2 +b3a%*%Mods2)), name="A12"), # a^2 mxAlgebra( expression=(c+b1c%*%Modp1 +b2c%*%Moda1 +b3c%*%Mods1) %*% t((c+b1c%*%Modp2 +b2c%*%Moda2 +b3c%*%Mods2)), name="C12"), # c^2 mxAlgebra( expression=(e+b1e%*%Modp1 +b2e%*%Moda1 +b3e%*%Mods1) %*% t((e+b1e%*%Modp2 +b2e%*%Moda2 +b3e%*%Mods2)), name="E12"), # e^2 # Matrix expCovDZ for expected covariance matrix for DZ twins # mxAlgebra( expression= rbind( cbind(A11+C11+E11, .5%x%A12+C12), cbind(t(.5%x%A12+C12), A22+C22+E22)),name="expCovDZ"), mxData( observed=dataDZ, type="raw"), mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVarsDZ), mxFitFunctionML() ) # twinACEModel2 <- mxModel("twinACE", ACEsumstatMZ, ACEsumstatDZ, mxFitFunctionMultigroup( c("ACEMZ","ACEDZ")) ) # # fit the model twinACEFit2 <- mxRun(twinACEModel2) summary(twinACEFit2) #----------------------------------------------------------------------------------------------------------------- #-------------------------------------------------SIGNIFICANCE TESTS---------------------------------------------- #----------------------------------------------------------------------------------------------------------------- # #test e moderation twinACEModel2fixe <- omxSetParameters(twinACEFit2, labels=c('be22'), free=FALSE, values=0) twinACEFit2fixe <- mxTryHard(twinACEModel2fixe,10) (testbe<-mxCompare(twinACEFit2,twinACEFit2fixe)) summary(twinACEFit2fixe) #test c moderation #twinACEModel2fixc <- omxSetParameters(twinACEFit2, labels=c('bc22'), free=FALSE, values=0) #twinACEFit2fixc <- mxTryHard(twinACEModel2fixc,10) #(testbc=mxCompare(twinACEFit2,twinACEFit2fixc)) #summary(twinACEFit2fixc) #