lmer_partial_pooling.md

library(tidyverse)
library(broom)
library(lme4)

set.seed(1234)

# Make fake data
# 2x2x2 factorial design
df <- tibble(id = 1:1000, 
             outcome = rnorm(1000),
             main_treatment = sample(c("A", "B"), 1000, replace = TRUE),
             sub_treatment1 = sample(c("J", "K"), 1000, replace = TRUE),
             sub_treatment2 = sample(c("X", "Y"), 1000, replace = TRUE))

# Group means
model_simpler <- lm(outcome ~ 0 + main_treatment, data = df)
tidy(model_simpler)
#> # A tibble: 2 x 5
#>   term            estimate std.error statistic p.value
#>   <chr>              <dbl>     <dbl>     <dbl>   <dbl>
#> 1 main_treatmentA  -0.0244    0.0453    -0.539   0.590
#> 2 main_treatmentB  -0.0286    0.0440    -0.651   0.515

# Complicated group means - combine all these interaction terms with different
# intercepts to get nested group means
model_3_ways <- lm(outcome ~ main_treatment:sub_treatment1:sub_treatment2,
                   data = df)
tidy(model_3_ways)
#> # A tibble: 8 x 5
#>   term                                 estimate std.error statistic p.value
#>   <chr>                                   <dbl>     <dbl>     <dbl>   <dbl>
#> 1 (Intercept)                           0.00864    0.0952    0.0907   0.928
#> 2 main_treatmentA:sub_treatment1J:sub… -0.0702     0.134    -0.525    0.600
#> 3 main_treatmentB:sub_treatment1J:sub…  0.00884    0.128     0.0693   0.945
#> 4 main_treatmentA:sub_treatment1K:sub…  0.0670     0.135     0.496    0.620
#> 5 main_treatmentB:sub_treatment1K:sub… -0.0602     0.129    -0.467    0.641
#> 6 main_treatmentA:sub_treatment1J:sub… -0.130      0.129    -1.01     0.315
#> 7 main_treatmentB:sub_treatment1J:sub… -0.0924     0.128    -0.720    0.472
#> 8 main_treatmentA:sub_treatment1K:sub…  0.0139     0.129     0.108    0.914

# Avoid multiple comparisons issue by partially pooling
# See p. 9 in http://www.stat.columbia.edu/~gelman/research/unpublished/multiple2.pdf
model_fancy_simple <- lmer(outcome ~ 0 + main_treatment + (1 + main_treatment | sub_treatment1),
                           data = df)
#> singular fit
summary(model_fancy_simple)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: 
#> outcome ~ 0 + main_treatment + (1 + main_treatment | sub_treatment1)
#>    Data: df
#> 
#> REML criterion at convergence: 2839.8
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -3.3797 -0.6611 -0.0130  0.6276  3.1878 
#> 
#> Random effects:
#>  Groups         Name            Variance Std.Dev. Corr 
#>  sub_treatment1 (Intercept)     0.005771 0.07597       
#>                 main_treatmentB 0.005104 0.07144  -1.00
#>  Residual                       0.994268 0.99713       
#> Number of obs: 1000, groups:  sub_treatment1, 2
#> 
#> Fixed effects:
#>                 Estimate Std. Error t value
#> main_treatmentA -0.02385    0.07026  -0.339
#> main_treatmentB -0.02848    0.04406  -0.647
#> 
#> Correlation of Fixed Effects:
#>             mn_trA
#> man_trtmntB 0.056 
#> convergence code: 0
#> singular fit

^{Created on 2019-02-26 by the reprex package (v0.2.1)}