tbates / umx

Making Structural Equation Modeling (SEM) in R quick & powerful
https://tbates.github.io/
44 stars 17 forks source link

umxACE with non-twin groups #108

Closed okeep closed 1 year ago

okeep commented 4 years ago

I was wondering if it would be possible to use the umxACE functions with groups other than MZ and DZ twins? As the function is currently designed it doesn't look like it can accommodate any groups except those two (e.g., adopted siblings, half-siblings, cousins etc.) which is a little unfortunate. Relatedly, would it be possible to include more than two kinship groups overall?

mcneale commented 4 years ago

It's possible to fit models to other types of relatives with OpenMx, which is the backbone of umx. As an intermediate step, you could modify the OpenMx model that umx has generated for you. Based on the help page for umxACE (though I would recommend umxACEv instead):

data(twinData) # ?twinData from Australian twins.
# Pick the variables
# 1. Height has a tiny variance, and this makes solution finding very hard.
# We'll scale height up by 10x to make the Optimizer's task easier.
twinData[, c("ht1", "ht2")] = twinData[, c("ht1", "ht2")] * 10

# 2. umxACE picks the variables it needs from the data.
mzData = twinData[twinData$zygosity %in% "MZFF", ]
dzData = twinData[twinData$zygosity %in% "DZFF", ]

# 3. umxACE can figure out variable names using sep: 
#    e.g. selVars = "wt" + sep= "_T" -> "wt_T1" "wt_T2"
m1 = umxACE(selDVs = "ht", sep = "", dzData = dzData, mzData = mzData)
# Grab the DZ model
halfSibs <- m1$DZ
# Replace DZ data with halfSib data
halfSibs$data <- mxData(mySibData, type='raw')
mxRename(halfSibs,"Half Sibs")
# Add half sibs to m1
m2 <- mxModel(m1, halfSibs)
# Tell the multiple group function that there's another group by overwriting the fit function
m2$top <- mxModel(m1$top,mxFitFunctionMultigroup(c("MZ","DZ","HalfSibs ")))
# Fit the model
summary(m2Run <- mxRun(m2))
tbates commented 4 years ago

Thanks @mcneale !

There are many variations like twins + half-sibs. The notion has been that anyone with the investment in raw data will invest in adding back to umx model extensions that handle their data.

In this case, something like a umxAddHalfSibs() helper could implement this "add group, add expcov, switchdata and expcov" workflow for umxACE, umxACEv, umxCP etc. (they all have the same pattern of algebras neatly stored in top, groups for each dataset).

Adding new groups as @mcneale neatly outlines, you would need also to be careful to use the correct algebra in your additional relatedness groups. The DZ group has:

m1$top$expCovDZ
$formula:  rbind(cbind(ACE, hAC), cbind(hAC, ACE)) 

But (sharing 1/2 one parent + none of another parent = .25 of their variable DNA) you would need to modify the between-twin additive genetic covariance to be .25 A+C right?

Need to be careful that the age covariation was handled correctly (as half sibs don't share a birthday)

So at minimum, for half sibs you'd need to modify model$top to include an algebra for the expected covariance expCovHS, and use this as top.expCovHS to set the expected covariance in your new half sib group. So for top, something like this:

newTop = mxModel(m1$top, 
   umxMatrix("quarter", "Full", 1, 1, values= .25, free=FALSE),
   mxAlgebra(name= "expCovHS", rbind(
         cbind(ACE, quarter %x% A +  C ),
         cbind(quarter %x% A +  C, ACE)
   )
)

m1 = mxModel(m1, newTop) # replace old top with new top

note too that many of the safeguards of umx are stripped away in this approach like whether the C matrix was being used to model dominance... (thousands of lines of umx are checking that all requests and settings are kept coherent and valid)

okeep commented 4 years ago

Thanks, I think your answers @mcneale and @tbates will help.