# NOT RUN {
# ============================
# = How heritable is height? =
# ============================
require(umx)
data(twinData) # ?twinData from Australian twins.
# Pick the variables
selDVs = c("ht1", "ht2")
mzData <- twinData[twinData$zygosity %in% "MZFF", ]
dzData <- twinData[twinData$zygosity %in% "DZFF", ]
m1 = umxACE(selDVs = selDVs, dzData = dzData, mzData = mzData) # -2ll= 9659
umxSummary(m1, std = FALSE) # unstandardized
# tip: with report = "html", umxSummary can print the table to your browser!
plot(m1)
# ========================================================
# = Evidence for dominance ? (DZ correlation set to .25) =
# ========================================================
m2 = umxACE("ADE", selDVs = selDVs, dzData = dzData, mzData = mzData, dzCr = .25)
umxCompare(m2, m1) # ADE is better
umxSummary(m2, comparison = m1) # nb: though this is ADE, columns are labeled ACE
# ==============================
# = Univariate model of weight =
# ==============================
# Things to note:
# 1. This variable has a large variance, but umx picks good starts.
# 2. umxACE can figure out variable names: provide "sep" and "wt" -> "wt1" "wt2"
# 3. umxACE picks the variables it needs from the data.
# 4. You can use boundDiag to lbound a, c, and e at 0 (prevents mirror-solutions).
m1 = umxACE(selDVs = "wt", dzData = dzData, mzData = mzData, sep = "", boundDiag = 0)
# We can modify this model, dropping shared environment, and see a comparison
m2 = umxModify(m1, update = "c_r1c1", comparison = TRUE)
# =====================================
# = Bivariate height and weight model =
# =====================================
data(twinData)
mzData <- twinData[twinData$zygosity %in% c("MZFF", "MZMM"),]
dzData <- twinData[twinData$zygosity %in% c("DZFF", "DZMM", "DZOS"), ]
mzData <- mzData[1:80,] # quicker run to keep CRAN happy
dzData <- dzData[1:80,]
selDVs = c("ht", "wt") # umx will add suffix (in this case "") + "1" or '2'
m1 = umxACE(selDVs = selDVs, dzData = dzData, mzData = mzData, suffix = '')
umxSummary(m1)
# =========================================================
# = Well done! Now you can make modify twin models in umx =
# =========================================================
# ===================
# = Ordinal example =
# ===================
require(umx)
data(twinData)
# Cut BMI column to form ordinal obesity variables
ordDVs = c("obese1", "obese2")
selDVs = c("obese")
obesityLevels = c('normal', 'overweight', 'obese')
cutPoints <- quantile(twinData[, "bmi1"], probs = c(.5, .2), na.rm = TRUE)
twinData$obese1 <- cut(twinData$bmi1, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels)
twinData$obese2 <- cut(twinData$bmi2, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels)
# Make the ordinal variables into mxFactors (ensure ordered is TRUE, and require levels)
twinData[, ordDVs] <- mxFactor(twinData[, ordDVs], levels = obesityLevels)
mzData <- twinData[twinData$zyg == 1, umx_paste_names(selDVs, "", 1:2)]
dzData <- twinData[twinData$zyg == 3, umx_paste_names(selDVs, "", 1:2)]
mzData <- mzData[1:80,] # just top 80 pairs to run fast
dzData <- dzData[1:80,]
str(mzData) # make sure mz, dz, and t1 and t2 have the same levels!
m1 = umxACE(selDVs = selDVs, dzData = dzData, mzData = mzData, suffix = '')
umxSummary(m1)
# ============================================
# = Bivariate continuous and ordinal example =
# ============================================
data(twinData)
selDVs = c("wt", "obese")
# Cut BMI column to form ordinal obesity variables
ordDVs = c("obese1", "obese2")
obesityLevels = c('normal', 'overweight', 'obese')
cutPoints <- quantile(twinData[, "bmi1"], probs = c(.5, .2), na.rm = TRUE)
twinData$obese1 <- cut(twinData$bmi1, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels)
twinData$obese2 <- cut(twinData$bmi2, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels)
# Make the ordinal variables into mxFactors (ensure ordered is TRUE, and require levels)
twinData[, ordDVs] <- mxFactor(twinData[, ordDVs], levels = obesityLevels)
mzData <- twinData[twinData$zyg == 1,] # umxACE can trim out unused variables on its own
dzData <- twinData[twinData$zyg == 3,]
mzData <- mzData[1:80,] # just top 80 so example runs in a couple of secs
dzData <- dzData[1:80,]
m1 = umxACE(selDVs = selDVs, dzData = dzData, mzData = mzData, suffix = '')
plot(m1)
# =======================================
# = Mixed continuous and binary example =
# =======================================
require(umx)
data(twinData)
# Cut to form category of 20% obese subjects
# and make into mxFactors (ensure ordered is TRUE, and require levels)
cutPoints <- quantile(twinData[, "bmi1"], probs = .2, na.rm = TRUE)
obesityLevels = c('normal', 'obese')
twinData$obese1 <- cut(twinData$bmi1, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels)
twinData$obese2 <- cut(twinData$bmi2, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels)
ordDVs = c("obese1", "obese2")
twinData[, ordDVs] <- mxFactor(twinData[, ordDVs], levels = obesityLevels)
selDVs = c("wt", "obese")
mzData <- twinData[twinData$zygosity %in% "MZFF", umx_paste_names(selDVs, "", 1:2)]
dzData <- twinData[twinData$zygosity %in% "DZFF", umx_paste_names(selDVs, "", 1:2)]
# }
# NOT RUN {
m1 = umxACE(selDVs = selDVs, dzData = dzData, mzData = mzData, suffix = '')
umxSummary(m1)
# }
# NOT RUN {
# ===================================
# Example with covariance data only =
# ===================================
require(umx)
data(twinData)
selDVs = c("wt1", "wt2")
mz = cov(twinData[twinData$zyg == 1, selDVs], use = "complete")
dz = cov(twinData[twinData$zyg == 3, selDVs], use = "complete")
m1 = umxACE(selDVs = selDVs, dzData = dz, mzData = mz, numObsDZ=569, numObsMZ=351)
umxSummary(m1)
plot(m1)
# }
Run the code above in your browser using DataLab