#' --- #' title: | #' Simulating Complex Data: #' Start with a DAG #' author: "MATH 4939" #' date: | #' January 2026 #' Last update: `r format(Sys.time(), '%B %d %Y %H:%M')` #' geometry: top=.4in,left=.4in,right=.4in,bottom=.45in,paperwidth=8in,paperheight=4.5in #' output: #' bookdown::pdf_document2: #' keep_tex: true #' latex_engine: xelatex #' pandoc_args: [ #' "-V", "classoption=twocolumn" #' ] #' toc: false #' header-includes: #' - \usepackage{pdfpages} #' - \pagenumbering{arabic} #' # - \usepackage{setspace} #' # - \AtBeginDocument{\singlespacing\let\oldcontentsline\contentsline\renewcommand{\contentsline}[4]{\oldcontentsline{#1}{#2}{#3}{\singlespacing #4}}} #' # - \AtEndDocument{\doublespacing} #' # bibliography: ../4939.bib #' --- #' \onecolumn #' \pagenumbering{gobble} #' \pagenumbering{arabic} #' \raggedright #' \tableofcontents #' \listoffigures #' \listoftables #' \large #+ setup, include=FALSE ------ knitr::opts_chunk$set(comment = ' ',fig.height = 3.2,fig.width=5,out.width="3.8in") library(spida2) library(png) library(knitr) #' #' \newpage #' # Setting up a DAG #' \newpage #' **Start with a DAG** #+ dag1, fig.cap='A DAG' include_graphics('simdag1.png') #' \newpage #' **Identify variables (nodes) without parents???** #+ dag3, fig.cap='A DAG with parent-less variables' include_graphics('simdag3.png') #' \newpage #' **Decide on unique variance for each variable** #+ dag4, fig.cap='A DAG with unique independent "error" variances for each variable' include_graphics('simdag4.png') #' \newpage #' **Decide on causal coefficients** #+ dag5, fig.cap='A DAG with causal coefficients' include_graphics('simdag5.png') #' #' # Write a program to generate data #' #' If you plan to vary some parameters, include them among parameters of your function. #' You should probably always make the sample size a parameter. #' gen <- function(N, seed = NULL) { # N: sample size # seed: optional random seed for reproducibility, default: current seed if(!is.null(seed)) set.seed(seed) # create a data frame with one of the parent-less variables dd <- data.frame(z4 = rnorm(N)) # generate remaining variables making sure, of course, to generate # all of a variable's parents before generating each variable itself. # You can use 'within' in base R or 'mutate' in 'dplyr'. dd <- within( dd, { z2 <- 4 * z4 + 0.1 * rnorm(N) z1 <- 3 * z4 + .5 * z2 + 0.1 * rnorm(N) z6 <- rnorm(N) z5 <- z4 + z6 + rnorm(N) z8 <- z6 + rnorm(N) z7 <- rnorm(N) x <- 1 * z7 + .5 * z1 + 0.1 * rnorm(N) z9 <- x + rnorm(N) m1 <- 2 * x + 0.5 * rnorm(N) z3 <- .2 * z2 + rnorm(N) y <- z8 + .4 * z3 + x + .3 * m1 + 0.1 * rnorm(N) } ) return(dd) # not essential to use 'return' but good style # according to the Google style sheet. } #' \newpage #' **Test it** gen(2) gen(3) gen(2, 874432) gen(4, 874432) # first 2 same as previous #' \newpage #' With a huge data set, your estimated parameters will be close #' to the population parameter, i.e. the coefficients you used #' in your simulation. #' #' If we regress y on its parents (immediate ancestors) we should get something #' close to the parameters we used to generate y. #' simdata <- gen(10000) head(simdata) dim(simdata) system.time( summary(lm(y ~ z8 + z3 + x + m1, data = simdata)) %>% print ) #' #' What do you think would happen to the SEs if you were to use a #' sample size of 1,000,000 instead of 10,000? #' \newpage #' **What models should give an unbiased estimate of the causal effect of X?** #' #' Here are 4 of them. Can you list the others? #' fit1 <- lm(y ~ x + z1, simdata) fit2 <- lm(y ~ x + z3, simdata) fit3 <- lm(y ~ x + z3 + z8, simdata) fit4 <- lm(y ~ x + z3 + z7, simdata) #' #' Here are two that might not. Why? #' fit5 <- lm(y ~ x, simdata) fit6 <- lm(y ~ x + z3 + m1, simdata) #' summary(fit1) summary(fit2) summary(fit3) summary(fit4) summary(fit5) summary(fit6) #' \newpage #' **Create a table with estimated values and SEs of $\hat{\beta}_x$** #' library(spida2) waldf(fit1) #' \newpage #' The `lapply` function applies its second argument (a function) to #' each element of the first argument, which should be a **vector** or #' a **list**. Note that a **data.frame** is a list whose elements #' are the variables. fitlist <- list(fit1, fit2, fit3, fit4, fit5, fit6) coeflist <- lapply(fitlist, function(fit) { ret <- waldf(fit)['x', c('coef','se','p-value')] ret$rse <- sqrt(mean(resid(fit)^2)) rownames(ret) <- as.character(formula(fit))[3] ret } ) #' \newpage #' This created a list whose elements are one-row data frames. coeflist #' #' We'd like to make them into rows of one data frame. #' #' The `rbind` function might do this. "rbind" 'binds' its arguments #' into rows of the same data.frame or matrix. #' #' The `cbind` function does the same thing to columns. #' #' We could do this with #' rbind(coeflist[[1]], coeflist[[2]], coeflist[[3]], coeflist[[4]]) # etc. #' but that wouldn't be efficient if there are many arguments. #' Also, you wouldn't be able to use such code within a function #' in which the number of arguments can vary. #' #' The solution is the `do.call` function, which takes two arguments: #' #' 1) a function #' 2) a list whose elements are used as individual arguments to the function #' do.call(rbind, coeflist) #' #' Let's fit a few more models that provide an estimate of #' the causal effect and compare SEs of $\hat{\beta}_x$. #' fit7 <- lm(y ~ x + z1 + z7, simdata) fit8 <- lm(y ~ x + z2, simdata) fit9 <- lm(y ~ x + z1 + z2 + z3, simdata) fit10 <- lm(y ~ x + z1 + z2 + z3 + z8, simdata) fitlist2 <- c(fitlist, list(fit7, fit8, fit9, fit10)) summ_table <- lapply(fitlist2, function(fit) { ret <- waldf(fit)['x', c('coef','se','p-value')] ret$sy <- sqrt(mean(resid(fit)^2)) rownames(ret) <- as.character(formula(fit))[3] ret } ) %>% do.call(rbind,.) summ_table #' Let's keep track of the models that are causal summ_table$causal <- c('yes','no')[c(1,1,1,1,2,2,1,1,1,1)] # when adding only one variable I don't bother using `within` summ_table <- sortdf(summ_table, ~ se) # use idempotent programming # when modifying an object #' `sortdf` is in **spida2**. The official way to order a data.frame in R: #' `summ_table <- summ_table[order(summ_table$se),]` summ_table summ_table$`p-value` <- pfmt(summ_table$`p-value`) summ_table #+ table kable(summ_table, caption = r"(Estimates of $\hat{\beta}_x$)") #+ library(kableExtra) #+ table2 summ_table2 <- cbind(formula = rownames(summ_table), summ_table) rownames(summ_table2) <- NULL kable(summ_table2, caption = r"(Estimates of $\hat{\beta}_x$)", booktabs = TRUE) %>% kable_styling #' #' # Graphic presentation #' gtab <- summ_table2[rep(1:nrow(summ_table2), each = 3),] gtab <- within( gtab, { factor <- c(-1,0,1) # replicates to number of rows in gtab val <- coef + factor * se formula <- as.factor(formula) formula <- reorder(formula, se) pch <- c(40, 16, 41)[factor + 2] col <- dplyr::if_else(causal == 'yes', 'blue','red') } ) gtab library(lattice) library(latticeExtra) library(latex2exp) #' \newpage #+ fig.cap=r"($\hat{\beta}_x \pm SE(\hat{\beta}_x)$ with $\beta_x = 1.6$)",fig.height=3.5,fig.width=4.5 p <- xyplot(formula ~ val, gtab, type = 'l', groups = formula, lwd = 3, ylab = '', col = c("red",'blue')[c(2,2,1,2,2,2,2,2,1,2)], xlab = list(label = TeX(r"(Interval estimates of $\beta_x$)"))) + xyplot(formula ~ val, type = 'p', pch = 16, groups = formula, col = c("red",'blue')[c(2,2,1,2,2,2,2,2,1,2)], data = subset(gtab, factor == 0)) + layer_(panel.grid(h=-1,v=-1)) + layer(panel.abline(v = 1.6, lwd = 2)) p #+ fig.cap=r"($\hat{\beta}_x \pm SE(\hat{\beta}_x)$ with $\beta_x = 1.6$, focusing on 1.6)",fig.height=3.5,fig.width=4.5 update(p,xlim = c(1.6-.05,1.6+.05) ) #' #' Why is there such a large variability in the SEs? #' #' The answer is around the corner, when we see the **Frisch-Waugh-Lovell Theorem**. #' #' # Exercises #' #' 1) What would happen if you were to use a sample size 100 times larger? Try it. #' 2) Think of an example in everyday life with at least one #' backdoor path and at least one causal path in addition to the #' direct path from X to Y. Generate plausible data. Note that #' the means of your variables don't have to be zero, nor the #' standard deviations one. #' Play with different models, classify whether they are causal #' or not, and explore the differences in estimates and standard errors. #' 3) If you are familiar with **ggplot2**, use it to reproduce graphics like those #' in this scripts. #' #' #+ knitr::knit_exit() #+ library(p3d) Init3d(cex = 2) summ_table$causal <- factor(summ_table$causal) Plot3d(se ~ coef + sy, summ_table) summ_table Id3d() as.data.frame(coeflist) coeflist$causal <- rep(c('yes','no'),c(4,2)) library(lattice) library(latticeExtra) xyplot( fitlist2 %>% c(fitlist[1:4], lapply( function(fit) { ret <- summary(fit)['x',c('coef','SE','p-value')] rownames(ret) <- formula(fit)[3] } ) ) %>% do.call(rbind, .) )] } ) knit_exit() #+ setup2 ------ knitr::opts_chunk$set(comment = ' ',fig.height = 3.2,fig.width=5) library(spida2) #+ # # For interactive use, the following sets the # working directory to this file's location. # setwd(this.path::here()) # library(lattice) library(latticeExtra) library(spida2) # remotes::install_github('gmonette/spida2') #' #' This is made-up data designed to be consistent with #' a study reported on CBC radio some years ago. #' #' A physician was interviewed who reported concerns #' over the use of sunscreen lotion because it was found #' that users of sunscreen lotion had a higher incidence #' of skin cancer than non-users. It was suggested that #' sun screen lotion might contain some harmful chemicals. #' #' The following simulated data is intended to help explore #' possible explanations. Do your own research to find out #' whether there is a current consensus on risk factors #' associated with the use of sunscreen lotion. #' \newpage # Direct text input of a data frame: ddin <- read.table(header = TRUE, text = " SkinDamage.mean SunScreen SunExposure n 6 1 1 50 9 0 1 10 3 1 0 10 4 0 0 70 ") dim(ddin) # Creating a big data frame by repeating rows of a samll one: dd <- ddin[rep(1:nrow(ddin), ddin$n), ] # What does this do? dim(dd) # Generate a random random seed sample(10000000,1) set.seed(9607680) # and use the same one for reproducibility dd$eps <- rnorm(nrow(dd)) # generate random epsilons dd$SkinDamage <- with(dd, SkinDamage.mean + ave(eps, SunExposure, SunScreen, FUN = scale)) # What does this do? #' #' The function 'ave' is a very powerful tool that applies the argument #' of FUN to each chunk of `eps` within each combination of values #' of `SunExposure` and `SunScreen`. The function `scale` rescales #' its input (with a location-scale transformation) so the result #' has mean 0 and standard deviation 1. Thus, the #' **sample average** of `SkinDamage` is equal to `SkinDamage.mean`. #' Otherwise, the sample average would differ by a small random amount. #' #' The only purpose in doing this here is to make the results come out #' as nice round numbers for presentation in class. #' \newpage xyplot(SkinDamage ~ SunScreen, dd) #' \newpage xyplot(SkinDamage ~ jitter(SunScreen, .1), dd) #' \newpage xyplot(SkinDamage ~ jitter(SunScreen, .1), dd, groups = SunExposure, auto.key = T) #' \newpage xyplot(SkinDamage ~ jitter(SunScreen, .1), dd, groups = SunExposure, auto.key = list(reverse.rows = T)) # why do this? #' \newpage xyplot(SkinDamage ~ jitter(SunScreen, .1), dd, groups = SunExposure, auto.key = list(reverse.rows = T), par.settings = list(superpose.symbol = list(pch = c(16,17),alpha = .5))) # why do this? #' \newpage p <- xyplot(SkinDamage ~ jitter(SunScreen, .1), dd, groups = SunExposure, auto.key = list(reverse.rows = T, title = 'Sun\nExposure',cex.title = 1), par.settings = list(superpose.symbol = list(pch = c(16,17),alpha = .5)), xlab = 'SunScreen') #' \newpage p #' \newpage #' #' How helpful is SunScreen? #' fit <- lm(SkinDamage ~ SunScreen, dd) summary(fit) dd$fit <- predict(fit) #' \newpage p + xyplot(fit ~ SunScreen, dd, type = 'b', lwd = 3) #' \newpage #' It looks like using SunScreen results in higher levels of Skin Damage #' What if we 'control' for SunExposure #' fit2 <- lm(SkinDamage ~ SunScreen * SunExposure, dd) summary(fit2) # What does this mean?? #' #' Interpreting output: Just think of partial derivatives! #' dd$fit2 <- predict(fit2) spida2::up(subset(dd, select = c(fit2, SunExposure, SunScreen)), ~ SunExposure + SunScreen) xyplot(fit2 ~ SunScreen, dd, type = 'b', groups = SunExposure, ylim = c(2,10), auto.key = list(title = 'Sun\nExposure')) p + xyplot(fit2 ~ SunScreen, dd, type = 'b', groups = SunExposure, lwd = 3) pall <- update( # creates the plot p, legend = NULL, auto.key= list(lines = T, title = 'Sun\nExposure', reverse.rows = TRUE), par.settings = list( plot.line = list(lwd = 3), superpose.line = list(lwd = 3) ) ) + xyplot(fit2 ~ SunScreen, dd, type = 'b', groups = SunExposure) + xyplot(fit ~ SunScreen, dd, type = 'b', col = 'darkgreen') #' \newpage pall # shows the plot #' \newpage #' What does a model look like without interactions? #' fitni <- lm(SkinDamage ~ SunScreen + SunExposure, dd) dd$fitni <- predict(fitni) #+ fig.cap='The red lines show the fitted values using a model without an interaction term' update(p, legend = NULL, auto.key= list(title = 'Sun Exposure')) + xyplot(fit2 ~ SunScreen, dd, type = 'b', groups = SunExposure) + xyplot(fit ~ SunScreen, dd, type = 'b', col = 'darkgreen') + xyplot(fitni ~ SunScreen, dd, groups = SunExposure, type = 'b', col = 'red') #' \newpage #' #' # Estimating other parameters #' #' Two approaches: #' #' - Refit with a different but equivalent formula #' - Use linear transformations and a Wald test #' #' ## Refitting with equivalent formulas #' fit_conditional1 <- lm(SkinDamage ~ factor(SunExposure)/SunScreen - 1, dd) summary(fit_conditional1) fit_conditional2 <- lm(SkinDamage ~ factor(SunScreen)/SunExposure - 1, dd) summary(fit_conditional2) #' #' ## Using linear transformations and Wald tests #' #' Two functions for Wald tests #' #' - car::lht #' - Advantage: can test non-zero hypotheses #' - Disadvantage: Row of hypothesis matrix must be linearly independent #' - spida2::wald (and walddf) #' - Advantage: #' - can test linearly dependent rows #' - spida2::walddf creates a data frame that is easy to plot #' - can simultaneously test subsets of coefficient matched with a regular #' expression #' - Disadvantage: Hypotheses must be = 0 #' #' fit <- lm(SkinDamage ~ SunScreen * SunExposure, dd) summary(fit) L <- rbind( "Effect of X | Z = 0" = c(0, 1, 0, 0), "Effect of X | Z = 1" = c(0, 1, 0, 1), "Effect of Z | X = 0" = c(0, 0, 1, 0), "Effect of X | Z = 0" = c(0, 0, 1, 1), "E(Y | X = 1, Z = 1)" = c(1, 1, 1, 1) # you can keep going ) library(spida2) wald(fit, L) # overall test is an F-test for all hypotheses = 0 fit wald(fit, 'SunScreen') # test whether we can drop SunScreen altogether wald(fit, 'SunExposure') # test whether we can drop SunExposure altogether wald(fit, ':') # test whether we can drop all interaction terms (useful if many) #' #' #' # The Paik-Agresti Diagram #' #' Two diagrams that summarize all #' #' - The Paik-Agresti Diagram #' - Small project: Add the Liu-Meng marginal effect of Z line #' UCBAdmissions UCBAdmissions %>% as.data.frame UCBAdmissions %>% as.data.frame -> ducb fit <- glm(Admit ~ Gender, ducb, family = binomial, weights = Freq) summary(fit) fitbydept <- glm(Admit ~ Dept/Gender - 1, ducb, family = binomial, weights = Freq) summary(fitbydept) library(spida2) #' \newpage paik(Admit ~ Gender + Dept | Freq, data = ducb, ylab = 'Proportion Rejected') #' \newpage #' Adjusting limits for the y-axis paik(Admit ~ Gender + Dept | Freq, data = ducb, ylab = 'Proportion Rejected', ylim = c(0,1)) #' #' # Exercise 1: How a superficial model can hide inequity. #' #' Consider the `death.penalty` data in the **spida2** package. It details the #' frequencies with which defendants received the death penalty (not necessarily #' carried out) in 674 homicide trials in the state of Florida from #' 1976 to 1987. The cases are classified according the the defendant's and #' the victim's race. #' #' Using only the defendant's race as a predictor could suggest that the #' the proportion of defendants receiving death sentences is relative similar #' between the two races. #' #' Adjusting for the victims' race, however, suggests a very different and #' disturbing picture. This illustrates how a superficial analysis can fail to #' reveal important issues and lead to entirely different conclusions. #' #' # Exercise 2: Simpson's Paradox and Choosing a Restaurant: A hypothetical example (tricky!) #' #' Mary (a woman) is choosing between restaurants A and B to take her friend, #' John (a man), out for dinner. #' Restaurant A has an average rating of 4.1 and restaurant B of 4.3.\newline #' \newline #' But looking at ratings by gender, among men, restaurant A has a rating of 4.0 and #' restaurant B of 3.8. Among women, restaurant A has a rating of 4.6 and #' restaurant B of 4.4. It seems that men and women separately prefer restaurant #' A but together they prefer restaurant B! #' #' 1) Draw a Paik-Agresti or other suitable diagram conditioning on gender to represent this data. #' 2) Draw a causal graph describing a plausible relationship among the three variables: restaurant, #' gender and ratings. #' 3) Without reference to possible causal mechanisms, explain numerically #' how the discrepancy between relative overall ratings #' and relative ratings by gender occurs. #' 4) Assuming that there are no other relevant factors related to restaurant ratings, #' what kind of variable would gender be in this context (you may now consider possible causal mechanisms)? #' Which restaurant should Mary choose? Why? #' 5) Now suppose that you have the same data but the variables are different. #' Restaurant A has an average rating of 4.1 and restaurant B of 4.3, as before. #' Each restaurant has two waiters #' who, very coincidentally, happen to have the same two names in #' each restaurant: Mr. Good and Mr. Bad. #' Among patrons served Mr. Bad in restaurant A, the ratings #' are 4.0 and those served by Mr. Bad (different waiter but the #' same name) in restaurant B, the ratings are 3.8. #' Among patrons served Mr. Good in restaurant A, the ratings #' are 4.6 and those served by Mr. Good in restaurant B, the ratings are 4.4. #' Without reference to possible causal mechanisms, explain numerically #' how such a discrepancy arises. #' 6) Assuming that there are no other relevant factors related to restaurant ratings and #' that the propensity for diners to report a rating is independent of waiters, #' what kind of variable would 'waiter' be in this context (you may now consider possible causal mechanisms)? #' Which restaurant should Mary choose? Why? #' #' #' \newpage #' #' # The Liu-Meng Diagram: The main effect of Z #' #' You could always flip the roles of X and Z and draw a Paik-Agresti #' diagram for the main effect of Z. #' #' But that wouldn't show the two in the same diagram. #' #' The Liu-Meng diagram (personal communication over lunch with #' Xiao-Li Meng, Liu was his undergraduate student who worked on #' the idea with Xiao-Li. #' Unfortunately I don't know Liu's given name but someday I #' must find out. #' #' 1) Plot the points: #' - (mean of X, mean of Y)| Z = 0 #' - this is a weighted mean of the #' - (mean of X, mean of Y)| Z = 0, X = 0, and #' - (mean of X, mean of Y)| Z = 0, X = 1 #' Thus, it lies in the convex hull of these two points, which is the line joining the two points! #' - (mean of X, mean of Y)| Z = 1 #' - this is a weighted mean of the #' - (mean of X, mean of Y)| Z = 1, X = 0, and #' - (mean of X, mean of Y)| Z = 1, X = 1 #' Similarly, it also on the line joining these two points. #' - The vertical distance between these two points is the #' main effect of Z on Y. #' - The horizontal distance between these two points is the #' main effect of Z on X. #' #' Small project: Write a function that produces a Liu-Meng #' diagram. #' #' R makes complicated things easy .... and sometimes easy things complicated: #' ddZ <- spida2::up(dd, ~ SunExposure, agg = ~ SkinDamage + SunScreen) ddZ #' \newpage pall + xyplot(SkinDamage ~ SunScreen, ddZ, type = 'b', lwd = 2, pch = 16, cex = 2)