Hi!
I have a script from SGDP summer school, I can analyze using 2 continuous variables with this script. However, when I try using 3 variables (2 continuous 1 dummy coded), I get stuck while setting up the matrix. Before starting the analysis, I standardized my variables so that the SDs were 1 and their mean was 0.
Part of my code is as follows:
nv <- 3 # number of variables for a twin ntv <- 2*nv # number of variables for a pair Vars <- c('bmi', 'age','dummy_x') selVars <- c('bmi1','age1', 'dummy_x1', 'bmi2','age2', 'dummy_x2') # Select Data for Analysis mzData <- subset(TwinData, zygosity=='MZ', selVars) dzData <- subset(TwinData, zygosity=='DZ', selVars) # Descriptive Statistics of data by zygosity group describe(mzData) describe(dzData) colMeans(mzData,na.rm=TRUE) cov(mzData,use="complete") cor(mzData,use="complete") colMeans(dzData,na.rm=TRUE) cov(dzData,use="complete") cor(dzData,use="complete") # ----------------------------------------------------------------------------------------------------------------------------------- # Now we estimate the MZ and DZ covariances and Means by fitting a 'model' to the data # 1a) Specify Bivariate Saturated Model (Cholesky Decomposition) # ----------------------------------------------------------------------------------------------------------------------------------- # Matrix & Algebra for expected means and covariances MZlow <-mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=T, values=5, name="LowMZ" ) CovMZ <-mxAlgebra( expression=LowMZ %*% t(LowMZ), name="ExpCovMZ" ) MeanMZ <-mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=c(0,0,0,0,0,0), name="ExpMeanMZ" ) DZlow <-mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=T, values=5, name="LowDZ" ) CovDZ <-mxAlgebra( expression=LowDZ %*% t(LowDZ), name="ExpCovDZ" ) MeanDZ <-mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=c(0,0,0,0,0,0), name="ExpMeanDZ" ) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups objMZ <- mxExpectationNormal( covariance="ExpCovMZ", means="ExpMeanMZ", dimnames=selVars) objDZ <- mxExpectationNormal( covariance="ExpCovDZ", means="ExpMeanDZ", dimnames=selVars) fitFunction <- mxFitFunctionML() # Combine Groups modelMZ <- mxModel( MZlow, CovMZ, MeanMZ, dataMZ, objMZ, fitFunction, name="MZ") modelDZ <- mxModel( DZlow, CovDZ, MeanDZ, dataDZ, objDZ, fitFunction, name="DZ") minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) SatModel <- mxModel( "Sat", modelMZ, modelDZ, minus2ll, obj) # ------------------------------------------------------------------------------------------------------------------------------- # 1a) RUN Saturated Model (Cholesky Decomposition) SatFit <- mxRun(SatModel) (SatSum <- summary(SatFit)) # Generate SatModel Output round(SatFit@output$estimate,4) mxEval(MZ.ExpMeanMZ, SatFit) mxEval(MZ.ExpCovMZ, SatFit) mxEval(DZ.ExpMeanDZ, SatFit) mxEval(DZ.ExpCovDZ, SatFit) # ----------------------------------------------------------------------------------------------------------------------------------- # 1b) Specify Bivariate Saturated Model (Gaussian Decomposition) # We will use this specification to fit a constrained model # ----------------------------------------------------------------------------------------------------------------------------------- svM <-c(0,0,0,0,0,0) # these are start values svSD <-c(1,1,1,1,1,1) svMZ<-c(-.04,.2,-.4,-.1,.2,-.4) svDZ <-c(-.03,-.1,-.1,-.1,.5,-.3) # Matrix & Algebra for expected means and covariances MeanMZ <-mxMatrix( "Full", 1, ntv, free=T, values=svM, labels=c("MZm11", "MZm21","MZm31","MZm12", "MZm22",'Mzm32'), name="ExpMeanMZ" ) MZsd <-mxMatrix( "Diag", ntv, ntv, free=T, values=svSD,labels=c("MZs11", "MZs21","MZs31","MZs12", "MZs22",'Mzs32'), name="sdMZ" ) Cormz <-mxMatrix( "Stand", ntv, ntv, free=T, values=svMZ, labels=c("rPhMZ1","rMZ1","MZxtxt1","MZxtxt2","rMZ2","rPhMZ2"),name="MZCor") CovMZ <-mxAlgebra( expression=sdMZ %*% MZCor %*% t(sdMZ), name="ExpCovMZ" ) MeanDZ <-mxMatrix( "Full", 1, ntv, free=T, values=svM, labels=c("DZm11", "DZm21", "DZm12", "DZm22"), name="ExpMeanDZ" ) DZsd <-mxMatrix( "Diag", ntv, ntv, free=T, values=svSD, labels=c("DZs11", "DZs21", "DZs12", "DZs22"), name="sdDZ" ) Cordz <-mxMatrix( "Stand", ntv, ntv, free=T, values=svDZ, labels=c("rPhDZ1","rDZ1","DZxtxt1","DZxtxt2","rDZ2","rPhDZ2"),name="DZCor") CovDZ <-mxAlgebra( expression=sdDZ %*% DZCor %*% t(sdDZ), name="ExpCovDZ" )
First thing I don't understand is how do we determine the svMZ and svDZ values? How many values do we need to enter? Originally there were 6 for 2 variables, in this case do I have to enter 9 values for 3 variables?
Secondly, how do we write the elements of the matrix in the corMZ and corDZ part?
Thanks!