#' ---
#' title: |
#'   Simpson's Paradox:  
#'   The Trapezoid of Means    
#'   with Paik-Agresti and Liu-Meng Diagrams
#' 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
#' ---
#' \large
#' 
#' # Setup
#' 
#' This R script is designed to be run either interactively
#' one line at a time, or as a Rmarkdown script that creates
#' a pdf file.
#' 
#' In a clean version of a report, this portion of the 
#' Rmarkdown script would be hidden by using a
#' 'include=FALSE' option for the 'chunk header'. That is, the
#' header line for the chunk would read:
#' 
#' > #+ setup, include=FALSE
#' 
#+ setup ------
knitr::opts_chunk$set(comment = '   ',fig.height = 3.2,fig.width=5) 
#+
# 
# 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)