if (FALSE) {
# ===================================================
# = Power to detect correlation of .3 in 200 people =
# ===================================================
# 1 Make some data
tmp = umx_make_raw_from_cov(qm(1, .3| .3, 1), n=2000, varNames= c("X", "Y"), empirical= TRUE)
# 2. Make model of true XY correlation of .3
m1 = umxRAM("corXY", data = tmp,
umxPath("X", with = "Y"),
umxPath(var = c("X", "Y"))
)
# 3. Test power to detect .3 versus 0, with n= 90 subjects
umxPower(m1, "X_with_Y", n= 90)
# ####################
# # Estimating power #
# ####################
#
# method = ncp
# n = 90
# power = 0.83
# sig.level = 0.05
# statistic = LRT
# =================================================
# = Tabulate Power across a range of values of n =
# =================================================
umxPower(m1, "X_with_Y", explore = TRUE)
# =====================================
# = Examples with method = empirical =
# =====================================
# Power to detect r = .3 given n=90
umxPower(m1, "X_with_Y", n = 90, method = "empirical")
# power is .823
# Test using cor.test doing the same thing.
pwr::pwr.r.test(r = .3, n = 90)
# n = 90
# r = 0.3
# sig.level = 0.05
# power = 0.827
# alternative = two.sided
# Power search for detectable effect size, given n = 90
umxPower(m1, "X_with_Y", explore = TRUE)
umxPower(m1, "X_with_Y", n= 90, explore = TRUE)
umxPower(m1, "X_with_Y", n= 90, method = "empirical", explore = TRUE)
data(twinData) # ?twinData from Australian twins.
twinData[, c("ht1", "ht2")] = twinData[, c("ht1", "ht2")] * 10
mzData = twinData[twinData$zygosity %in% "MZFF", ]
dzData = twinData[twinData$zygosity %in% "DZFF", ]
m1 = umxACE(selDVs = "ht", selCovs = "age", sep = "", dzData = dzData, mzData = mzData)
# drop more than 1 path
umxPower(m1, update = c("c_r1c1", "age_b_Var1"), method = 'ncp', n=90, explore = TRUE)
# Specify only 1 parameter (not 'age_b_Var1' and 'c_r1c1' ) to search a parameter:power relationship
# note: Can't use method = "ncp" with search)
umxPower(m1, update = c("c_r1c1", "age_b_Var1"), method = 'empirical', n=90, explore = TRUE)
umxPower(m1, update = c("c_r1c1"), method = 'empirical', n=90, explore = TRUE)
}
Run the code above in your browser using DataLab