Skip to content

Commit b54d8e0

Browse files
tbatestbates
committed
better var names in doc data
more plot and summary work
1 parent 5db136b commit b54d8e0

File tree

7 files changed

+186
-138
lines changed

7 files changed

+186
-138
lines changed

‎R/datasets.R

Lines changed: 17 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -17,23 +17,23 @@
1717
# ==================================================
1818
#' Twin data for Direction of causation modelling
1919
#'
20-
#' A dataset containing indicators for two traits `a` and `b`, each measureed in MZ and DZ twins.
20+
#' A dataset containing indicators for two traits `varA` and `varB`, each measureed in MZ and DZ twins.
2121
#'
22-
#' It's designed for testing hypotheses about whether `a` causes `b`, `b` causes `a`, both cause each other.
22+
#' It is designed to show off [umxDoC()] testing the hypothesis `varA` causes `varB`, `varB` causes `varA`, both cause each other.
2323
#'
2424
#' * *zygosity* "MZFF", "DZFF", "MZMM", or "DZMM"
25-
#' * *a1_T1* Twin one's manifest 1 for trait a
26-
#' * *a2_T1* Twin one's manifest 2 for trait a
27-
#' * *a3_T1* Twin one's manifest 3 for trait a
28-
#' * *b1_T1* Twin one's manifest 1 for trait b
29-
#' * *b2_T1* Twin one's manifest 2 for trait b
30-
#' * *b3_T1* Twin one's manifest 3 for trait b
31-
#' * *a1_T2* Twin two's manifest 1 for trait a
32-
#' * *a2_T2* Twin two's manifest 2 for trait a
33-
#' * *a3_T2* Twin two's manifest 3 for trait a
34-
#' * *b1_T2* Twin two's manifest 1 for trait b
35-
#' * *b2_T2* Twin two's manifest 2 for trait b
36-
#' * *b3_T2* Twin two's manifest 3 for trait b
25+
#' * *varA1_T1* Twin one's manifest 1 for varA
26+
#' * *varA2_T1* Twin one's manifest 2 for varA
27+
#' * *varA3_T1* Twin one's manifest 3 for varA
28+
#' * *varB1_T1* Twin one's manifest 1 for varB
29+
#' * *varB2_T1* Twin one's manifest 2 for varB
30+
#' * *varB3_T1* Twin one's manifest 3 for varB
31+
#' * *varA1_T2* Twin two's manifest 1 for varA
32+
#' * *varA2_T2* Twin two's manifest 2 for varA
33+
#' * *varA3_T2* Twin two's manifest 3 for varA
34+
#' * *varB1_T2* Twin two's manifest 1 for varB
35+
#' * *varB2_T2* Twin two's manifest 2 for varB
36+
#' * *varB3_T2* Twin two's manifest 3 for varB
3737
#'
3838
#' @docType data
3939
#' @keywords datasets
@@ -52,56 +52,14 @@
5252
#' dzData = subset(docData, zygosity %in% c("DZFF", "DZMM"))
5353
#' par(mfrow = c(1, 2)) # 1 rows and 3 columns
5454
#' plot(a1_T2 ~a1_T1, ylim = c(-4, 4), data = mzData, main="MZ")
55-
#' tmp = round(cor.test(~a1_T1 + a1_T2, data = mzData)$estimate, 2)
55+
#' tmp = round(cor.test(~varA1_T1 + varA1_T2, data = mzData)$estimate, 2)
5656
#' text(x=-4, y=3, labels = paste0("r = ", tmp))
57-
#' plot(a1_T2 ~a1_T1, ylim = c(-4, 4), data = dzData, main="DZ")
58-
#' tmp = round(cor.test(~a1_T1 + a1_T2, data = dzData)$estimate, 2)
57+
#' plot(varA1_T2 ~varA1_T1, ylim = c(-4, 4), data = dzData, main="DZ")
58+
#' tmp = round(cor.test(~varA1_T1 + varA1_T2, data = dzData)$estimate, 2)
5959
#' text(x=-4, y=3, labels = paste0("r = ", tmp))
6060
#' par(mfrow = c(1, 1)) # back to as it was
6161
NULL
6262

63-
# How I coded this data from the Boulder example
64-
#
65-
# GFF = read.table("~/bin/umx/data/DHBQ_bs.dat", header = TRUE, sep = "\t", as.is = c(T), na.strings = -999)
66-
# x = umx_rename(GFF, old = "zyg2" , replace = "zyg_2grp"); names(x)
67-
# x = umx_rename(x , old = "zyg" , replace = "zyg_6grp"); names(x)
68-
# x = umx_rename(x , grep = "([12bs])$", replace = "_T\\1") ; names(x)
69-
# x$sex_T1 = factor(x$sex_T1, levels = 0:1, labels = c("male", "female"))
70-
# x$sex_T2 = factor(x$sex_T2, levels = 0:1, labels = c("male", "female"))
71-
# x$sex_Tb = factor(x$sex_Tb, levels = 0:1, labels = c("male", "female"))
72-
# x$sex_Ts = factor(x$sex_Ts, levels = 0:1, labels = c("male", "female"))
73-
# x$zyg_6grp = factor(x$zyg_6grp, levels = 1:6, labels = c("MZMM", "DZMM", "MZFF", "DZFF", "DZFM", "DZMF"))
74-
# GFF$zyg_2grp = factor(GFF$zyg_2grp, levels = 1:2, labels = c("MZ", "DZ"))
75-
76-
# GFF = GFF[, c("zyg_6grp", "zyg_2grp", "divorce", "sex_T1", "age_T1", "gff_T1", "fc_T1", "qol_T1", "hap_T1", "sat_T1", "AD_T1", "SOMA_T1", "SOC_T1", "THOU_T1", "sex_T2", "age_T2", "gff_T2", "fc_T2", "qol_T2", "hap_T2", "sat_T2", "AD_T2", "SOMA_T2", "SOC_T2", "THOU_T2", "sex_Tb", "age_Tb", "gff_Tb", "fc_Tb","qol_Tb", "hap_Tb", "sat_Tb", "AD_Tb","SOMA_Tb","SOC_Tb", "THOU_Tb","sex_Ts", "age_Ts", "gff_Ts", "fc_Ts", "qol_Ts", "hap_Ts", "sat_Ts", "AD_Ts","SOMA_Ts","SOC_Ts", "THOU_Ts")]
77-
78-
# save("GFF", file = "GFF.rda")
79-
# system(paste("open ",shQuote(getwd(), type = "csh")))
80-
# update_wordlist get_wordlist(pkg = "~/bin/umx")
81-
82-
# 1. Figure out what things are.
83-
84-
# table(x$sex_Tb) # all 0 so male = 0
85-
# table(x$sex_Ts) # all 1 so female = 1
86-
# umx_aggregate(sex_T2 ~ zyg_6grp, data = x)
87-
# |zyg_6grp |sex_T2 |
88-
# |:-----------|:------------------|
89-
# |1 (n = 448) |male 448; female 0 |
90-
# |2 (n = 389) |male 389; female 0 |
91-
# |3 (n = 668) |male 0; female 668 |
92-
# |4 (n = 484) |male 0; female 484 |
93-
# |5 (n = 504) |male 0; female 504 |
94-
# |6 (n = 407) |male 407; female 0 |
95-
# umx_aggregate(sex_T1 ~ zyg_6grp, data = x)
96-
# |zyg_6grp |sex_T1 |
97-
# |:-----------|:------------------|
98-
# |1 (n = 457) |male 457; female 0 |
99-
# |2 (n = 391) |male 391; female 0 |
100-
# |3 (n = 661) |male 0; female 661 |
101-
# |4 (n = 478) |male 0; female 478 |
102-
# |5 (n = 426) |male 426; female 0 |
103-
# |6 (n = 460) |male 0; female 460 |
104-
10563
# ===================================
10664
# = General Family Functioning data =
10765
# ===================================

‎R/umxDoC.R

Lines changed: 93 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -45,24 +45,25 @@
4545
#' # = Does Rain cause Mud? =
4646
#' # ========================
4747
#'
48+
#' # =======================================
49+
#' # = 2. Define manifests for var 1 and 2 =
50+
#' # =======================================
51+
#' var1 = paste0("varA", 1:3)
52+
#' var2 = paste0("varB", 1:3)
53+
#'
4854
#' # ================
4955
#' # = 1. Load Data =
5056
#' # ================
5157
#' data(docData)
52-
#' mzData = subset(docData, zygosity %in% c("MZFF", "MZMM"))
53-
#' dzData = subset(docData, zygosity %in% c("DZFF", "DZMM"))
54-
#'
55-
#' # =======================================
56-
#' # = 2. Define manifests for var 1 and 2 =
57-
#' # =======================================
58-
#' var1 = paste0("a", 1:3)
59-
#' var2 = paste0("b", 1:3)
58+
#' docData = umx_scale_wide_twin_data(c(var1, var2), docData, sep= "_T")
59+
#' mzData = subset(docData, zygosity %in% c("MZFF", "MZMM"))
60+
#' dzData = subset(docData, zygosity %in% c("DZFF", "DZMM"))
6061
#'
6162
#' # =======================================================
6263
#' # = 2. Make the non-causal (Cholesky) and causal models =
6364
#' # =======================================================
64-
#' Chol= umxDoC(var1= var1, var2= var2, mzData= mzData, dzData= dzData, causal= FALSE)
65-
#' DoC = umxDoC(var1= var1, var2= var2, mzData= mzData, dzData= dzData, causal= TRUE)
65+
#' Chol = umxDoC(var1= var1, var2= var2, mzData= mzData, dzData= dzData, causal= FALSE)
66+
#' DoC = umxDoC(var1= var1, var2= var2, mzData= mzData, dzData= dzData, causal= TRUE)
6667
#'
6768
#' # ================================================
6869
#' # = Make the directional models by modifying DoC =
@@ -83,7 +84,15 @@
8384
#'
8485
#' }
8586
#'
86-
umxDoC <- function(name = "DOC", var1Indicators, var2Indicators, mzData= NULL, dzData= NULL, sep = "_T", causal= TRUE, autoRun = getOption("umx_auto_run"), intervals = FALSE, tryHard = c("no", "yes", "ordinal", "search"), optimizer = NULL) {
87+
umxDoC <- function(name = "DoC", var1Indicators, var2Indicators, mzData= NULL, dzData= NULL, sep = "_T", causal= TRUE, autoRun = getOption("umx_auto_run"), intervals = FALSE, tryHard = c("no", "yes", "ordinal", "search"), optimizer = NULL) {
88+
# TODO: umxDoC add some name checking to avoid variables like "a1"
89+
if(name == "DoC"){
90+
if(causal){
91+
name = "DoC"
92+
} else {
93+
name = "Chol"
94+
}
95+
}
8796
tryHard = match.arg(tryHard)
8897
umx_check(is.logical(causal), "stop", "causal must be TRUE or FALSE")
8998
nSib = 2 # Number of siblings in a twin pair.
@@ -113,18 +122,18 @@ umxDoC <- function(name = "DOC", var1Indicators, var2Indicators, mzData= NULL, d
113122
top = mxModel("top", # (was "ACE")
114123
umxMatrix("dzAr" , "Full", nrow=nLat, ncol=nLat, free=FALSE, values= c(1,.5,.5,1) ), # Heredity Matrix for DZ
115124
umxMatrix("Ones" , "Full", nrow=nLat, ncol=nLat, free=FALSE, values= 1 ), # Unit Matrix - For Com Env and MZ
116-
umxMatrix("Diag1", "Iden", nrow=nSib, ncol=nSib), # Identity matrix (2by2: 1s on diag, 0 off diag)
125+
umxMatrix("Diag1", "Iden", nrow=nSib, ncol=nSib), # Identity matrix (2by2: 1s on diag, 0 off diag)
117126

118127
# Matrices for Cholesky (swapped out after if causal)
119-
umxMatrix("a", type="Lower", nrow=nLat, ncol=nLat, free=TRUE, values= .2), # Genetic effects on Latent Variables
120-
umxMatrix("c", type="Lower", nrow=nLat, ncol=nLat, free=TRUE, values= .2), # Common env effects on Latent Variables
128+
umxMatrix("a", type="Lower", nrow=nLat, ncol=nLat, free= TRUE, values= .2), # Genetic effects on Latent Variables
129+
umxMatrix("c", type="Lower", nrow=nLat, ncol=nLat, free= TRUE, values= .2), # Common env effects on Latent Variables
121130
umxMatrix("e", type="Lower", nrow=nLat, ncol=nLat, free= c(FALSE,TRUE,FALSE), values= 1), # Non-shared env effects on Latent Variables
122131

123132
# 4x4 Matrices for A, C, and E
124133
mxAlgebra(name="A" , Ones %x% (a %*% t(a))),
125134
mxAlgebra(name="Adz", dzAr %x% (a %*% t(a))),
126135
mxAlgebra(name="C" , Ones %x% (c %*% t(c))),
127-
mxAlgebra(name="E" , Diag1 %x% (e %*% t(e))),
136+
mxAlgebra(name="E" , Diag1 %x% (e %*% t(e))),
128137
mxAlgebra(name="Vmz", A + C + E),
129138
mxAlgebra(name="Vdz", Adz + C + E),
130139

@@ -164,24 +173,20 @@ umxDoC <- function(name = "DOC", var1Indicators, var2Indicators, mzData= NULL, d
164173
MZ = mxModel("MZ", mzData, mxExpectationNormal("top.expCovMZ", means= "top.expMean", dimnames= selVars), mxFitFunctionML() )
165174
DZ = mxModel("DZ", dzData, mxExpectationNormal("top.expCovDZ", means= "top.expMean", dimnames= selVars), mxFitFunctionML() )
166175

167-
if(!causal){
168-
# ========================
169-
# = Cholesky-based model =
170-
# ========================
171-
model = mxModel("Chol", top, MZ, DZ, mxFitFunctionMultigroup(c("MZ", "DZ")) )
172-
}else{
176+
if(causal){
173177
# ===================
174178
# = DOC-based model =
175179
# ===================
176180
# Replace lower ace Matrices with diag.
177181
# Because covariance between the traits is "caused", theses matrices are diagonal instead of lower
178182
top = mxModel(top,
179-
umxMatrix("a", "Diag", nrow=nLat, ncol=nLat, free=TRUE, values=0.2), # Genetic effects on Latent Variables
180-
umxMatrix("c", "Diag", nrow=nLat, ncol=nLat, free=TRUE, values=0.2), # Common env effects on Latent Variables
181-
umxMatrix("e", "Diag", nrow=nLat, ncol=nLat, free=FALSE, values=1) # Non-shared env effects on Latent Variables
183+
umxMatrix("a", "Diag", nrow=nLat, ncol=nLat, free=TRUE, values= 0.2), # Genetic effects on Latent Variables
184+
umxMatrix("c", "Diag", nrow=nLat, ncol=nLat, free=TRUE, values= 0.2), # Common env effects on Latent Variables
185+
umxMatrix("e", "Diag", nrow=nLat, ncol=nLat, free=FALSE, values= 1.0) # E@1
182186
)
183-
model = mxModel("DOC", top, MZ, DZ, mxFitFunctionMultigroup(c("MZ", "DZ")) )
184187
}
188+
model = mxModel(name, top, MZ, DZ, mxFitFunctionMultigroup(c("MZ", "DZ")) )
189+
185190
# Factor loading matrix of Intercept and Slope on observed phenotypes
186191
# SDt = mxAlgebra(name= "SDt", solve(sqrt(Diag1 *Rt))) # Standardized deviations (inverse)
187192
model = omxAssignFirstParameters(model)
@@ -219,9 +224,33 @@ umxDoC <- function(name = "DOC", var1Indicators, var2Indicators, mzData= NULL, d
219224
#' @seealso - [umxDoC()], [umxSummary.MxModelDoC()], [umxModify()]
220225
#' @md
221226
#' @examples
222-
#' #
223-
umxPlotDoC <- function(x = NA, means = FALSE, std = TRUE, digits = 2, showFixed = TRUE, file = "name", format = c("current", "graphviz", "DiagrammeR"), SEstyle = FALSE, strip_zero = TRUE, ...) {
224-
message("I'm sorry Dave, no plot for doc yet ;-(")
227+
#'
228+
#' # ================
229+
#' # = 1. Load Data =
230+
#' # ================
231+
#' data(docData)
232+
#' mzData = subset(docData, zygosity %in% c("MZFF", "MZMM"))
233+
#' dzData = subset(docData, zygosity %in% c("DZFF", "DZMM"))
234+
#'
235+
#' # =======================================
236+
#' # = 2. Define manifests for var 1 and 2 =
237+
#' # =======================================
238+
#' var1 = paste0("varA", 1:3)
239+
#' var2 = paste0("varB", 1:3)
240+
#'
241+
#' # =======================================================
242+
#' # = 2. Make the non-causal (Cholesky) and causal models =
243+
#' # =======================================================
244+
#' Chol= umxDoC(var1= var1, var2= var2, mzData= mzData, dzData= dzData, causal= FALSE)
245+
#' DoC = umxDoC(var1= var1, var2= var2, mzData= mzData, dzData= dzData, causal= TRUE)
246+
#'
247+
#' # ================================================
248+
#' # = Make the directional models by modifying DoC =
249+
#' # ================================================
250+
#' a2b = umxModify(DoC, "a2b", free = TRUE, name = "A2B")
251+
#' plot(a2b)
252+
umxPlotDoC <- function(x = NA, means = FALSE, std = TRUE, digits = 2, showFixed = TRUE, file = "name", format = c("current", "graphviz", "DiagrammeR"), SEstyle = FALSE, strip_zero = FALSE, ...) {
253+
message("I'm sorry Dave, no plot for DoC yet ;-(")
225254
# 1. ✓ draw latents
226255
# 2. ✓ draw manifests,
227256
# 3. ✓ draw ace to latents
@@ -234,25 +263,27 @@ umxPlotDoC <- function(x = NA, means = FALSE, std = TRUE, digits = 2, showFixed
234263
umx_check_model(model, "MxModelDoC", callingFn = "umxPlotDoC")
235264

236265
if(std){
237-
model = xmu_standardize_DoC(model)
266+
message("I'm sorry Dave, no std for DoC yet ;-(")
267+
# model = xmu_standardize_DoC(model)
238268
}
239269

240-
nFac = dim(model$top$a_cp$labels)[[1]]
241-
nVar = dim(model$top$as$values)[[1]]
242-
selDVs = dimnames(model$MZ$data$observed)[[2]]
243-
selDVs = selDVs[1:(nVar)]
244-
selDVs = sub("(_T)?[0-9]$", "", selDVs) # trim "_Tn" from end
245-
246-
out = list(str = "", latents = c(), manifests = c())
270+
nFac = dim(model$top$a_cp$labels)[[1]]
271+
nVar = dim(model$top$as$values)[[1]]
272+
selDVs = dimnames(model$MZ$data$observed)[[2]]
273+
selDVs = selDVs[1:(nVar)]
274+
selDVs = sub("(_T)?[0-9]$", "", selDVs) # trim "_Tn" from end
275+
out = list(str = "", latents = c(), manifests = c())
276+
selLat = c("a", "b")
247277

248278
# Process [ace] matrices
249279
# 1. Collect latents
250-
out = xmu_dot_mat2dot(model$top$a, cells = "diag", from = "rows", toLabel = c("a", "b"), fromType = "latent", showFixed = showFixed, p = out)
251-
out = xmu_dot_mat2dot(model$top$c, cells = "diag", from = "rows", toLabel = c("a", "b"), fromType = "latent", showFixed = showFixed, p = out)
252-
out = xmu_dot_mat2dot(model$top$e, cells = "diag", from = "rows", toLabel = c("a", "b"), fromType = "latent", showFixed = showFixed, p = out)
280+
out = xmu_dot_mat2dot(model$top$a, cells = "diag", from = "rows", toLabel = selLat, fromType = "latent", showFixed = showFixed, p = out)
281+
out = xmu_dot_mat2dot(model$top$c, cells = "diag", from = "rows", toLabel = selLat, fromType = "latent", showFixed = showFixed, p = out)
282+
out = xmu_dot_mat2dot(model$top$e, cells = "diag", from = "rows", toLabel = selLat, fromType = "latent", showFixed = showFixed, p = out)
253283

254284
# 2. Process "FacLoad" nVar * nFac matrix: latents into common paths.
255-
out = xmu_dot_mat2dot(model$top$FacLoad, cells= "any", toLabel= selDVs, from= "cols", fromLabel= c("a", "b"), fromType= "latent", showFixed = showFixed, p= out)
285+
out = xmu_dot_mat2dot(model$top$FacLoad, cells= "any", toLabel= selDVs, from= "cols", fromLabel= selLat, fromType= "latent", showFixed = showFixed, p= out)
286+
256287
# 3. Process "as" matrix
257288
out = xmu_dot_mat2dot(model$top$as, cells = "any", toLabel = selDVs, from = "rows", fromType = "latent", showFixed = showFixed, p = out)
258289
out = xmu_dot_mat2dot(model$top$cs, cells = "any", toLabel = selDVs, from = "rows", fromType = "latent", showFixed = showFixed, p = out)
@@ -264,7 +295,7 @@ umxPlotDoC <- function(x = NA, means = FALSE, std = TRUE, digits = 2, showFixed
264295
# [,1] [,2]
265296
# [1,] "a2a" "b2a"
266297
# [2,] "a2b" "b2b"
267-
298+
out = xmu_dot_mat2dot(model$top$beta, cells = "any", toLabel = selLat, from = "cols", fromType = "latent", showFixed = showFixed, p = out, fromLabel=selLat)
268299
# Process "expMean" 1 * nVar matrix
269300
if(means){
270301
# from = "one"; target = selDVs[c]
@@ -314,13 +345,30 @@ plot.MxModelDoC <- umxPlotDoC
314345
#' @seealso - [umxDoC()], [plot.MxModelDoC()], [umxModify()], [umxCP()], [plot()], [umxSummary()] work for IP, CP, GxE, SAT, and ACE models.
315346
#' @md
316347
#' @examples
317-
#' require(umx)
318-
#' data(twinData)
348+
#' # ================
349+
#' # = 1. Load Data =
350+
#' # ================
319351
#' umx_set_auto_plot(FALSE) # turn off autoplotting for CRAN
352+
#' data(docData)
320353
#' mzData = subset(docData, zygosity %in% c("MZFF", "MZMM"))
321354
#' dzData = subset(docData, zygosity %in% c("DZFF", "DZMM"))
322-
#' DoC = umxDoC(var1= paste0("a", 1:3), var2 = paste0("b", 1:3),
323-
#' mzData= mzData, dzData= dzData, causal= TRUE)
355+
#'
356+
#' # =======================================
357+
#' # = 2. Define manifests for var 1 and 2 =
358+
#' # =======================================
359+
#' var1 = paste0("varA", 1:3)
360+
#' var2 = paste0("varB", 1:3)
361+
#'
362+
#' # =======================================================
363+
#' # = 2. Make the non-causal (Cholesky) and causal models =
364+
#' # =======================================================
365+
#' Chol= umxDoC(var1= var1, var2= var2, mzData= mzData, dzData= dzData, causal= FALSE)
366+
#' DoC = umxDoC(var1= var1, var2= var2, mzData= mzData, dzData= dzData, causal= TRUE)
367+
#'
368+
#' # ================================================
369+
#' # = Make the directional models by modifying DoC =
370+
#' # ================================================
371+
#' A2B = umxModify(DoC, "a2b", free = TRUE, name = "A2B")
324372
#' A2B = umxModify(DoC, "a2b", free = TRUE, name = "A2B", comp=TRUE)
325373
#' B2A = umxModify(DoC, "b2a", free = TRUE, name = "B2A", comp=TRUE)
326374
#' umxCompare(B2A, A2B)

‎data/docData.rda

4 Bytes
Binary file not shown.

‎man/docData.Rd

Lines changed: 17 additions & 17 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

0 commit comments

Comments
 (0)