You are here

Adding covariates to Twin Ace model

38 posts / 0 new
Last post
Sabha's picture
Offline
Joined: 05/18/2010 - 15:46
Adding covariates to Twin Ace model

Hi ,

I am new to Twin data analysis and would like to get help on how to adjust for covariates ( Age and Gender) in TWIN ACE model to estimate the genetic and environmental variance from the total molecular phenotypic variance . What modifications are needed in the script UnivariateTwinAnalysis_MatrixRaw.R for age and gender adjusted estimation of A, C and E .

Any insight would be of great help.

Sanchita

trinewa's picture
Offline
Joined: 11/30/2009 - 06:35
I have the same question.

I have the same question. Anyone who can answer this?

neale's picture
Offline
Joined: 07/31/2009 - 15:14
It turns out that Hermine & I

It turns out that Hermine & I were both working on responding to this. Her script will probably be better than mine. Various pieces are available in the documentation (the ACE model script and the definition variable example script) and also with a moderated ACE model (for GxE interaction) in which the moderator was regressed out of the observed variable. To learn OpenMx scripting it would be a GREAT exercise for to try to work from these basic scripts and modify them yourself. There are a few bits and pieces that may not be immediately obvious, of which the most significant is that the algebra for computing the moderated means must appear in the mxModel()s that contain the data and the mxFIMLObjective() function calls.

Anyway, to get you going, let's step through it. Suppose we start with the ACE Model in matrix specification from the OpenMx documentation (you may want to try modifying the path specification version as an exercise - there is no better way to learn than to try stuff out).

http://openmx.psyc.virginia.edu/docs/OpenMx/latest/GeneticEpi_Matrix.html#ace-model-a-twin-analysis

I altered it up a bit to make it a bit more readable (to me), attached as UnivariateTwinAnalysisModeratedMatrix-1.R

The remaining modifications are:
i) simulate a continuous age variable and a binary sex variable for each twin in each zygosity group
ii) glue these fake data together with the observed variables into a data frame for each zygosity group
iii) set up a matrix in the twinACE mxModel part that will be the beta regression weights on each covariate (two in this case, we are assuming that the same regression exists for twin1 and twin2, and in both MZ and DZ groups)
iv) patch in the algebra to multiply the regression weights by the individual's

These are in the attached UnivariateTwinAnalysisModeratedMatrix-3.R

Now, I have not tested that these scripts actually work as advertised. They look like they do the right thing, but a simulation exercise, say using mvrnorm(), would be a good way to validate them. Another exercise for the OpenMx student :) (or :( depending on your perspective).

trinewa's picture
Offline
Joined: 11/30/2009 - 06:35
Thanks! I have tried to build

Thanks! I have tried to build the above explanaition into my script, but I don't seem to be able to work it out from the point that I stand right now.
I am working on a tutorial script from an OpenMx course, and in the data preparation script the following procedure is written like this:

No Age or sex data available, thus I'll simulate it

Age <- rnorm(n=length(IQ1),mean=30,sd=4)
hist(Age)
IQ1.Age <-lm(IQ1~Age)
summary(IQ1.Age)

Note: If there is a significant effect, save residuals to use in OpenMx

IQ1.res <- IQ1.Age$residuals
hist(IQ1.res)

Now, I do have age and sex data in my data file.
I know there are significant sex effects on the IQ1 and IQ2 scores. How do I use the IQ1.res variable into the saturated script for testing differences in means and variances across twins and zygosity groups? And can I utilize the same IQ1.res variable as a covariate in the univariate ACE script?

neale's picture
Offline
Joined: 07/31/2009 - 15:14
Um, with continuous data one

Um, with continuous data one can go either way - i) regress out the covariates first, or ii) specify them in the model. However, it is redundant to do both.

Under i) one would want to use IQ1.res as the observed variable, and don't add age as a covariate (sex perhaps yes).
Under ii) one would skip the lm() piece and just put IQ in the ACE model as observed variable and have age/sex as covariates

If one is interested in judging the statistical significance of the effects of the covariates, it is necessary to use some GEE type of correction to account for correlated observations if using lm() [I am not expert at doing this because I use the ML approach, i.e., ii) instead which corrects for the non-independence of twin1 and twin2 observations by treating the pair as a single data vector]

I'm not sure I follow your question, though.

trinewa's picture
Offline
Joined: 11/30/2009 - 06:35
Once again, thank you! I see

Once again, thank you! I see where I was confused. I have now managed to run alternative i) according to the script you suggested, it worked just fine.

Batouli's picture
Offline
Joined: 06/24/2010 - 01:34
Dear All, I have been trying

Dear All,
I have been trying to run a univariate ACE model, using OpenMx, in which I would like to have age and sex as covariates. I have used “UnivariateTwinAnalysis_MatrixRaw-3” script, and I changed it a little to fit my own data. (both scripts and data are attached). There is no error and the script is working well, but the beta for both age and sex are non-sense. The outputs are like:
1 aa XX 1 1 0.08062441 3.177188e-02
2 cc YY 1 1 -0.04613104 5.783706e-02
3 ee ZZ 1 1 -0.07737747 7.225128e-03
4 mean expMean 1 1 0.66396822 1.638469e-02
5 betaAge beta 1 1 0.10000000 1.078024e+22
6 betaSex beta 1 2 0.10000000 1.078024e+22

The Beta of 0.1 is not changed at all, with any kind of change you apply to your input data. Could anybody tell me please what the problem of my script or my data is? I have only used same-sex MZ and DZ, but I do not think if that might be the source of error in estimating sex-beta.
Thank you in advance,
Amir

neale's picture
Offline
Joined: 07/31/2009 - 15:14
data.Variable not dataFileName.Variable

The optimizer did not get access to the definition variables, because they were labeled as "mzData.ageT1mz" whereas they should have been labeled as data.ageT1mz - the data. tag is critical for OpenMx to realize that you mean a definition variable. Hence your parameter estimates were stuck at their initial values.

mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
Err, this smells like a bug.

Err, this smells like a bug. At least it should have complained if it didn't recognize "mzData.ageT1mz". I'll look into this.

neale's picture
Offline
Joined: 07/31/2009 - 15:14
Not a bug from where I sit

This was expected behavior IMO. mzData.blah is a legal parameter label, just like any other (except it has a . in it which is unusual but not illegal afaik) so the optimizer took it to be a parameter label (fixed as it happens) not a definition variable. Possibly, there is a design flaw; if we want to reject parameter labels that have a . in them, except those that are correctly referencing definition variables, then I guess an error could be thrown and a suitable error message provided: Aha! It looks like you were trying to reference a definition variable because your parameter label has a '.' in it. However, the bit before the '.' doesn't say 'data' or the bit after it doesn't reference a variable in the data frame supplied to the mxData() function call. So please use another label without a '.' in it if you want an ordinary parameter label, or correct the bits before and after the '.' so that they reference a definition variable. Or just go down the pub and have a drink to drown your sorrows...

mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
No, we had agreed that "." is

No, we had agreed that "." is a reserved character and can only be used in specific contexts. The following are all legal uses of the dot character:

  • "data.ageT1mz"
  • "MZ.data.ageT1mz"
  • "model2.B[1,1]"

Apparently we have an error condition that is falling through the cracks, when the "." is used and is neither a valid definition variable nor a valid matrix (row,col) reference.

Batouli's picture
Offline
Joined: 06/24/2010 - 01:34
Thank you very much dear Mike

Thank you very much dear Mike and Michael
Now my script is working. The problem was that the name of the age column in the MZ group (and DZ) was “age1” not “ageT1mz”, so it did not know what “ageT1mz” is. I think it uses the “header” of the input data, not the ones allocated later.
Thank you again for your help.
Regards,
Amir

mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
Two comments: You can use

Two comments:

<

ol>

  • You can use names(dzData) to determine the column names of a data frame. The definition variables need to match up with the column names.
  • Make sure your definition variables are of the form "data.age1" instead of "mzData.age1". The former is a valid definition variable, the latter is not.
  • Batouli's picture
    Offline
    Joined: 06/24/2010 - 01:34
    Thank you again dear Michael.

    Thank you again dear Michael. Yes, you are right. My script worked when it was "data.age1" not the "mzData.age1". I appreciate your help again. Regards, Amir

    mspiegel's picture
    Offline
    Joined: 07/31/2009 - 15:24
    The next binary release will

    The next binary release will report the following error for this script,

    Error: The reference 'mzData.ageT1mz' is illegal because it contains the '.' character in mxMatrix(type = "Full", nrow = 2, ncol = 2, free = F, label = c("mzData.ageT1mz", "mzData.sexT1mz", "mzData.ageT2mz", "mzData.sexT2mz"), name = "MZDefVars") . To write a definition variable use 'data.ageT1mz'

    trinewa's picture
    Offline
    Joined: 11/30/2009 - 06:35
    I have run into trouble

    I have run into trouble trying to enter 'sex ' as covariate into a multivariate ADE Cholesky script (three informants), based on the principles given in the univariate script UnivariateTwinAnalysis_MatrixRaw-3.R above.

    I keep getting the error message :Error in mxModel("ADE", mxMatrix(type = "Lower", nrow = nvar, ncol = nvar, :
    argument is missing, with no default

    Very grateful for ideas to what is wrong in the following script:

    mzData <- data.frame(subset(RSmv_data, ZYGEND==1, c(T1MRS, T1FRS, T1TRS, T2MRS, T2FRS, T2TRS, ASEX,BSEX)))
    dzData <- data.frame(subset(RSmv_data, ZYGEND==2, c(T1MRS, T1FRS, T1TRS, T2MRS, T2FRS, T2TRS, ASEX,BSEX)))

    ADE_Cholesky2_Model <- mxModel("ADE_Cholesky2",
    mxModel("ADE",
    # Matrices a, d, and e to store a, d, and e path coefficients
    mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.4, name="a" ),
    mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE,values=.4, name="d" ),
    mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.4, name="e" ),
    # Matrices A, D, and E compute variance components
    mxAlgebra( a %% t(a), name="A" ),
    mxAlgebra( d %
    % t(d), name="D" ),
    mxAlgebra( e %% t(e), name="E" ),
    # Algebra to compute total variances and standard deviations (diagonal only)
    mxAlgebra( A+D+E, name="V" ),
    mxMatrix( type="Iden", nrow=nvar, ncol=nvar, name="I"),
    mxAlgebra( solve(sqrt(I
    V)), name="iSD"),

    # Matrix &amp; Algebra for expected means vector
        mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= 3.6, name="Mean" ),
        mxAlgebra(  cbind(Mean,Mean, Mean, Mean, Mean, Mean), name="expMean"),
    

    Adapted from Neale's univ script: Declare a matrix for the definition variable regression parameters, called beta

    mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= 0, label=c("betaSex"), name="beta"),
    
    # Algebra for expected variance/covariance matrix in MZ
        mxAlgebra(  rbind  ( cbind(A+D+E , A+D),
                                        cbind(A+D   , A+D+E)), name="expCovMZ" ),
    # Algebra for expected variance/covariance matrix in DZ
        mxAlgebra(  rbind  ( cbind(A+D+E     , 0.5%x%A+0.25%x%D),
                                        cbind(0.5%x%A+0.25%x%D, A+D+E)),  name="expCovDZ" ), 
    # From Neale univ.: Algebra for making the means a function of the definition variable sex
      mxModel("MZ", mxData( observed=mzData, type="raw" ),
        mxMatrix( type="Full", nrow=1, ncol=2, free=F, label=c("data.ASEX","data.BSEX"), name="MZDefVars"),
        mxAlgebra( expression=ADE.expMean + ADE.beta %*% MZDefVars, name="expMeanMZ"),
        mxFIMLObjective( covariance="ADE.expCovMZ", means="expMeanMZ", dimnames=selVars )   ),
    mxModel("DZ", mxData( observed=dzData, type="raw" ),
        mxMatrix( type="Full", nrow=1, ncol=2, free=F, label=c("data.ASEX","data.BSEX"), name="DZDefVars"),
        mxAlgebra( expression=ADE.expMean + ADE.beta %*% DZDefVars, name="expMeanDZ"),
        mxFIMLObjective( covariance="ADE.expCovDZ", means="expMeanDZ", dimnames=selVars )   ),
    
    ),
    mxModel("MZ",
        mxData( observed=mzData, type="raw" ),
        mxFIMLObjective( covariance="ADE.expCovMZ", means="ADE.expMean", dimnames=selVars )
    ),
    mxModel("DZ", 
        mxData( observed=dzData, type="raw" ),
        mxFIMLObjective( covariance="ADE.expCovDZ", means="ADE.expMean", dimnames=selVars ) 
    ),
    mxAlgebra( MZ.objective + DZ.objective, name="modelfit" ),
    mxAlgebraObjective("modelfit")
    

    )
    ADE_Cholesky2_Fit <- mxRun(ADE_Cholesky2_Model)
    ADE_Cholesky2_Summ <- summary(ADE_Cholesky2_Fit)
    ADE_Cholesky2_Summ

    tbates's picture
    Offline
    Joined: 07/31/2009 - 14:25
    that error nearly always

    that error nearly always means a trailing ","

    in this case it is

    mxFIMLObjective( covariance="ADE.expCovDZ", means="expMeanDZ", dimnames=selVars )   ),
    
    ),
    

    delete the comma and your model enters fine.

    But you also seem to have duplicated groups for MZ and DZ: Once inside modelADE and once inside ADE_Cholesky2

    trinewa's picture
    Offline
    Joined: 11/30/2009 - 06:35
    I have deleted the comma and

    I have deleted the comma and tried to keep the groups for MZ and DZ either inside modelADE or inside the ADE_Cholesky2. Either way I still get error messages. Below you see my atttempt to write the script connecting the MZ and DZ groups within the modelADE. This produces the following error message:
    Running ADE_Cholesky2
    Error: Trying to evaluate 'MZ.expMeanMZ' in model 'ADE_Cholesky2' generated the error message: non-conformable arrays

    I have attached an example data file - SO happy if you could help me out and make this script work.

    nvar <- 3 # number of variables

    tnvar <- 6 # number of variablesmax family size
    nlower <-tnvar
    (tnvar+1)/2 # number of free elements in a lower matrix tnvar*tnvar

    Vars <- c('Mex','Fex','Tex')
    selVars <- paste("T",c(rep(1,nvar),rep(2,nvar)),Vars,sep="")
    print(selVars)

    mzData <- data.frame(subset(Exmv_data, Zyg==1, c(T1Mex, T1Fex, T1Tex, T2Mex, T2Fex, T2Tex, T1sex,T2sex)))
    dzData <- data.frame(subset(Exmv_data, Zyg==2, c(T1Mex, T1Fex, T1Tex, T2Mex, T2Fex, T2Tex, T1sex,T2sex)))

    ADE_Cholesky2_Model <- mxModel("ADE_Cholesky2",
    mxModel("ADE",
    # Matrices a, d, and e to store a, d, and e path coefficients
    mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.4, name="a" ),
    mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE,values=.4, name="d" ),
    mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.4, name="e" ),
    # Matrices A, D, and E compute variance components
    mxAlgebra( a %% t(a), name="A" ),
    mxAlgebra( d %
    % t(d), name="D" ),
    mxAlgebra( e %% t(e), name="E" ),
    # Algebra to compute total variances and standard deviations (diagonal only)
    mxAlgebra( A+D+E, name="V" ),
    mxMatrix( type="Iden", nrow=nvar, ncol=nvar, name="I"),
    mxAlgebra( solve(sqrt(I
    V)), name="iSD"),

    # Matrix &amp; Algebra for expected means vector
        mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= 3.6, name="Mean" ),
        mxAlgebra(  cbind(Mean,Mean, Mean), name="expMean"),
    

    From univariate script: Declare a matrix for the definition variable regression parameters, called beta

    mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= 0, label=c("betaSex"), name="beta"),
    
    # Algebra for expected variance/covariance matrix in MZ
        mxAlgebra(  rbind  ( cbind(A+D+E , A+D),
                                        cbind(A+D   , A+D+E)), name="expCovMZ" ),
    # Algebra for expected variance/covariance matrix in DZ
        mxAlgebra(  rbind  ( cbind(A+D+E     , 0.5%x%A+0.25%x%D),
                                        cbind(0.5%x%A+0.25%x%D, A+D+E)),  name="expCovDZ" ), 
    # Algebra for making the means a function of the definition variable sex
      mxModel("MZ", mxData( observed=mzData, type="raw" ),
        mxMatrix( type="Full", nrow=1, ncol=2, free=F, label=c("data.T1sex","data.T1sex"), name="MZDefVars"),
        mxAlgebra( expression=ADE.expMean + ADE.beta %*% MZDefVars, name="expMeanMZ"),
        mxFIMLObjective( covariance="ADE.expCovMZ", means="expMeanMZ", dimnames=selVars )
    

    ),
    mxModel("DZ", mxData( observed=dzData, type="raw" ),
    mxMatrix( type="Full", nrow=1, ncol=2, free=F, label=c("data.T1sex","data.T1sex"), name="DZDefVars"),
    mxAlgebra( expression=ADE.expMean + ADE.beta %*% DZDefVars, name="expMeanDZ"),
    mxFIMLObjective( covariance="ADE.expCovDZ", means="expMeanDZ", dimnames=selVars ))

    ),
    mxAlgebra( MZ.objective + DZ.objective, name="modelfit" ),
    mxAlgebraObjective("modelfit")
    )
    ADE_Cholesky2_Fit <- mxRun(ADE_Cholesky2_Model)
    ADE_Cholesky2_Summ <- summary(ADE_Cholesky2_Fit)
    ADE_Cholesky2_Summ

    mspiegel's picture
    Offline
    Joined: 07/31/2009 - 15:24
    Here's some help in debugging

    Here's some help in debugging a MxAlgebra error. The argument 'compute=TRUE' will instruct mxEval() to attempt to perform the algebra operations. The default behavior of mxEval() is to simply return the current value stored in the algebra. The argument 'show=TRUE' will instruct mxEval() to print out what it is evaluating. You can see below that ADE.expMean is a 1 x 9 matrix, ADE.beta is 1 x 1 matrix, and MZ.MZDefVars is a 1 x 2 matrix.

    > mxEval(MZ.expMeanMZ, ADE_Cholesky2_Model, compute=TRUE, show=TRUE)
    mxEval(ADE.expMean + ADE.beta %*% MZ.MZDefVars, ADE_Cholesky2_Model, TRUE) 
    Error: Trying to evaluate 'MZ.expMeanMZ' in model 'ADE_Cholesky2' 
        generated the error message: non-conformable arrays
    > mxEval(ADE.expMean, ADE_Cholesky2_Model, compute=TRUE, show=TRUE)
    mxEval(cbind(ADE.Mean, ADE.Mean, ADE.Mean), ADE_Cholesky2_Model, TRUE) 
         [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
    [1,]  3.6  3.6  3.6  3.6  3.6  3.6  3.6  3.6  3.6
    > mxEval(ADE.beta, ADE_Cholesky2_Model, compute=TRUE, show=TRUE)
    ADE_Cholesky2_Model[["ADE.beta"]]@values 
         [,1]
    [1,]    0
    > mxEval(MZ.MZDefVars, ADE_Cholesky2_Model, compute=TRUE, show=TRUE)
    ADE_Cholesky2_Model[["MZ.MZDefVars"]]@values 
         [,1] [,2]
    [1,]    0    0
    
    neale's picture
    Offline
    Joined: 07/31/2009 - 15:14
    Hi There are quite a few

    Hi

    There are quite a few issues with the script that stop it from working. Most of these center on the lack of agreement between the dimensions of matrices that you are trying to add or multiply. I note that beta is still a 1x1 matrix but there are three observed variables here and they probably need different beta weights. So if beta was nrow=1 and ncol=nvar it would be a start, but then in the expMean statement you would need something like ADE.beta * MZDefVars, or alternatively just make MZDefVars 1x1 and use MZDefVars %*% ADE.beta

    One of the problems here is the shortage of information about why an algebra is not conformable. That is something OpenMx should improve. Another is that the scripting style makes it difficult to check things out 'along the way'. So if, for example, you were to construct mxMatrices before the mxModel, you could inspect them and make sure they are the right dimensions for computing the algebras.

    Another thing that can really help is to use a pencil & paper and write out the matrix formulae you are trying to evaluate using notation of the form iAj %*% jBk where i and j are the number of rows and columns in A. In this case the matrices are conformable for matrix multiplication because the number of cols in A is the same as the number of rows in B.

    trinewa's picture
    Offline
    Joined: 11/30/2009 - 06:35
    Thanks both for valuable

    Thanks both for valuable hints! In my struggle to get the matrixes into compatible sizes, I have a question concerning the sex variable. Sex can vary between twins in the DOS group. Thus, I cannot use Twin 1 sex as the only basis for setting the sex variable for each pair, nor for each individual in a pair (unlike an AGE variable, e.g.). Am I right in assuming that I will need a means matrix with 6 means (3 variables, and 2 twins), something like a nrow=ntwin and ncol=nvar matrix? And equally with the ADE.beta matrix? Or do I need separate matrixes for each twin in each zygosity group?

    neale's picture
    Offline
    Joined: 07/31/2009 - 15:14
    With a binary variable such

    With a binary variable such as sex, it is reasonably easy to use it directly as a dummy indicator variable, assuming that it's coded 0/1. You are right that you would want a 1x6 mean vector at the end, so it might be obtained by:

    u %x% mu + defsex %x% sexbetas

    where defsex is 1x2 containing data.sex1 and data.sex2, and sexbetas is a 1x3 matrix. mu is for the population means, and u is a 1x2 unit vector to replicate mu for twin 1 and twin 2.

    If the variance components are allowed to vary according to sex, I might be tempted to go for a 6 mxModel set up - one "holdall" model that sets up the matrices etc, and 5 mxModels within it which fit the data from each of the 5 zygosity groups. Thus it will be possible to allow for a variety of sex-limitation models. An alternative might be to use the sex definition variables to adjust the ACE or ADE parameters according to sex. This might turn out to be more succinct (only need one data group if you work with zygosity coded 1 or .5 in a similar fashion). However, what I have outlined is more likely consistent with, e.g., Boulder 2010 workshop examples.

    trinewa's picture
    Offline
    Joined: 11/30/2009 - 06:35
    Hi again! What you say at the

    Hi again! What you say at the end of your last reply seems to be connected to an earlier discussion in the general SEM forum: http://openmx.psyc.virginia.edu/thread/250

    Please help me check out if I have understood your answers:

    So far, what I believe I have been trying to do has been what Bates in the earlier discussion calls alternative 1. I am not sure if I can find an example script at the Boulder 2010 workshop that applies to alternative 1, but adding sex as a covariate is what you have been helping me to do so far in the present discussion. Altough the script with your suggested additions semed to run well in the univariate case, I am not sure if it included the dizygotic opposite sex group correctly (I would appreciate your advice on this). I am still working on adapting the covariate script to the multivariate ADE case. I find it hard to know whether the resulting script is Ok, and I would be extremely happy if someone would be willing to share a well tested multivariate script with sex as covariate to validiate mine up against.

    Alternative 2, running the ACE/ADE models with raw data residualized for sex, is equivalent to alternative 1, is that right? This is very simple in that one could run analyses, univariate as well as multivariate, with residualized dependent variables, and no other changes would need to be made in the no covariates scripts. I have tried to save residualized variables for two twins and three informants into a new data file in R, following this procedure:

    T1Mres <- lm(T1M ~ ASEX, na.action = na.exclude)
    T2Mres <- lm(T2M ~ BSEX, na.action = na.exclude)
    T1Fres <- lm(T1F ~ ASEX, na.action = na.exclude)
    T2Fres <- lm(T2F ~ BSEX, na.action = na.exclude)
    T1Tres <- lm(T1T ~ ASEX, na.action = na.exclude)
    T2Tres <- lm(T2T ~ BSEX, na.action = na.exclude)

    T1Mresid <- residuals(T1Mres)
    T2Mresid <- residuals(T2Mres)
    T1Fresid <- residuals(T1Fres)
    T2Fresid <- residuals(T2Fres)
    T1Tresid <- residuals(T1Tres)
    T2Tresid <- residuals(T2Tres)

    frame1 <- data.frame(T1Mresid=T1Mresid, T2Mresid=T2Mresid,T1Fresid=T1Fresid, T2Fresid=T2Fresid,T1Tresid=T1Tresid, T2Tresid=T2Tresid,ZYG=ZYG, ZYGsex=ZYGsex)

    write.table(x=frame1, file="MFTresiduals.DAT",
    quote=FALSE,sep='\t',na='-1.000',row.names=FALSE)

    ResidMFT <- read.table("MFTresidualS.DAT",header=TRUE, na.strings="-1.000")
    summary (ResidMFT)
    names (ResidMFT)

    I read from the Neale/Bates discussion that there is a problem with missing data in the sex variable with this procedure. I really do not get the conclusion of that discussion, though. Thus, my question conserning alternative 2 is: given that alternative 2 seems to be such an easy way out, is there any reason why I should drop this and prefer to go through the (for me)complicated scripting in alternative 1?

    The alternative with 6 mxModels that you suggest in your latest reply, allowing the variance components to vary according to sex, would that be a sex as a moderator model? That is, alternative 3 as mentioned in your discussion with Bates? Then, would it be a good idea to use the univariate heterogeneity script provided at the Boulder 2010 workshops site as a starting point for alternative 3? Are you aware of any example script that covers the alternative you mention here in the multivariate case?

    Sorry about the length of this comment, but I need to get this right.

    tbates's picture
    Offline
    Joined: 07/31/2009 - 14:25
    One critical thing with

    One critical thing with residualising, is that all subjects need to be residualised together, otherwise you are applying different residualisation functions to different groups. You don't want to do the different kinds of twin (either birth order or zygosity) separately.

    So do the residualisation while the data are in long format (one row per person), then move to a wide format (one family per row). ?reshape can help with this.

    Also, not sure from the snippet, but it looks like you are residualizing each sex separately ? If so, that definitely won't do what you want: You need both sexes in the data for residualization by sex to be valid (I'm interpreting things like "T1M" as "twin 1Male")

    trinewa's picture
    Offline
    Joined: 11/30/2009 - 06:35
    Thank you for important

    Thank you for important points. T1M means Twin 1 rated by Mother, så any sex can be represented in that variable.

    trinewa's picture
    Offline
    Joined: 11/30/2009 - 06:35
    Comment deleted, redefined as

    Comment deleted, redefined as new thread.

    martonandko's picture
    Offline
    Joined: 02/19/2015 - 06:55
    Warning messages when usinig openMX v2.0.1

    Thanks for your great UnivariateTwinAnalysis_MatrixRaw-3.R script! Adding a third continuous covariate (openMX_ACE_3cov.R) Running the script I get warning messages:
    1-2: In mxFIMLObjective(covariance = "twinACE.expCovMZ", means = "expMeanMZ", :
    Objective functions like mxFIMLObjective() have been deprecated in favor of expectation and fit functions.
    Please use mxExpectationNormal(covariance= , means = , ...) instead, and add a call to mxFitFunctionML(). See examples at help(mxExpectationNormal)
    3: In mxAlgebraObjective("twin") :
    Objective functions like mxAlgebraObjective() have been deprecated in favor of expectation and fit functions.
    Please use mxFitFunctionAlgebra(algebra = ...). See examples at help(mxFitFunctionAlgebra)

    Are these of any concern? I tried to implement the above suggested modifications (openMX_ACE_3cov_modif.R) and got another warning message:

    > twinACEFit <- mxRun(twinACEModel)
    Running twinACE
    Error: The references 'MZ.fitfunction' and 'DZ.fitfunction' are unknown in the algebra named 'twinACE.twin'
    In addition: Warning messages:
    1: In is.na(x) : is.na() applied to non-(list or vector) of type 'symbol'
    2: In is.na(x) : is.na() applied to non-(list or vector) of type 'symbol'

    And the script won't run. Is there anything I need to add to the script?

    My initial goal would be to get bootstrapped CI for the A, C, E estimates, using this model definition using the umxCI_boot funtion, but it won’t take in the Model.

    mhunter's picture
    Offline
    Joined: 07/31/2009 - 15:26
    No worries

    Note the difference between a WARNING and an ERROR. Warnings mean the code worked, but you might want to check that is worked as you intended. Errors mean the code did not work. The warning message says that the code worked, but there's a newer way to do the thing you're doing.

    See http://openmx.psyc.virginia.edu/2012/11/new-feature-trunk-mxexpect-and-mxfit-functions about the fit function expectation split.

    For your modified code, change

    mxExpectationNormal( covariance="twinACE.expCovDZ", means="expMeanDZ", dimnames=selVars )

    to

    mxExpectationNormal( covariance="twinACE.expCovDZ", means="expMeanDZ", dimnames=selVars ),
    mxFitFunctionML()

    and similarly throughout. Objective functions have been replaced with a pair of functions: expectations which map model matrices to a probability distribution, and fit functions which compare the model expected distribution with the observed data in a single measure of fit.

    Also where you have

    mxAlgebra( expression=MZ.objective + DZ.objective, name="twin" ),
    mxFitFunctionAlgebra("twin")

    you can (but do not have to) replace it with

    mxFitFunctionMultigroup( c("MZ.fitfunction",  "DZ.fitfunction"))
    martonandko's picture
    Offline
    Joined: 02/19/2015 - 06:55
    It's working

    Thank You, it seems to be working fine!

    Cindy.s's picture
    Offline
    Joined: 07/04/2015 - 18:52
    Covariates be added in the expected covariance structure?

    Hi openMX community! I’m trying to work out some classical twin models where I’m interested if some factor (for example: Hypertension) may cause difference in the estimated variance components of ACE models.
    Can you add additional binary and continuous variable as covariates the same way as Sex and Age was added in Neale’s above posted script? Could I add Hypertension/Weight as another covariate and then compare the results with a model where only Age and Sex are covariates based on chi-square test to assess if the A, C and E components are significantly different between the two models? Should the covariates also be included in the expected covariance structure, not only the mean structure if I would like to “cancel out” the effect of a covariate?

    AdminRobK's picture
    Offline
    Joined: 01/24/2014 - 12:15
    Yes, you can condition the

    Yes, you can condition the expected covariance structure on covariates, in a way that is similar to what Professor Neale does in his script. Such a model is called a moderation model, and in this context, the covariate is called a "moderator." The idea is that the path coefficients from the latent A, C, and E factors each reflect a main effect of the latent factor plus an interaction between the latent factor and the moderator.

    You might want to check out this other thread for details. Also, I might be able to dig up an example script.

    Cindy.s's picture
    Offline
    Joined: 07/04/2015 - 18:52
    Methods for adding covariates

    Thank you RobK for your help, I really appreciate it! Just to get it strait:
    Covariates effects can be regressed on the mean structure using linear regression if there is suspicion that they have an effect on the dependent variables value. If I assume a covariate also effects the dependent variables variability, then I need to control for the effect of the moderator on the covariance structure.
    If my moderator is binary then basically I have a multiple group case with 6 groups (GxE model), if discordant MZ pairs are not present, then I have a sex limitation model, but basically they both are the same.
    If my moderator is continuous, then I can add it to the covariance matrix if the correlation is 1 between the pairs, for example age, otherwise it’s easier to go with multivariate models such as bivariate Cholesky to assess the connection between the heritability of two continuous variables.
    I would really appreciate if you had any sample scripts for GxE models, especially ones with all the 6 groups.

    AdminRobK's picture
    Offline
    Joined: 01/24/2014 - 12:15
    Right

    I wouldn't word it exactly the way you do, but you've definitely got the right idea. I have been working recently on an example script for the so-called "bivariate" moderation model, but it's not ready yet.

    zhanglab's picture
    Offline
    Joined: 06/17/2014 - 10:29
    A question about Dr. Neale's code

    Dear OpenMx experts,

    I have a question about Dr. Neale's code in the file "UnivariateTwinAnalysis_MatrixRaw-3.R". Is this code using the ACE-XYZ-M model, which uses age and sex as the moderators (Purcell 2002, Variance Component Models for Gene-Environment Interaction in Twin Analysis)?

    If not, is this code just simply regressing out age and sex using linear regression before classic ACE model fitting?

    Thank you.

    AdminRobK's picture
    Offline
    Joined: 01/24/2014 - 12:15
    The latter

    It's just regressing out age and sex in the model for the means. There's no moderation of biometric variance components in UnivariateTwinAnalysis_MatrixRaw-3.R .

    jhonkevin123's picture
    Offline
    Joined: 11/13/2011 - 04:27
    Helpful Post

    Thanks For Helpful Posts.
    I had the same problem as well.

    Jhon Kevin

    cvanhulle's picture
    Offline
    Joined: 05/07/2010 - 17:28
    covariates with ordinal data

    Hi all,

    I have an ordinal model and I'm trying to include covariates in the means model. I tried the suggestions in this thread, but I keep getting an error - The job for model 'univACEOrd' exited abnormally with the error message: Objective function returned an infinite value.

    If I comment out sections related to the covariates the model fits fine. The script is below. Any suggestions would be very much appreciated.

    Thanks,

    Carol

    univACEOrdModel <- mxModel("univACEOrd",
    mxModel("ACE",
    # Matrices a, c, and e to store a, c, and e path coefficients
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="a" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="c11", name="c" ),
    mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", name="e" ),
    # Matrices A, C, and E compute variance components
    mxAlgebra( expression=a %% t(a), name="A" ),
    mxAlgebra( expression=c %
    % t(c), name="C" ),
    mxAlgebra( expression=e %% t(e), name="E" ),
    mxMatrix( type="Zero", nrow=1, ncol=nv, name="M" ),
    mxAlgebra( expression= cbind(M,M), name="expMean" ),
    mxMatrix(type="Full", nrow=1, ncol=2, free=T, values=0.5, label="thre", name="expThre", dimnames=list('th1',selVars)),
    # Algebra to compute total variances and standard deviations (diagonal only)
    mxAlgebra( expression=A+C+E, name="V" ),
    mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
    mxAlgebra( expression=solve(sqrt(I
    V)), name="sd"),
    mxAlgebra( expression=cbind(A/V,C/V,E/V),name="stndVCs"),
    # confidence interval
    mxCI("ACE.stndVCs"),
    # Constraint on variance of ordinal variables
    mxConstraint( V==I, name="Var1"),
    #Declare matrix of definition variables
    mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=-.15, label="betaSex", name="beta"),

    # Algebra for expected variance/covariance matrix in MZ
        mxAlgebra( expression= rbind  ( cbind(A+C+E , A+C),
                                        cbind(A+C   , A+C+E)), name="expCovMZ" ),
    # Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix
        mxAlgebra( expression= rbind  ( cbind(A+C+E     , 0.5%x%A+C),
                                        cbind(0.5%x%A+C , A+C+E)),  name="expCovDZ" ) 
                                           ),
    mxModel("MZ",
    

    Matrix & Algebra for expected means vector

      mxMatrix(type="Full", nrow=1, ncol=2, free=F, label=c("data.gender1","data.gender2"), name="MZDefVars"),
        mxAlgebra(expression=ACE.expMean+ACE.beta %*% MZDefVars, name="expMeanMZ"),
      mxData( observed=mzData, type="raw" ),
        mxFIMLObjective( covariance="ACE.expCovMZ", means="expMeanMZ", dimnames=selVars, thresholds="ACE.expThre" )
    ),
    mxModel("DZ", 
    

    Matrix & Algebra for expected means vector

      mxMatrix(type="Full", nrow=1, ncol=2, free=F,label=c("data.gender1","data.gender2"), name="DZDefVars"),
      mxAlgebra(expression=ACE.expMean + ACE.beta %*% DZDefVars, name="expMeanDZ"),
        mxData( observed=dzData, type="raw" ),
        mxFIMLObjective( covariance="ACE.expCovDZ", means="expMeanDZ", dimnames=selVars, thresholds="ACE.expThre" )
    ),
    mxAlgebra( expression=MZ.objective + DZ.objective, name="2sumll" ),
    mxAlgebraObjective("2sumll")
    

    )

    cvanhulle's picture
    Offline
    Joined: 05/07/2010 - 17:28
    covariates with ordinal data

    Ok, the script is sort of working.