Simulating the data

In the first part of this lab we are going to replicate some of the ideas presented in a post written by Michael Clark.

We are going to use the faux package to simulate nested data (observations nested within clusters). First, we define the number of clusters n_cluster and the number of observations in each cluster, obs_per_cluster. For example, we might want to have 8 clusters with 100 observations in each cluster (to keep things simple, let’s stick to a “balanced” data-set).

obs_per_cluster <-  100
n_cluster <-  8

Let’s generate some data, where (level 1) observed outcomes are related to predictors in the following manner: Y=ai+bi⋅X+ϵwhere: ϵ∼(0,σ2) But rather than being regular population parameters, ai and bi are random intercepts and slopes, such that each cluster has its own, unique intercept and slope. The random intercept is a level 2 random effect: ai=α0+uiwhere: ui∼(0,σ2u) And the random slope is a level 2 random effect: bi=β0+viwhere: vi∼(0,σ2v)

In this stage we have only one predictor at level one X. We have no predictor at level 2. We also have to define six different population parameters. We choose α0=1,β0=0.5 and three different standard deviations: σ=1,σu=0.5,σv=0.25,ρu,v=0.5

# Fixed intercept
alpha_0 <-1 # choose a number you like, e.g. 1  

# Fixed slope
beta_0 <- 0.5 #  e.g. 0.5

# Residuals st.dev (level 1) 
sig  <- 1 #  e.g. 1

# Random intercept st.dev (level 2) 
sig_u <- 0.5 #  e.g. 0.5

# Random slope st.dev (level 2)
sig_v <- 0.25 #  e.g. 0.25

# Independent random effects
cor <- 0.5 #  e.g. 0.5

Simulate our data

# Create a dataset that can be reproduced
set.seed(888)  
# step 1: create the units 
df <- add_random(
    cluster = n_cluster,  # first the clusters
    obs = obs_per_cluster # second, the observations
    ) 

# Print the result 
# try to figure out what just happened
df |> 
  print(n = 15)

## # A tibble: 800 × 2
##    cluster  obs   
##    <chr>    <chr> 
##  1 cluster1 obs001
##  2 cluster1 obs002
##  3 cluster1 obs003
##  4 cluster1 obs004
##  5 cluster1 obs005
##  6 cluster1 obs006
##  7 cluster1 obs007
##  8 cluster1 obs008
##  9 cluster1 obs009
## 10 cluster1 obs010
## 11 cluster1 obs011
## 12 cluster1 obs012
## 13 cluster1 obs013
## 14 cluster1 obs014
## 15 cluster1 obs015
## # ℹ 785 more rows

# step 2: Add anything at level 2
# Here we have random intercept/slope 
df <- df |>
  add_ranef(
    "cluster",
    u_i = sig_u,
    v_i = sig_v,
    .cors = cor
    ) |> mutate(
      a_i = alpha_0 + u_i, 
      b_i = beta_0  + v_i
    )

# Print the result 
# try to figure out what just happened
df |> 
  print(n = 15)

## # A tibble: 800 × 6
##    cluster  obs      u_i   v_i   a_i   b_i
##    <chr>    <chr>  <dbl> <dbl> <dbl> <dbl>
##  1 cluster1 obs001  1.01 0.167  2.01 0.667
##  2 cluster1 obs002  1.01 0.167  2.01 0.667
##  3 cluster1 obs003  1.01 0.167  2.01 0.667
##  4 cluster1 obs004  1.01 0.167  2.01 0.667
##  5 cluster1 obs005  1.01 0.167  2.01 0.667
##  6 cluster1 obs006  1.01 0.167  2.01 0.667
##  7 cluster1 obs007  1.01 0.167  2.01 0.667
##  8 cluster1 obs008  1.01 0.167  2.01 0.667
##  9 cluster1 obs009  1.01 0.167  2.01 0.667
## 10 cluster1 obs010  1.01 0.167  2.01 0.667
## 11 cluster1 obs011  1.01 0.167  2.01 0.667
## 12 cluster1 obs012  1.01 0.167  2.01 0.667
## 13 cluster1 obs013  1.01 0.167  2.01 0.667
## 14 cluster1 obs014  1.01 0.167  2.01 0.667
## 15 cluster1 obs015  1.01 0.167  2.01 0.667
## # ℹ 785 more rows

# step 3: add anything at level 1
df <- df |>
  mutate(
    x       = rnorm(n_cluster * obs_per_cluster),
    epsilon = rnorm(n_cluster * obs_per_cluster, 0, sig),
    y = a_i + b_i * x + epsilon
    )


# Print the result 
# try to figure out what just happened
df |> 
  print(n = 20)

## # A tibble: 800 × 9
##    cluster  obs      u_i   v_i   a_i   b_i       x epsilon       y
##    <chr>    <chr>  <dbl> <dbl> <dbl> <dbl>   <dbl>   <dbl>   <dbl>
##  1 cluster1 obs001  1.01 0.167  2.01 0.667  1.18    1.76    4.56  
##  2 cluster1 obs002  1.01 0.167  2.01 0.667 -0.191   0.0103  1.89  
##  3 cluster1 obs003  1.01 0.167  2.01 0.667 -2.19    0.0497  0.593 
##  4 cluster1 obs004  1.01 0.167  2.01 0.667 -0.338   0.663   2.44  
##  5 cluster1 obs005  1.01 0.167  2.01 0.667 -0.0519 -1.88    0.0916
##  6 cluster1 obs006  1.01 0.167  2.01 0.667 -0.280   0.819   2.64  
##  7 cluster1 obs007  1.01 0.167  2.01 0.667 -0.749  -0.220   1.29  
##  8 cluster1 obs008  1.01 0.167  2.01 0.667  0.777  -2.67   -0.148 
##  9 cluster1 obs009  1.01 0.167  2.01 0.667 -0.223   1.55    3.41  
## 10 cluster1 obs010  1.01 0.167  2.01 0.667 -0.142  -0.898   1.01  
## 11 cluster1 obs011  1.01 0.167  2.01 0.667 -0.173  -1.66    0.233 
## 12 cluster1 obs012  1.01 0.167  2.01 0.667  2.97   -1.58    2.41  
## 13 cluster1 obs013  1.01 0.167  2.01 0.667  0.235   0.487   2.65  
## 14 cluster1 obs014  1.01 0.167  2.01 0.667  0.328  -0.598   1.63  
## 15 cluster1 obs015  1.01 0.167  2.01 0.667  1.64   -2.26    0.839 
## 16 cluster1 obs016  1.01 0.167  2.01 0.667 -1.00   -2.12   -0.777 
## 17 cluster1 obs017  1.01 0.167  2.01 0.667  0.196   0.102   2.24  
## 18 cluster1 obs018  1.01 0.167  2.01 0.667  1.12    1.29    4.04  
## 19 cluster1 obs019  1.01 0.167  2.01 0.667 -0.0120  0.524   2.52  
## 20 cluster1 obs020  1.01 0.167  2.01 0.667 -1.60   -0.0359  0.906 
## # ℹ 780 more rows

# For the models, we really only need 
# the cluster id, the y and the x

df  <- df |> select(cluster, x, y)

# Print the result 
# try to figure out what just happened
df |> 
  print(n = 20)

## # A tibble: 800 × 3
##    cluster        x       y
##    <chr>      <dbl>   <dbl>
##  1 cluster1  1.18    4.56  
##  2 cluster1 -0.191   1.89  
##  3 cluster1 -2.19    0.593 
##  4 cluster1 -0.338   2.44  
##  5 cluster1 -0.0519  0.0916
##  6 cluster1 -0.280   2.64  
##  7 cluster1 -0.749   1.29  
##  8 cluster1  0.777  -0.148 
##  9 cluster1 -0.223   3.41  
## 10 cluster1 -0.142   1.01  
## 11 cluster1 -0.173   0.233 
## 12 cluster1  2.97    2.41  
## 13 cluster1  0.235   2.65  
## 14 cluster1  0.328   1.63  
## 15 cluster1  1.64    0.839 
## 16 cluster1 -1.00   -0.777 
## 17 cluster1  0.196   2.24  
## 18 cluster1  1.12    4.04  
## 19 cluster1 -0.0120  2.52  
## 20 cluster1 -1.60    0.906 
## # ℹ 780 more rows

Exercise 1

Adjust the parameters, simulate your clustered data and present the results of the three types of models using the modelsummary package. Interpret your models and compare the estimations with each other, and with the six parameters used to simulate your data. Is there one single model that is superior to all the rest in terms of estimating the parameters? If so which is it, and what makes it better? If not, Which model works best for what purpose? In your answer, please pay attention to the standard errors and the Intraclass Correlation Coefficient (ICC) of your multilevel model(s). What does the ICC mean? How much does it deviate from what you would expect? For more about the ICC, please see below.

1.1 Creating the models

# Complete pooling: Using the formula y ~ x
model1 <- lm(y ~ x, df)
summary(model1)

## 
## Call:
## lm(formula = y ~ x, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3620 -0.7753  0.0580  0.7489  3.0002 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.44286    0.04180   34.52   <2e-16 ***
## x            0.57455    0.04413   13.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.182 on 798 degrees of freedom
## Multiple R-squared:  0.1752, Adjusted R-squared:  0.1742 
## F-statistic: 169.5 on 1 and 798 DF,  p-value: < 2.2e-16

# No pooling
# Without interaction: y ~ x + cluster
model2 <- lm(y ~ x + cluster, df)
summary(model2)

## 
## Call:
## lm(formula = y ~ x + cluster, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3430 -0.7171  0.0069  0.7030  3.2091 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.08949    0.10373  20.144  < 2e-16 ***
## x                0.53890    0.03897  13.827  < 2e-16 ***
## clustercluster2 -0.21004    0.14709  -1.428    0.154    
## clustercluster3 -1.46952    0.14671 -10.016  < 2e-16 ***
## clustercluster4 -0.84613    0.14670  -5.768 1.15e-08 ***
## clustercluster5 -0.19608    0.14669  -1.337    0.182    
## clustercluster6 -0.91967    0.14673  -6.268 6.02e-10 ***
## clustercluster7 -0.05306    0.14670  -0.362    0.718    
## clustercluster8 -1.48028    0.14671 -10.090  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.037 on 791 degrees of freedom
## Multiple R-squared:  0.3706, Adjusted R-squared:  0.3642 
## F-statistic: 58.21 on 8 and 791 DF,  p-value: < 2.2e-16

# With interaction: y ~ x + cluster + x:cluster
model3 <- lm(y ~ x + cluster + x:cluster, df)
summary(model3)

## 
## Call:
## lm(formula = y ~ x + cluster + x:cluster, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1839 -0.6748 -0.0029  0.7242  3.1635 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.088158   0.102495  20.373  < 2e-16 ***
## x                  0.483323   0.098917   4.886 1.25e-06 ***
## clustercluster2   -0.164792   0.147537  -1.117  0.26436    
## clustercluster3   -1.473608   0.145238 -10.146  < 2e-16 ***
## clustercluster4   -0.837578   0.144951  -5.778 1.09e-08 ***
## clustercluster5   -0.189039   0.145010  -1.304  0.19274    
## clustercluster6   -0.917547   0.145798  -6.293 5.17e-10 ***
## clustercluster7   -0.061670   0.144960  -0.425  0.67064    
## clustercluster8   -1.469061   0.145190 -10.118  < 2e-16 ***
## x:clustercluster2 -0.118099   0.147322  -0.802  0.42301    
## x:clustercluster3 -0.006334   0.146260  -0.043  0.96547    
## x:clustercluster4 -0.262204   0.145179  -1.806  0.07129 .  
## x:clustercluster5  0.179757   0.143592   1.252  0.21099    
## x:clustercluster6  0.062611   0.171976   0.364  0.71590    
## x:clustercluster7  0.429776   0.147513   2.913  0.00368 ** 
## x:clustercluster8  0.177528   0.145676   1.219  0.22335    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.025 on 784 degrees of freedom
## Multiple R-squared:  0.3912, Adjusted R-squared:  0.3795 
## F-statistic: 33.58 on 15 and 784 DF,  p-value: < 2.2e-16

# Multilevel model (partial pooling)
# Using random intercepts
model4 <- lmer(y ~ x + (1 | cluster), df)
summary(model4)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: y ~ x + (1 | cluster)
##    Data: df
## 
## REML criterion at convergence: 2361.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2093 -0.6995  0.0045  0.6744  3.0720 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  cluster  (Intercept) 0.3678   0.6065  
##  Residual             1.0759   1.0373  
## Number of obs: 800, groups:  cluster, 8
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   1.44265    0.21754   6.99863   6.632 0.000296 ***
## x             0.53993    0.03897 791.58664  13.856  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##   (Intr)
## x 0.001

# Using random intercepts and slopes
model5 <- lmer(y ~ x + (x | cluster), df)
summary(model5)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: y ~ x + (x | cluster)
##    Data: df
## 
## REML criterion at convergence: 2350.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.12721 -0.66833 -0.01454  0.70420  3.06881 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  cluster  (Intercept) 0.37050  0.6087       
##           x           0.03285  0.1812   0.21
##  Residual             1.04987  1.0246       
## Number of obs: 800, groups:  cluster, 8
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  1.44662    0.21826 7.00104   6.628 0.000296 ***
## x            0.54128    0.07501 7.10418   7.216 0.000163 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##   (Intr)
## x 0.182

1.2 Comparing the models

modelsummary(list("Complete pooling" = model1, 
                  "No pooling" = model2, 
                  "No pooling with interaction" = model3, 
                  "Multilevel model with random intercepts" = model4, 
                  "Multilevel model with random intercepts and slopes" = model5),
             stars = TRUE)

	Complete pooling	No pooling	No pooling with interaction	Multilevel model with random intercepts	Multilevel model with random intercepts and slopes
(Intercept)	1.443***	2.089***	2.088***	1.443***	1.447***
	(0.042)	(0.104)	(0.102)	(0.218)	(0.218)
x	0.575***	0.539***	0.483***	0.540***	0.541***
	(0.044)	(0.039)	(0.099)	(0.039)	(0.075)
clustercluster2		−0.210	−0.165
		(0.147)	(0.148)
clustercluster3		−1.470***	−1.474***
		(0.147)	(0.145)
clustercluster4		−0.846***	−0.838***
		(0.147)	(0.145)
clustercluster5		−0.196	−0.189
		(0.147)	(0.145)
clustercluster6		−0.920***	−0.918***
		(0.147)	(0.146)
clustercluster7		−0.053	−0.062
		(0.147)	(0.145)
clustercluster8		−1.480***	−1.469***
		(0.147)	(0.145)
x × clustercluster2			−0.118
			(0.147)
x × clustercluster3			−0.006
			(0.146)
x × clustercluster4			−0.262+
			(0.145)
x × clustercluster5			0.180
			(0.144)
x × clustercluster6			0.063
			(0.172)
x × clustercluster7			0.430**
			(0.148)
x × clustercluster8			0.178
			(0.146)
SD (Intercept cluster)				0.606	0.609
SD (x cluster)					0.181
Cor (Intercept~x cluster)					0.214
SD (Observations)				1.037	1.025
Num.Obs.	800	800	800	800	800
R2	0.175	0.371	0.391
R2 Adj.	0.174	0.364	0.380
R2 Marg.				0.154	0.154
R2 Cond.				0.369	0.387
AIC	2542.0	2339.8	2327.1	2369.2	2362.8
BIC	2556.1	2386.6	2406.8	2387.9	2390.9
ICC				0.3	0.3
Log.Lik.	−1268.015	−1159.889	−1146.574
F	169.541	58.215	33.584
RMSE	1.18	1.03	1.01	1.03	1.02
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Answer

Using this simulated data, we see that the complete pooling and multilevel models can estimate the intercept and slope fairly well. The complete pooling model being a good estimation in this case might be explained by little variation between clusters in the simulated data and high number of observations. For the no pooling models, they can estimate the global slope very well but the intercept is quite far from the set parameter in the simulation. However, the errors are not so big in this case, since we have a high number of observations and not so many clusters, not so many degrees of freedom are being consumed in the calculations. Finally, the multilevel models perform very well when estimating the parameters and the estimated values are pretty close from the set values. Nevertheless, the multilevel model adding random slopes doesn’t seem to improve a lot when compared to the multilevel model with only random intercepts. The final decision of which model to use should take into account the nature of the data, if a correlation exists between the level 1 and 2 effects (rho), then the model with random slopes and intercepts should be preferred. From the ICC value (0.3 in both multilevel models), we can see that this dataset doesn’t have a large variability between groups and therefore, the multilevel structure is not that critical since most of the variability is within groups. This might also explain why the complete pooling gave good results in the case as well.

Exercise 2

Use your models to create a caterpillar plot for the intercepts and lattice and/or spaghetti plots to show the features of your multilevel models. Describe and discuss the plots: what new light do they shed? What are the advantages or disadvantages of presenting the models using these plots?

2.1 Caterpillar plots

# First, we need the list of clusters
clusters <- df |> 
  pull(cluster) |> 
  unique()

# predictions of y for each cluster at x = 0: multilevel random intercept
predict_y_model4 <- model4 |> 
  ranef() |> 
  # Used the broom.mixed::augment function 
  # to calculate the confidence intervals
  augment() |> 
  # convert to a tibble
  as_tibble() |> 
  # Include only the random intercept
  filter(variable == "(Intercept)") |> 
  # Now we order the results by ascending estimates
  arrange(estimate) |> 
  mutate(cluster = fct_inorder(level)) 

# predictions of y for each cluster at x = 0: multilevel random intercept and slope
predict_y_model5 <- model5 |> 
  ranef() |> 
  # Used the broom.mixed::augment function 
  # to calculate the confidence intervals
  augment() |> 
  # convert to a tibble
  as_tibble() |> 
  # Include only the random intercept
  filter(variable == "(Intercept)") |> 
  # Now we order the results by ascending estimates
  arrange(estimate) |> 
  mutate(cluster = fct_inorder(level)) 

# Plotting the caterpillars
p4 <- predict_y_model4 |> 
  ggplot(aes(xmin = lb, xmax = ub, y = cluster, x = estimate)) + 
    geom_pointrange( ) + 
    labs(title = "Multilevel Random Intercept", x = "Estimate of the intercept (Y_hat when X=0)", y = NULL) +
    theme_classic()

p5 <- predict_y_model5 |> 
  ggplot(aes(xmin = lb, xmax = ub, y = cluster, x = estimate)) + 
    geom_pointrange( ) + 
    labs(title = "Multilevel Random Intercept and Slope", x = "Estimate of the intercept (Y_hat when X=0)", y = NULL) +
    theme_classic()

library(patchwork)
print(p4 + p5)

2.2 Spaghetti plots

# Spaghetti plots
p4_spa <- predictions(model4) |>  
  ggplot(aes(x, y, color = cluster)) + 
  geom_point() + 
  geom_line(aes(y = estimate)) +
  theme_classic() +
  labs(title = "Multilevel Random Intercept")

p5_spa <- predictions(model5) |>  
  ggplot(aes(x, y, color = cluster)) + 
  geom_point() + 
  geom_line(aes(y = estimate)) +
  theme_classic() +
  labs(title = "Multilevel Random Intercept and Slope")

print(p4_spa + p5_spa + plot_layout(guides = 'collect'))

2.3 Lattice plots

predictions(model4) |>  
  ggplot(aes(x, y, color = cluster)) + 
  geom_point() + 
  geom_line(aes(y = estimate)) +
  theme_classic() +
  labs(title = "Multilevel Random Intercept") +
  facet_wrap(~cluster)

predictions(model5) |>  
  ggplot(aes(x, y, color = cluster)) + 
  geom_point() + 
  geom_line(aes(y = estimate)) +
  theme_classic() +
  labs(title = "Multilevel Random Intercept and Slope") +
  facet_wrap(~cluster)

Answer

These plots help to verify how much variation exists in each cluster and between clusters. From all the plots, we can see that in this simulated dataset there is not much inter or intravariation between groups, which we had already concluded in Exercise 1 by analyzing the models. Nevertheless, these plots are very helpful to have a visual representation of that and helping us choosing the best model. It is also an extra verification that the model adding random slopes doesn’t seem to be critical in this case (this can easily be seen in the lattice plots, where the cluster all have a very similar slope).

Exercise 3

3.1 Load data

loading the radon data and select measurements from the state of Minnesota

radon_mn <- read_rds("https://bit.ly/BMLR-Radon") |> 
  filter(state == "MN") |>
  # drop counties without measurements
  mutate(county = fct_drop(county))

radon_mn %>% 
  mutate(first_floor_text = case_when(
    first_floor == 0 ~ "Basement",
    TRUE ~ "First floor"
  )) -> radon_mn

# Defining pretty names for the labels
label(radon_mn$radon) <- "Randon"
label(radon_mn$log.radon) <- "Log of randon"
label(radon_mn$u) <- "Uranium"

3.2 Explore data

Create a table1 and use GGally::ggpairs to explore your data

# table
table1(~ radon + log.radon + u | first_floor_text,
       data = radon_mn)

	Basement (N=766)	First floor (N=153)	Overall (N=919)
Randon
Mean (SD)	5.07 (4.54)	3.28 (3.83)	4.77 (4.48)
Median [Min, Max]	3.90 [0.100, 48.2]	2.10 [0, 25.4]	3.60 [0, 48.2]
Log of randon
Mean (SD)	1.33 (0.783)	0.713 (0.999)	1.22 (0.853)
Median [Min, Max]	1.36 [-2.30, 3.88]	0.742 [-2.30, 3.23]	1.28 [-2.30, 3.88]
Uranium
Mean (SD)	-0.139 (0.366)	-0.0970 (0.360)	-0.132 (0.365)
Median [Min, Max]	-0.0965 [-0.882, 0.528]	-0.0900 [-0.882, 0.504]	-0.0965 [-0.882, 0.528]

# ggpairs
p1 <- ggpairs(radon_mn,
        columns = c("radon", "log.radon", "u"),
        aes(color = first_floor_text),
        title = 'Radon in the state of Minnesota',
        columnLabels = c(label(radon_mn$radon), label(radon_mn$log.radon),
                         label(radon_mn$u)),
        legend = c(3,3),
        lower = list(continuous = "density", combo = "box_no_facet")
        ) + scale_fill_discrete(name = "")

p1

Excercise 4

You will now visualize the estimated and confidence interval for the average log Radon across different counties, against the (jittered) number of measurements taken in said county. We will use the estimation to compare between two models: the no pooling model and the multilevel model with partial pooling, both without predictors. The horizontal line in each plot represents an estimate of the average radon level across all counties.

Replicate the figures below and explain why the counties with fewer measurements have more variable estimates and larger higher standard errors. What is the main advantage and disadvantage of no-pooling analysis, as illustrated in the plot on the left compared with the multilevel model? Explain.

# count
radon_mn_count <- radon_mn |>
  count(county, name = "n.obs")

# no pooling model
m.no_pool <- radon_mn |>
  lm(log.radon ~ 0 + county, data = _)

# multilevel model
m.partial_pool <- radon_mn |>
  lmer(log.radon ~ (1|county), data = _)

# the fixed effects (no pooling) model
m.no_pool |> 
  predictions(newdata = datagrid(county = unique)) |> 
  as_tibble() |> inner_join(radon_mn_count, by = "county") |> 
  select(county, estimate, conf.low, conf.high, n.obs) |> 
  mutate(mean_estimate = mean(estimate)) %>% 
  ggplot(
    aes(
      x    = n.obs, 
      y    = estimate, 
      ymin = conf.low, 
      ymax = conf.high
      )
    ) + 
  labs(title = "Fixed effects model", 
       y = "Estimated log randon in the county",
       x = "Number of observations") + 
  geom_pointrange(
    position = position_jitter(width = .5, height = 0)
    ) +
  geom_hline(yintercept = 1.34) + 
  scale_x_continuous(trans = "log10") + 
  theme_bw() -> p_nopool

# the multilevel (partial pooling) model
m.partial_pool |> 
  predictions(newdata = datagrid(county = unique)) |> 
  as_tibble() |> inner_join(radon_mn_count, by = "county") |> 
  select(county, estimate, conf.low, conf.high, n.obs) |> 
  ggplot(
    aes(
      x    = n.obs, 
      y    = estimate, 
      ymin = conf.low, 
      ymax = conf.high
      )
    ) + 
  labs(title = "Multilevel model", 
       y = "Estimated log randon in the county",
       x = "Number of observations") + 
  geom_pointrange(
    position = position_jitter(width = .5, height = 0)
    ) +
  geom_hline(yintercept = 1.313) + 
  scale_x_continuous(trans = "log10") + 
  theme_bw() -> p_partial

p_nopool / p_partial

Answer

The main advantage of the no-pooling analysis is that it is very simple to implement and understand. However, it has the disadvantages: it is not very robust to outlines and it is not very efficient. We can see the confidence intervals are much wider in the no-pooling model than in the multilevel model. This is because the no-pooling model is not very efficient, as it is not taking advantage of the information from the other counties and it is a completely within county analysis (it does not contain information across counties and just stratify them). The multilevel model is more efficient because it borrows information across counties. The counties with fewer measurements have more variable estimates and larger higher standard errors because they have less information. The counties with more measurements have more information and therefore, the estimates are more precise and the standard errors are smaller.

Excercise 5

In this exercise, you will replicate a set of models as shown below. This illustration is based on the ideas shown in figure 12.4 in the book. Please replicate this illustration with a twist: for example, you may choose a different set of counties than the one shown below or the one shown in the book, or you may try something different. Please explain the graphs. What do the dashed red and black lines show? What is the blue line? What determines whether the blue line would be more similar to the red line or the black line? How is this figure compare to the one shown in chapter 12 of the book?

5.1 Complete pooling model

mdl_comple <- radon_mn |>
  lm(log.radon ~ first_floor, data = _)

# Collect the predictions
pred_comple <- mdl_comple |>
  predictions(
    newdata = datagrid(
      first_floor = c(0, 1)
    )
  ) |> 
  select(first_floor, estimate, conf.low, conf.high) |> 
  as_tibble( )|> 
  mutate(Model = "Complete pooling model")

# Create the dataset
data.comple <- expand_grid(first_floor = c(0, 1),
                           county = unique(radon_mn$county))

# tidy the results
pred_comple <- data.comple %>% 
  left_join(pred_comple, by = "first_floor") |> 
  filter(county %in% c("CLAY", "RICE", "SCOTT","AITKIN", 
                       "DAKOTA", "GOODHUE", "STEARNS","WASHINGTON",
                       "REDWOOD",
                       "HENNEPIN", "ST LOUIS"))

5.2 No pooling model

mdl_nopool <- radon_mn |>
  lm(log.radon ~ first_floor + county, data = _)

# Collect the predictions
pred_nopool <- mdl_nopool |>
  predictions(
    newdata = datagrid(
      first_floor = c(0, 1),
      county = unique
    )
  ) |> 
  select(county, first_floor, estimate, conf.low, conf.high) |> 
  as_tibble( )|> 
  mutate(Model = "No pooling model") |> 
  filter(county %in% c("CLAY", "RICE", "SCOTT","AITKIN", 
                       "DAKOTA", "GOODHUE", "STEARNS","WASHINGTON",
                       "REDWOOD",
                       "HENNEPIN", "ST LOUIS"))

5.3 Multilevel model

mdl_multilevel <- radon_mn |>
  lmer(log.radon ~ first_floor + (first_floor|county), data = _)

# Collect the predictions
pred_multi <- mdl_multilevel |>
  predictions(
    newdata = datagrid(
      first_floor = c(0, 1),
      county = unique
    )
  ) |> 
  select(county, first_floor, estimate, conf.low, conf.high) |> 
  as_tibble( )|> 
  mutate(Model = "Multilevel model") |> 
  filter(county %in% c("CLAY", "RICE", "SCOTT","AITKIN", 
                       "DAKOTA", "GOODHUE", "STEARNS","WASHINGTON",
                       "REDWOOD",
                       "HENNEPIN", "ST LOUIS"))

5.4 plot

radon_mn |> 
  filter(county %in% c("CLAY", "RICE", "SCOTT","AITKIN", 
                       "DAKOTA", "GOODHUE", "STEARNS","WASHINGTON", "REDWOOD",
                       "HENNEPIN", "ST LOUIS")) |>
  ggplot(aes(x = first_floor, y = log.radon)) +
  geom_point() + 
  geom_line(data = pred_nopool, aes(first_floor, estimate), color = "Black")  + 
  geom_line(data = pred_comple, aes(first_floor, estimate), color = "Red")  + 
  geom_line(data = pred_multi, aes(first_floor, estimate), color = "Blue")  + 
  facet_wrap(~county)

Answer

The dashed red line shows the complete pooling model, the black line shows the no pooling model and the blue line shows the multilevel model. The complete pooling model is the same for all counties because it does not take into account the county level. The no pooling model is completely different for each county because it does not take into account the information from the other counties. The multilevel model is a combination of the complete pooling model and the no pooling model. The multilevel model is more similar to the complete pooling model when there is fewer number of measurements per county because it relies on the information from the other counties which is similar to the complete pooling model. The multilevel model is more similar to the no pooling model when there is a higher number of measurements per county because it relies on the information from the county itself which is similar to the no pooling model.

Excercise 6

In this exercise you are going to add a fixed effect at the county level (level 2). This illustration is based on figure 12.6 in the book, and it shows the estimated county level coefficients, plotted vs. county level uranium levels, along with estimated multilevel regression line at level 2. After replicating this graph, please explain it briefly.

# The simple linear model for the general relationship between 
# Radon measurements and county level uranium
md.lm <- radon_mn |>
  lm(log.radon ~ u , data = _)

# The multilevel model for the radon as a function of level-2 uranium
md.lvl12 <- radon_mn |>
  lmer(log.radon ~ u +  (1|county), data = _)


# Predictions for the linear model
pred.lm <- predictions(
  md.lm,
  )

# predictions for the multilevel model
# stratified by the level of uranium
pred.lvl12 <- predictions(
  md.lvl12,
  by = "u"
)


pred.lvl12 |>
  ggplot(aes(x = u, xend = u, y = conf.low, yend = conf.high)) +
  geom_segment() +
  geom_point(aes(y = estimate)) +
  geom_line(
    data = pred.lm,
    aes(x = u, y = estimate),
    ) +
  labs(y = "ECounty level estimates",
       x = "Uranium (county level)")

Answer

The black line illustrates the changes in log.Radon concerning the change in uranium levels based on the complete pooling model. The point of each segment represents the mean of log.radon for a county with a specific level of uranium, in the multilevel model. The county-level coefficients approximately follow this line, but not exactly. The deviations of these coefficients from the line are roughly reflected by the standard errors of random intercepts at level 2.

Excercise 7

Using the simulation techniques in questions 1 and 2, generate a synthetic dataset that would be similar to the radon data, in the sense that if you fitted the three models to your simulated dataset, you would get similar parameter estimations to the ones you got when fitting your models to the real Radon dataset. Explore your synthetic data using similar visualizations, and comment on similarities or differences between the Radon dataset and your simulated dataset.

# Print model coefficients
print(mdl_multilevel)

## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: log.radon ~ first_floor + (first_floor | county)
##    Data: radon_mn
## REML criterion at convergence: 2168.325
## Random effects:
##  Groups   Name        Std.Dev. Corr 
##  county   (Intercept) 0.3487        
##           first_floor 0.3436   -0.34
##  Residual             0.7462        
## Number of obs: 919, groups:  county, 85
## Fixed Effects:
## (Intercept)  first_floor  
##      1.4628      -0.6811

# Fixed intercept
alpha_0 <- 1.46277

# Fixed slope
beta_0 <- -0.6810982

# Residuals st.dev (level 1) 
sig  <- 0.7461543

# Random intercept st.dev (level 2) 
sig_u <- 0.3487

# Random slope st.dev (level 2)
sig_v <- 0.3436

# corr for random effects
cor <- -0.34

# Simulate the data
fun_simulation <- function(){

df <- add_random(county = 85) %>% 
  add_random(radon = radon_mn_count$n.obs, .nested_in = "county") |>
  add_ranef(
    "county",
    u_i = sig_u,
    v_i = sig_v,
    .cors = cor
    ) |> mutate(
      a_i = alpha_0 + u_i, 
      b_i = beta_0  + v_i
    ) |>
  mutate(
    first_floor = rbinom(n(), 1, sum(radon_mn$first_floor)/length(radon_mn$first_floor)),
    epsilon = rnorm(n(), 0, sig),
    y = a_i + b_i * first_floor + epsilon
    ) |>
  mutate(
    radon = exp(y)) |> 
  select(county, first_floor, radon)


# Fit the models
## Complete pooling
mdl_comple <- df |>
  lm(log(radon) ~ first_floor, data = _)

## No pooling
mdl_nopool <- df |>
  lm(log(radon) ~ first_floor + county, data = _)

## Multilevel
mdl_multilevel <- df |>
  lmer(log(radon) ~ first_floor + (first_floor|county), data = _)

# Extract the coefficients
## Complete pooling
beta_comple <- mdl_comple |>
  tidy() |>
  filter(term == "first_floor") |>
  pull(estimate)
## No pooling
beta_nopool <- mdl_nopool |>
  tidy() |>
  filter(term == "first_floor") |>
  pull(estimate)
## Multilevel
beta_multilevel <- mdl_multilevel |>
  tidy() |>
  filter(term == "first_floor") |>
  pull(estimate)

# Return the coefficients
return(
  tibble(
    beta_comple = beta_comple,
    beta_nopool = beta_nopool,
    beta_multilevel = beta_multilevel
  )
)

}

# Run the simulation
sim <- map(1:500, ~ fun_simulation()) |> 
  bind_rows(.id = "simulation")

## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')

# Bias and RMSE
sim |> 
  reframe(
    bias_comple = mean(abs(beta_comple - beta_0)),
    bias_nopool = mean(abs(beta_nopool - beta_0)),
    bias_multilevel = mean(abs(beta_multilevel - beta_0)),
    rmse_comple = sqrt(mean((beta_comple - beta_0)^2)),
    rmse_nopool = sqrt(mean((beta_nopool - beta_0)^2)),
    rmse_multilevel = sqrt(mean((beta_multilevel - beta_0)^2))
  )

# Boxplots
sim |> 
  pivot_longer(
    cols = starts_with("beta"),
    names_to = "model",
    values_to = "beta"
  ) |> 
  ggplot(aes(x = model, y = beta, fill = model)) +
  geom_boxplot(alpha = 0.6) +
  scale_fill_viridis_d() + 
  geom_hline(yintercept = beta_0, color = "darkred", linetype = "dashed", linewidth = 1) +
  labs(y = "Beta", x = "Model") +
  theme_bw()

Answer

After running the simulation 500 times, we can see that the bias and RMSE are smaller for the multilevel model than for the other two models. The boxplots show that the multilevel model is the one that is closest to the true value of the slope, and also the one with the smallest variability. On the other hand, the complete pooling model and the no pooling model have the similar bias and RMSE indicating that they perform similarly under the data generating process of multilevel structure. However, this conclusion does not hold for all the possible data generating processes and does not mean the multilevel model is always the best model in terms of bias and RMSE.

Team 13 Lab Work

Qingyang Shi, Yiling Zhou, Jaroslav Krc, Yao Lin, and Nathan Ribeiro

Simulating the data

Exercise 1

1.1 Creating the models

1.2 Comparing the models

Answer

Exercise 2

2.1 Caterpillar plots

2.2 Spaghetti plots

2.3 Lattice plots

Answer

Exercise 3

3.1 Load data

3.2 Explore data

Excercise 4

Answer

Excercise 5

5.1 Complete pooling model

5.2 No pooling model

5.3 Multilevel model

5.4 plot

Answer

Excercise 6

Answer

Excercise 7

Answer