Demostración de regresión con pooling completo, sin pooling y con pooling parcial en R

Question

Gelman y Hill (pp. 255-259) demostrar en R cómo conseguir un "regresión completa" , "regresión sin piscina" y "regresión parcial" .

No tengo sus datos para replicar lo que hicieron. Pero utilizando los datos a continuación, me preguntaba si mi comprensión de estos 3 tipos de regresión es correcta?

library(lme4)                                   # needed for partial-pooling
group <- gl(2, 50, labels = c("Ctl","Trt"))     # group indicator
    y <- c(Ctl = rnorm(50), Trt = rnorm(50, 1)) # dependent variable

complete_pooling <- lm(y ~ 1)
      no_pooling <- lm(y ~ group)
 partial_pooling <- lmer(y ~ 1 + (1|group))

Answer 1

Así que seguí adelante y generé algunos datos para demostrar que funcionan como se esperaba.

library(tidyverse)
library(lme4)

if(!require(modelr)){
  install.packages('modelr')
}
library(modelr)

pop_mean<-10
n_groups<-4
groups<-gl(n_groups, 20)
Z<-model.matrix(~groups-1)
group_means<-rnorm(n_groups, 0, 2.5)

y<- pop_mean + Z%*%group_means + rnorm(length(groups), 0, 0.5)
d<-tibble(y, groups)

El mecanismo de generación de datos de arriba abajo es el siguiente...

$$ \theta_i \sim \mathcal{N}(10, 2.5) $$

$$y_{i,j} \sim \mathcal{N}(\theta_i, 0.5) $$

Echemos un vistazo a la agrupación completa, la no agrupación y la agrupación parcial.

Puesta en común completa

Esto supone que todos los datos se generan a partir de una única distribución normal, con cierta media y varianza. La agrupación completa utiliza todos los datos para estimar esa media.

complete_pooling<-lm(y~1, data = d)
summary(complete_pooling)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    9.264      0.214   43.29   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.914 on 79 degrees of freedom

Sin agrupación

En este escenario, estamos de acuerdo en que los grupos son distintos, pero estimamos sus medias utilizando los datos de esos grupos.

no_pooling<-lm(y~groups-1, data = d) #remove the intercept from the model
summary(no_pooling)

Coefficients:
        Estimate Std. Error t value Pr(>|t|)    
groups1   6.2116     0.1045   59.44   <2e-16 ***
groups2  10.9183     0.1045  104.48   <2e-16 ***
groups3  10.5156     0.1045  100.63   <2e-16 ***
groups4   9.4088     0.1045   90.04   <2e-16 ***
---

group_means + pop_means # pretty close
>>> 6.311974 10.878787 10.354225  9.634138

Así que estimamos bastante bien las medias de los grupos.

Puesta en común parcial

partial_pooling<-lmer(y~ 1 + 1|groups, data = d)

summary(partial_pooling)

Random effects:
 Groups   Name        Variance Std.Dev.
 groups   (Intercept) 4.5362   2.1298  
 Residual             0.2184   0.4673  
Number of obs: 80, groups:  groups, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept)    9.264      1.066   8.688

modelr::data_grid(d, groups) %>% modelr::add_predictions(partial_pooling)

# A tibble: 4 x 2
  groups  pred
  <fct>  <dbl>
1 1       6.22
2 2      10.9 
3 3      10.5 
4 4       9.41

Como puede ver, las estimaciones de los grupos están parcialmente agrupadas hacia la media de la población (son ligeramente menos extremas que el modelo de agrupación completa).

He aquí un código para reproducir estos resultados. Los resultados no son exactamente lo mismo porque no puse la semilla aleatoria cuando escribí esto.

library(tidyverse)
library(lme4)

if(!require(modelr)){
  install.packages('modelr')
}
library(modelr)

#Generate data
set.seed(123)
pop_mean<-10
n_groups<-4
groups<-gl(n_groups, 20)
Z<-model.matrix(~groups-1)
group_means<-rnorm(n_groups, 0, 2.5)

y<- pop_mean + Z%*%group_means + rnorm(length(groups), 0, 0.5)

d = tibble(y, groups)

complete_pooling<-lm(y~1, data = d)
no_pooling<-lm(y~groups-1, data = d)
partial_pooling<-lmer(y~ 1 + 1|groups, data = d)

modelr::data_grid(d, groups) %>% modelr::add_predictions(partial_pooling)

EDITAR:

He aquí un ejemplo con un efecto fijo.

library(tidyverse)
library(lme4)

if(!require(modelr)){
  install.packages('modelr')
}
library(modelr)

#Generate data
set.seed(123)
pop_mean<-10
n_groups<-4
groups<-gl(n_groups, 20)
x<-rnorm(length(groups))
Z<-model.matrix(~groups-1)
group_means<-rnorm(n_groups, 0, 2.5)

y<- pop_mean + 2*x + Z%*%group_means + rnorm(length(groups), 0, 0.5)

d = tibble(y, groups,x)

complete_pooling<-lm(y~x, data = d)
no_pooling<-lm(y~groups + x -1, data = d)
partial_pooling<-lmer(y~ x + 1 + 1|groups, data = d)

modelr::data_grid(d, groups,x=0) %>% modelr::add_predictions(partial_pooling)

Observará que las estimaciones del efecto en el modelo de agrupación parcial se agrupan hacia las estimaciones de agrupación completa. Se aproximan ligeramente.