Demonstrate Lord's Paradox example, following Michael Clarke

clark-lord-paradox.R

# Demonstrate Lord's Paradox example, following Michael Clarke
# https://m-clark.github.io/docs/lord/index.html

library(ggplot2)
library(dplyr)

#' ## Data example: Initial and final weight for boys & girls
#' 
#' Make a plot showing the distribution of gain scores perpendicular to the line of slope = 1
#' 

#' Generate the data
set.seed(1234)
N <- 200
group   <- rep(c(0, 1), each = N/2)
initial <- .75*group + rnorm(N, sd=.25)
final   <- .4*initial + .5*group + rnorm(N, sd=.1)
change  <- final - initial
df <- data.frame(id = factor(1:N), 
                group = factor(group, 
                               labels = c('Female', 'Male')), 
                initial, 
                final, 
                change)

#' ## plot, with regression lines and data ellipses for each group
#' Add the line of unit slope, which is orthogonal to an axis for the gain score
ggplot(df, aes(x = initial, y = final, color = group)) +
	geom_point() + 
	geom_smooth(method = "lm", formula  =  y~x) +
  stat_ellipse(size = 1.2) +
  geom_abline(slope  =  1, color = "black", size = 1.2) +
  coord_fixed(xlim = c(-.6, 1.2), 
              ylim = c(-.6, 1.2)) +
  theme_bw() +
  theme(legend.position = c(.15, .85)) 

# code from: https://stackoverflow.com/questions/71001954/ggplot2-projecting-points-or-distribution-on-a-non-orthogonal-eg-45-degree

#' ## Create a distribution plot of the gain scores.
(my_hist <- df %>%
    mutate(gain = final - initial) %>%
    ggplot(aes(gain)) +
    geom_density())


#' Now we can extract the guts of that plot, and transform the coordinates to where 
#' we want them to appear in the combined plot:
  
a <- ggplot_build(my_hist)
rot = pi * 3/4
diag_hist <- tibble(
  x = a[["data"]][[1]][["x"]],
  y = a[["data"]][[1]][["y"]]
  ) %>%
  # squish
  mutate(y = y*0.2) %>%
  # rotate 135 deg CCW
  mutate(xy =   x*cos(rot) - y*sin(rot),
         dens = x*sin(rot) + y*cos(rot)) %>%
  # slide toward new origin
  mutate(xy = xy - 0.7,  #  origin shift based on plot range below
         dens = dens - 0.7)

#' ## Replot, with annotation for gain score distribution
#' 
ggplot(df, aes(x = initial, y = final, color = group)) +
  geom_point() + 
  geom_smooth(method = "lm", formula  =  y~x) +
  stat_ellipse(size = 1.2) +
  geom_abline(slope  =  1, color = "black", size = 1.2) +
  coord_fixed(clip = "off", 
              xlim = c(-0.7,1.6),
              ylim = c(-0.7,1.6), 
              expand = expansion(0)) +
  annotate("segment", x = -1.4, xend = 0, y = 0, yend = -1.4) +
  annotate("path", x = diag_hist$xy, y = diag_hist$dens) +
  theme_bw() +
  theme(legend.position = c(.15, .85), 
        plot.margin = unit(c(.1,.1,2,2), "cm"))