phenotypic correlation

11 posts / 0 new

Last post

Mon, 06/18/2018 - 09:37

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

phenotypic correlation

Hello,

I am a little bit confused about how to calculate the phenotypic correlation between variable 1 and variable 2, in a bivariate model with one binary variable and one continuous variable.

I mean if I use this: mxEval(cov2cor(V), fitADE, T). I get a value of 0.200 whereas if I calculate the correlation in SPSS (also in R) I get a value of 0.150.

I do not know if I am doing something wrong. Should the results be the same right?
I have tried with different methods, with and without covariables etc… But I get the same value.

Can I get this value from the saturated model?

Thank you so much in advance.

Mon, 06/18/2018 - 11:21

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

biserial vs point-biserial

What you're calculating from OpenMx is a model-expected biserial correlation, which is the correlation between the continuous variable and the normally distributed continuum underlying the binary variable. In contrast, directly calculating the correlation using the default method in most statistical software will give you a point-biserial correlation, which is simply the Pearson correlation between a continuous variable and a binary variable numerically coded with zeroes and ones. In general, the two values won't be the same.

Mon, 06/18/2018 - 13:28

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Thank you so much.

I see, it is the difference between biserial correlation and point-biserial correlation.

Also, Open Mx uses tetrachoric for the correlation between two binary variables and polychoric for ordinal variables right?

What result should I report? I mean what is more correct I know the difference is not too big but I would like to know what is the correct correlation.

Mon, 06/18/2018 - 15:07

(Reply to #3) #4

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

biserial

I see, it is the difference between biserial correlation and point-biserial correlation.

Also, Open Mx uses tetrachoric for the correlation between two binary variables and polychoric for ordinal variables right?

Essentially, yes. In completely saturated model and with no missing data, the values you get from OpenMx for all four of those types of correlation should match what you get when calculating those correlations a different way. But in general, that won't be true if you're fitting a non-saturated model with OpenMx.

What result should I report? I mean what is more correct I know the difference is not too big but I would like to know what is the correct correlation.

I think reporting the biserial correlation makes more sense, since it reflects what your OpenMx analyses assume about a normally distributed continuum underlying the binary trait. If you want to report it merely as a descriptive statistic, then report it from the saturated model. If it's somehow relevant to your substantive conclusions, then you might want to report it from the "best" model.

Mon, 06/18/2018 - 15:46

(Reply to #4) #5

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Thank you!

And how can I get the result from the saturated model?

I mean I can get the correlations with:

mxEval(cov2cor(DZ.covDZ), fit, T) and mxEval(cov2cor(MZ.covMZ), fit, T)

But I do not know how to get the phenotypic correlation for the whole sample at the same time in the saturated model.

Thank you so much in advance.

Mon, 06/25/2018 - 11:39

(Reply to #5) #6

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

highest-dimensional I guess

Oh, I see. Then I guess I would report the correlation as calculated from the highest-dimensional (most parameters) model that constrains MZ and DZ twin means and variances.

Mon, 06/25/2018 - 12:15

(Reply to #6) #7

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Perfect!!

Thank you so much Rob!

Tue, 08/14/2018 - 14:26

(Reply to #7) #8

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Hello again,

I am a little bit confused since I am trying to get the phenotypic correlation from the saturated model but I am not sure if I am doing it correctly.

I am trying to get the phenotypic correlation from the highest-dimensional model as you told me but unless I equate all the means, variances and covariances. I get different phenotypic correlation for MZ and DZ. I mean when I run these lines:
mxEval(cov2cor(DZ.covDZ), fitEMVZ, T)
mxEval(cov2cor(MZ.covMZ), fitEMVZ, T)

I get different values (for MZ and DZ) for the phenotypic correlations. The values are similar but not the same

Here my script:

# Select Continuous Variables
vars      <- c('DP_Ln',"AN_Ln", 'EX_Ln')                   # list of continuous variables names
nvc       <- 3                         # number of continuous variables
ntvc      <- nvc*2                     # number of total continuous variables
conVars   <- paste(vars,c(rep(1,nvc),rep(2,nvc)),sep="")
 
# Select Ordinal Variables
nth       <- 1                         # number of thresholds
vars      <- c('APN')                   # list of ordinal variables names
nvo       <- 1                         # number of ordinal variables
ntvo      <- nvo*2                     # number of total ordinal variables
ordVars   <- paste(vars,c(rep(1,nvo),rep(2,nvo)),sep="")
ordData   <- twinData
 
# Select Variables for Analysis
vars      <- c('DP_Ln',"AN_Ln", 'EX_Ln','APN')             # list of variables names
nv        <- nvc+nvo                   # number of variables
ntv       <- nv*2                      # number of total variables
selVars   <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
 
# Select Covariates for Analysis
 
covVars   <- c('age1',"age2", "Sex1" , "Sex2")
nc        <- 4                         # number of covariates   
 
# Select Data for Analysis
mzData    <- subset(ordData, Zyg==1|Zyg==3, c(selVars, covVars))
dzData    <- subset(ordData, Zyg==2|Zyg==4| Zyg==5 , c(selVars, covVars))
#mzDataF   <- cbind(mzData[,conVars],mxFactor( x=mzData[,ordVars], levels=c(0:nth)), age=mzData[,covVars] ) 
#dzDataF   <- cbind(dzData[,conVars],mxFactor( x=dzData[,ordVars], levels=c(0:nth)), age=dzData[,covVars] ) 
mzDataF   <- mzData 
dzDataF   <- dzData
mzDataF$APN1 <- mxFactor(mzDataF$APN1, levels =c(0,1))
mzDataF$APN2 <- mxFactor(mzDataF$APN2, levels =c(0,1))
dzDataF$APN1 <- mxFactor(dzDataF$APN1, levels =c(0,1))
dzDataF$APN2 <- mxFactor(dzDataF$APN2, levels =c(0,1))
 
 
# Set Starting Values
frMV      <- c(TRUE,TRUE, TRUE, FALSE)             # free status for variables
frCvD     <- diag(frMV,ntv,ntv)        # lower bounds for diagonal of covariance matrix
frCvD[lower.tri(frCvD)] <- TRUE        # lower bounds for below diagonal elements
frCvD[upper.tri(frCvD)] <- TRUE        # lower bounds for above diagonal elements
frCv      <- matrix(as.logical(frCvD),4)
svMe      <- c(11,18,11,0.2)                   # start value for means
svVa      <- c(2,1,1,1)                    # start value for variance
lbVa      <- .0001                     # start value for lower bounds
svVas     <- diag(svVa,ntv,ntv)        # assign start values to diagonal of matrix
lbVas     <- diag(lbVa,ntv,ntv)        # assign lower bounds to diagonal of matrix
svLTh     <- -2                      # start value for first threshold
svITh     <- 0.1                         # start value for increments
svTh      <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv)     # start value for thresholds
lbTh      <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv)     # lower bounds for thresholds
 
# Create Labels
labMeMZ   <- paste("meanMZ",selVars,sep="_")
labMeDZ   <- paste("meanDZ",selVars,sep="_")
labMeZ    <- paste("meanZ",selVars,sep="_")
labCvMZ   <- labLower("covMZ",ntv)
labCvDZ   <- labLower("covDZ",ntv)
labCvZ    <- labLower("covZ",ntv)
labVaMZ   <- labDiag("covMZ",ntv)
labVaDZ   <- labDiag("covDZ",ntv)
labVaZ    <- labDiag("covZ",ntv)
 
labThMZ   <- paste(paste("t",1:nth,"MZ",sep=""),rep(ordVars,each=nth),sep="_")
labThDZ   <- paste(paste("t",1:nth,"DZ",sep=""),rep(ordVars,each=nth),sep="_")
labThZ    <- paste(paste("t",1:nth,"Z",sep=""),rep(ordVars,each=nth),sep="_")
labBe     <- labFull("beta",nc,1)
 
# ------------------------------------------------------------------------------
# PREPARE MODEL
 
# Saturated Model
# Create Matrices for Covariates and linear Regression Coefficients
defAge    <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age1", "data.age1", "data.age1","data.age1","data.age2", "data.age2", "data.age2","data.age2"), name="Age" )
pathB1     <- mxMatrix( type="Full", nrow=4, ncol=1, free=TRUE, values=.01, label=c("b11", "b12","b13","b14", "b11", "b12", "b13","b14"), name="b1" )
 
defSex    <- mxMatrix( type="Full", nrow=1, ncol=8, free=FALSE, labels=c("data.Sex1", "data.Sex1","data.Sex1", "data.Sex1","data.Sex2", "data.Sex2", "data.Sex2","data.Sex2"), name="Sex" )
pathB2     <- mxMatrix( type="Full", nrow=1, ncol=8, free=TRUE, values=.01, label=c("b21", "b22","b23","b24", "b21", "b22", "b23","b24"), name="b2" )
 
# Create Algebra for expected Mean & Threshold Matrices
meanMZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeMZ, name="meanMZ" )
meanDZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeDZ, name="meanDZ" )
expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanMZ" )
expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanDZ" )
thinMZ    <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=0, name="thinMZ" )
thinDZ    <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=0, name="thinDZ" )
inc       <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" )
threMZ    <- mxAlgebra( expression= inc %*% thinMZ, name="threMZ" )
threDZ    <- mxAlgebra( expression= inc %*% thinDZ, name="threDZ" )
 
# Create Algebra for expected Covariance Matrices
covMZ     <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" )
covDZ     <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" )
 
# Create Data Objects for Multiple Groups
dataMZ    <- mxData( observed=mzDataF, type="raw" )
dataDZ    <- mxData( observed=dzDataF, type="raw" )
 
# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="covMZ", means="expMeanMZ", dimnames=selVars, thresholds="threMZ", threshnames=ordVars )
expDZ     <- mxExpectationNormal( covariance="covDZ", means="expMeanDZ", dimnames=selVars, thresholds="threDZ", threshnames=ordVars )
funML     <- mxFitFunctionML()
 
# Create Model Objects for Multiple Groups
pars      <- list(pathB1,pathB2, inc )
defs      <- list(defAge,defSex)
modelMZ   <- mxModel( name="MZ", pars, defs, meanMZ, expMeanMZ, covMZ, thinMZ, threMZ, dataMZ, expMZ, funML )
modelDZ   <- mxModel( name="DZ", pars, defs, meanDZ, expMeanDZ, covDZ, thinDZ, threDZ, dataDZ, expDZ, funML )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
 
# Create Confidence Interval Objects
ciCor     <- mxCI( c('MZ.covMZ','DZ.covDZ' ))
ciThre    <- mxCI( c('MZ.threMZ','DZ.threDZ' ))
MZZZ       <-mxAlgebra(expression =cov2cor(MZ.covMZ), name="MZZ")
DZZZ       <-mxAlgebra(expression =cov2cor(DZ.covDZ), name="DZZ")
ci2     <- mxCI( c('MZZ','DZZ') )
 
# Build Saturated Model with Confidence Intervals
model     <- mxModel( "twoSATja", pars, modelMZ, modelDZ, multi, MZZZ, DZZZ, ci2 )
 
# ------------------------------------------------------------------------------
# RUN MODEL
 
# Run Saturated Model
fit       <- mxTryHardOrdinal( model, intervals=T, extraTries=31 )
sum       <- summary( fit )
round(fit$output$estimate,4)
fitGofs(fit)
 
# ------------------------------------------------------------------------------
# RUN SUBMODELS
 
 
# Constrain expected Thresholds to be equal across twin order
modelEMTVO  <- mxModel( fit, name="eqThresTwin" )
 
 
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("DZ.thinDZ[1,1]","DZ.thinDZ[1,2]"), newlabels=c("THINDZ"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("MZ.thinMZ[1,1]","MZ.thinMZ[1,2]"), newlabels=c("THINMZ"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("meanMZ_DP1","meanMZ_DP2"), newlabels=c("meanMZ_DP"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("meanDZ_DP1","meanDZ_DP2"), newlabels=c("meanDZ_DP"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("meanMZ_EX1","meanMZ_EX2"), newlabels=c("meanMZ_EX"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("meanDZ_EX1","meanDZ_EX2"), newlabels=c("meanDZ_EX"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("meanMZ_AN1","meanMZ_AN2"), newlabels=c("meanMZ_AN"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("meanDZ_AN1","meanDZ_AN2"), newlabels=c("meanDZ_AN"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covMZ_1_1","covMZ_5_5"), newlabels=c("covMZDP"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covDZ_1_1","covDZ_5_5"), newlabels=c("covDZDP"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covMZ_2_1","covMZ_6_5"), newlabels=c("covMZDPAN"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covDZ_2_1","covDZ_6_5"), newlabels=c("covDZDPAN"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covMZ_2_2","covMZ_6_6"), newlabels=c("covMZAN"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covDZ_2_2","covDZ_6_6"), newlabels=c("covDZAN"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covMZ_3_1","covMZ_7_5"), newlabels=c("covMZDPEX"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covDZ_3_1","covDZ_7_5"), newlabels=c("covDZDPEX"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covMZ_3_2","covMZ_7_6"), newlabels=c("covMZANEX"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covDZ_3_2","covDZ_7_6"), newlabels=c("covDZANEX"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covMZ_3_3","covMZ_7_7"), newlabels=c("covMZEX"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covDZ_3_3","covDZ_7_7"), newlabels=c("covDZEX"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covMZ_4_1","covMZ_8_5"), newlabels=c("covMZDPAPN"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covDZ_4_1","covDZ_8_5"), newlabels=c("covDZDPAPN"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covMZ_4_2","covMZ_8_6"), newlabels=c("covMZANAPN"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covDZ_4_2","covDZ_8_6"), newlabels=c("covDZANAPN"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covMZ_4_3","covMZ_8_7"), newlabels=c("covMZEXAPN"))
modelEMTVO <- omxSetParameters( modelEMTVO, label=c("covDZ_4_3","covDZ_8_7"), newlabels=c("covDZEXAPN"))
modelEMTVO  <- omxAssignFirstParameters( modelEMTVO )
fitEMTVO   <- mxRun( modelEMTVO, intervals=F )
mxCompare(fit, subs <- list( modelEMTVO) )
round(fitEMVO$output$estimate,4)
fitGofs(fitEMTVO)
mxCompare (fit, fitEMTVO) 
 
# Constrain expected Means and Variances to be equal across twin order and zygosity
modelEMVZ <- mxModel(fitEMTVO, name="eqMVarsZyg" )
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("meanMZ_DP","meanDZ_DP"), newlabels=c("MEANDP"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("meanMZ_AN","meanDZ_AN"), newlabels=c("MEANAN"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("meanMZ_EX","meanDZ_EX"), newlabels=c("MEANEX"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZDP","covDZDP"), newlabels=c("COVDP"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZDEAN","covDZDE
AN"), newlabels=c("COVDPAN"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZAN","covDZAN"), newlabels=c("COVAN"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZDPEX","covDZDPEX"), newlabels=c("COVDPEX"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZANEX","covDZANEX"), newlabels=c("COVANEX"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZEX","covDZEX"), newlabels=c("COVEX"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZDPAPN","covDZDPAPN"), newlabels=c("COVDPAPN"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZANAPN","covDZANAPN"), newlabels=c("COVANAPN"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZEXAPN","covDZEXAPN"), newlabels=c("COVEXAPN"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_5_1","covDZ_5_1"), newlabels=c("COVDP1DP2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_5_2","covDZ_5_2"), newlabels=c("COVAN1DP2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_5_3","covDZ_5_3"), newlabels=c("COVEX1DP2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_5_4","covDZ_5_4"), newlabels=c("COVAPN1DP2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_6_1","covDZ_6_1"), newlabels=c("COVDP1AN2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_6_2","covDZ_6_2"), newlabels=c("COVAN1AN2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_6_3","covDZ_6_3"), newlabels=c("COVEX1AN2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_6_4","covDZ_6_4"), newlabels=c("COVAPNAN2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_7_1","covDZ_7_1"), newlabels=c("COVDP1EX2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_7_2","covDZ_7_2"), newlabels=c("COVAN1EX2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_7_3","covDZ_7_3"), newlabels=c("COVEX1EX2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_7_4","covDZ_7_4"), newlabels=c("COVAPN1EX2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_8_1","covDZ_8_1"), newlabels=c("COVDP1APN2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_8_2","covDZ_8_2"), newlabels=c("COVAN1APN2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_8_3","covDZ_8_3"), newlabels=c("COVEX1APN2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("covMZ_8_4","covDZ_8_4"), newlabels=c("COVAPN1APN2"))
modelEMVZ <- omxSetParameters( modelEMVZ, label=c("THINMZ","THINDZ"), newlabels=c("THIN"))
modelEMVZ  <- omxAssignFirstParameters( modelEMVZ )
fitEMVZ   <- mxRun( modelEMVZ, intervals=F )
 
 
mxCompare (fit, fitEMVZ) 
 
 
mxEval(cov2cor(DZ.covDZ), fitEMVZ, T) 
mxEval(cov2cor(MZ.covMZ), fitEMVZ, T)

But I do not know if equate all the covariances is correct. What parameters should I equate in models EMTVO and EMVZ to do it correctly? And to get the phenotypic correlation for the whole sample?

By the way, should in the univariate models the correlations for MZ (also DZ) should be the same across the models ACE AE and saturated if there is no missing data? I mean in the univariate model I usually get a value for rMZ (for example 0.70) and in the ACE a value slightly different (0.65) I do not know why the value is different form one model to another.
Thank you so much and sorry for the inconveniences

Tue, 08/14/2018 - 15:37

(Reply to #8) #9

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Mistake

I have seen that probably my script was wrong since the ordinal variable is dichotomous and I have fixed its variance.

I have changed the script (I am not sure if correctly). But I still have different values for MZ y DZ.
Here the script modified:

# Select Continuous Variables
vars      <- c('DP_Ln',"AN_Ln", 'EX_Ln')                   # list of continuous variables names
nvc       <- 3                         # number of continuous variables
ntvc      <- nvc*2                     # number of total continuous variables
conVars   <- paste(vars,c(rep(1,nvc),rep(2,nvc)),sep="")
 
# Select Ordinal Variables
nth       <- 1                         # number of thresholds
vars      <- c('APN')                   # list of ordinal variables names
nvo       <- 1                         # number of ordinal variables
ntvo      <- nvo*2                     # number of total ordinal variables
ordVars   <- paste(vars,c(rep(1,nvo),rep(2,nvo)),sep="")
ordData   <- twinData
 
# Select Variables for Analysis
vars      <- c('DP_Ln',"AN_Ln", 'EX_Ln','APN')             # list of variables names
nv        <- nvc+nvo                   # number of variables
ntv       <- nv*2                      # number of total variables
selVars   <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
 
# Select Covariates for Analysis
#ordData[,'age'] <- ordData[,'age']/100
#ordData   <- ordData[-which(is.na(ordData$age)),]
covVars   <- c('age1',"age2", "Sex1" , "Sex2")
nc        <- 4                         # number of covariates   
 
# Select Data for Analysis
mzData    <- subset(ordData, Zyg==1|Zyg==3, c(selVars, covVars))
dzData    <- subset(ordData, Zyg==2|Zyg==4| Zyg==5 , c(selVars, covVars))
#mzDataF   <- cbind(mzData[,conVars],mxFactor( x=mzData[,ordVars], levels=c(0:nth)), age=mzData[,covVars] ) 
#dzDataF   <- cbind(dzData[,conVars],mxFactor( x=dzData[,ordVars], levels=c(0:nth)), age=dzData[,covVars] ) 
mzDataF   <- mzData 
dzDataF   <- dzData
mzDataF$APN1 <- mxFactor(mzDataF$APN1, levels =c(0,1))
mzDataF$APN2 <- mxFactor(mzDataF$APN2, levels =c(0,1))
dzDataF$APN1 <- mxFactor(dzDataF$APN1, levels =c(0,1))
dzDataF$APN2 <- mxFactor(dzDataF$APN2, levels =c(0,1))
 
 
# Set Starting Values
frMV      <- c(TRUE,TRUE, TRUE, FALSE)             # free status for variables
frCvD     <- diag(frMV,ntv,ntv)        # lower bounds for diagonal of covariance matrix
frCvD[lower.tri(frCvD)] <- TRUE        # lower bounds for below diagonal elements
frCvD[upper.tri(frCvD)] <- TRUE        # lower bounds for above diagonal elements
frCv      <- matrix(as.logical(frCvD),4)
svMe      <- c(11,18,11,0.2)                   # start value for means
svVa      <- c(2,1,1,1)                    # start value for variance
lbVa      <- .0001                     # start value for lower bounds
svVas     <- diag(svVa,ntv,ntv)        # assign start values to diagonal of matrix
lbVas     <- diag(lbVa,ntv,ntv)        # assign lower bounds to diagonal of matrix
svLTh     <- -2                      # start value for first threshold
svITh     <- 0.1                         # start value for increments
svTh      <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv)     # start value for thresholds
lbTh      <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv)     # lower bounds for thresholds
 
# Create Labels
labMeMZ   <- paste("meanMZ",selVars,sep="_")
labMeDZ   <- paste("meanDZ",selVars,sep="_")
labMeZ    <- paste("meanZ",selVars,sep="_")
labCvMZ   <- labLower("covMZ",ntv)
labCvDZ   <- labLower("covDZ",ntv)
labCvZ    <- labLower("covZ",ntv)
labVaMZ   <- labDiag("covMZ",ntv)
labVaDZ   <- labDiag("covDZ",ntv)
labVaZ    <- labDiag("covZ",ntv)
 
labThMZ   <- paste(paste("t",1:nth,"MZ",sep=""),rep(ordVars,each=nth),sep="_")
labThDZ   <- paste(paste("t",1:nth,"DZ",sep=""),rep(ordVars,each=nth),sep="_")
labThZ    <- paste(paste("t",1:nth,"Z",sep=""),rep(ordVars,each=nth),sep="_")
labBe     <- labFull("beta",nc,1)
 
# ------------------------------------------------------------------------------
# PREPARE MODEL
 
# Saturated Model
# Create Matrices for Covariates and linear Regression Coefficients
defAge    <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age1", "data.age1", "data.age1","data.age1","data.age2", "data.age2", "data.age2","data.age2"), name="Age" )
pathB1     <- mxMatrix( type="Full", nrow=4, ncol=1, free=TRUE, values=.01, label=c("b11", "b12","b13","b14", "b11", "b12", "b13","b14"), name="b1" )
 
defSex    <- mxMatrix( type="Full", nrow=1, ncol=8, free=FALSE, labels=c("data.Sex1", "data.Sex1","data.Sex1", "data.Sex1","data.Sex2", "data.Sex2", "data.Sex2","data.Sex2"), name="Sex" )
pathB2     <- mxMatrix( type="Full", nrow=1, ncol=8, free=TRUE, values=.01, label=c("b21", "b22","b23","b24", "b21", "b22", "b23","b24"), name="b2" )
 
# Create Algebra for expected Mean & Threshold Matrices
meanMZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeMZ, name="meanMZ" )
meanDZ    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeDZ, name="meanDZ" )
expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanMZ" )
expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanDZ" )
thinMZ    <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ    <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
inc       <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" )
threMZ    <- mxAlgebra( expression= inc %*% thinMZ, name="threMZ" )
threDZ    <- mxAlgebra( expression= inc %*% thinDZ, name="threDZ" )
 
# Create Algebra for expected Covariance Matrices
covfree <- matrix(TRUE,8,8)
covfree[4,4] <- FALSE
covfree[8,8] <- FALSE
covMZ     <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=covfree, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" )
covDZ     <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=covfree, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" )
 
# Create Data Objects for Multiple Groups
dataMZ    <- mxData( observed=mzDataF, type="raw" )
dataDZ    <- mxData( observed=dzDataF, type="raw" )
 
# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="covMZ", means="expMeanMZ", dimnames=selVars, thresholds="threMZ", threshnames=ordVars )
expDZ     <- mxExpectationNormal( covariance="covDZ", means="expMeanDZ", dimnames=selVars, thresholds="threDZ", threshnames=ordVars )
funML     <- mxFitFunctionML()
 
# Create Model Objects for Multiple Groups
pars      <- list(pathB1,pathB2, inc )
defs      <- list(defAge,defSex)
modelMZ   <- mxModel( name="MZ", pars, defs, meanMZ, expMeanMZ, covMZ, thinMZ, threMZ, dataMZ, expMZ, funML )
modelDZ   <- mxModel( name="DZ", pars, defs, meanDZ, expMeanDZ, covDZ, thinDZ, threDZ, dataDZ, expDZ, funML )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
 
# Create Confidence Interval Objects
ciCor     <- mxCI( c('MZ.covMZ','DZ.covDZ' ))
ciThre    <- mxCI( c('MZ.threMZ','DZ.threDZ' ))
MZZZ       <-mxAlgebra(expression =cov2cor(MZ.covMZ), name="MZZ")
DZZZ       <-mxAlgebra(expression =cov2cor(DZ.covDZ), name="DZZ")
ci2     <- mxCI( c('MZZ','DZZ') )
 
# Build Saturated Model with Confidence Intervals
model     <- mxModel( "twoSATja", pars, modelMZ, modelDZ, multi, MZZZ, DZZZ, ci2 )
 
# ------------------------------------------------------------------------------
# RUN MODEL
 
# Run Saturated Model
fit       <- mxTryHardOrdinal( model, intervals=F, extraTries=11 )
sum       <- summary( fit )
round(fit$output$estimate,4)
fitGofs(fit)
 
# ------------------------------------------------------------------------------
# RUN SUBMODELS
 
 
# Constrain expected Thresholds to be equal across twin order
modelEMTVO  <- mxModel( fit, name="eqThresTwin" )
modelEMTVO <- omxSetParameters(model=fit, labels=c("tauMZ1","tauMZ2","tauDZ1","tauDZ2"), newlabels=c("tauMZ","tauMZ","tauDZ","tauDZ"), name="eqThreshTwin")
modelEMTVO <- omxAssignFirstParameters(modelEMTVO)
fitEMTVO   <- mxRun( modelEMTVO, intervals=F )
mxCompare(fit, subs <- list( modelEMTVO) )
round(fitEMVO$output$estimate,4)
fitGofs(fitEMTVO)
mxCompare (fit, fitEMTVO) 
 
# Constrain expected Means and Variances to be equal across twin order and zygosity
modelEMTVZ <- mxModel(fitEMTVO, name="eqMVarsZyg" )
modelEMTVZ <-  omxSetParameters(model=fitEMTVO, labels=c("tauMZ","tauDZ"), newlabels=c("tau","tau"), name="eqThresZyg")
modelEMTVZ <- omxAssignFirstParameters(modelEMTVZ)
 
fitEMTVZ   <- mxTryHardOrdinal( modelEMTVZ, intervals=F )
mxCompare ( fit, subs <- list( fitEMTVO, fitEMTVZ) )
 
mxCompare (fit, fitEMTVZ) 
 
 
mxEval(cov2cor(DZ.covDZ), fitEMTVZ, T) 
mxEval(cov2cor(MZ.covMZ), fitEMTVZ, T)

Should I equate more parameters? Am I doing something wrong?

Thanks

Fri, 08/17/2018 - 15:53

#10

AdminNeale

Offline

Joined: 03/01/2013 - 14:09

Completely Saturated not ideal

Hi Juan

First, I think you should either constrain the variance of your ordinal phenotype to 1, or fix the residual E to some non-zero value. The second approach is a bit kinder to the optimizer, but the first yields standardized estimates, which is likely the metric that you prefer (correlation not covariance).

Second, a saturated model is usually one where every statistic - mean/threshold/variance/covariance is estimated with a different free parameter. In your case, this model is 'too saturated' because it will yield a different estimated phenotypic correlation for twin 1 and twin 2, and for MZ and DZ twins (4 different phenotypic correlations altogether). What you seem to need is a sort of semi-saturated model where all the variances, covariances, means and thresholds are equated between twin 1 and twin 2, and across MZ and DZ. The covariance across twins is not relevant here, unless you want to control for their non-independence so as to get a better estimate of the standard error of the correlation. If the SE is not of interest, then simply set up 2 columns of data - the continuous and ordinal measures - and obtain the polychoric correlation by estimating the saturated model (and then standardizing as you have with cov2cor).

Controlling for non-independence effects on the standard error of the correlation could be done several ways. You might set up the same block of phenotypic covariances for all 4 blocks (T1, T2, MZ and DZ), but let the off-diagonal blocks differ between MZ and DZ pairs. I'm thinking something like
rbind(cbind(P,M),cbind(M,P)) for the MZs and rbind(cbind(P,D),cbind(D,P)) for the DZs, where the same parameter labels are in place for the P in both groups. Does this make sense to you?

HTH
Mike

Thu, 08/23/2018 - 05:55

(Reply to #10) #11

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Hi Michael,

Thank you so much for your response. Your information is really helpful for me. I know what I have to do now.

However, I am not sure if I know how to do it. I am trying to constrain the variance of the ordinal phenotype to 1. I had thought to do it with mxconstraint but covMZ_4_4, covMZ_8_8, covDZ_4_4, covDZ_8_8 are not present in my model.
I have also tried it with these lines:

matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
matOc <- mxMatrix( type="Full", nrow=1, ncol=nv, free=FALSE, values=oc, name="Oc" )
var1 <- mxConstraint( expression=diag2vec(Oc%&%covMZ)==Unv1, name="Var1" )
var2 <- mxConstraint( expression=diag2vec(Oc%&%covDZ)==Unv1, name="Var2" )

But that does not work either.

I have also changed these lines to do it (fixing the threshold to 0) but I am not sure if it is correct:

thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )

How could I fix the variance for the ordinal phenotype?

Also, I do not know how to set up 2 columns of data and put it in to the model. Should I create 2 different matrices?

Thank you so much again. I have been trying to solve it but I am not able to do it (sorry)

Best wishes.

Main menu

Navigation

You are here

phenotypic correlation

Main menu

User login

Navigation

You are here

phenotypic correlation

Search form