6
$\begingroup$

Gelman & Hill (pp. 255-259) demonstrate in R how to achieve a "complete-pooling regression", "no-pooling regression", and "partial-pooling regression".

I don't have their data to replicate what they did. But using the data below, I was wondering if my understanding of these 3 types of regression is correct?

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))
$\endgroup$
4
$\begingroup$

So I went ahead and generated some data to demonstrate that these work as expected.

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)

The data generating mechanism from the top down is as follows...

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

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

Let's take a look at complete, no, and partial pooling.

Complete Pooling

This should return the same as the sample mean of y. This assumes that all the data are generated from a single normal distribution, with some mean and variance. The complete pooling uses all the data to estimate that one mean.

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

No Pooling

In this scenario, we agree that the groups are distinct, but we estimate their means using the data from those groups.

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

So we estimate the group means fairly well.

Partial Pooling

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

As you can see, the estimates for the groups are partially pooled towards the population mean (they are slightly less extreme than the complete pooling model).

Here is some code to reproduce these results. The results are not exactly the same because I didn't set the random seed when I wrote this.

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)

EDIT:

Here is an example with a fixed effect.

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)


You will note that the effect estimates in the partial pooling model are pooled towards the complete pooling estimates. They are ever so slightly closer.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ @rnorouzian The bug is fixed. If the book your referencing has a particular example, please post a minimal example in your question. I can't access the pages from you link. The model formula in the lm call will depend on what is being modelled. $\endgroup$ – Demetri Pananos Jul 5 at 4:47
  • 1
    $\begingroup$ I'll update with a fill script. $\endgroup$ – Demetri Pananos Jul 5 at 4:55
  • 1
    $\begingroup$ @rnorouzian I removed the intercept just as a means of more easily demonstrating what the group means were. Leaving the intercept in the model is equally valid, but changes the coefficients from the estimated group means to the difference between a reference group's mean and the indicated group's mean. $\endgroup$ – Demetri Pananos Jul 5 at 5:07
  • $\begingroup$ I think that is the correct expression for the ICC. So far as fixed effects go, let me see if I have time later today $\endgroup$ – Demetri Pananos Jul 7 at 14:00
  • $\begingroup$ @rnorouzian I'm not particualtly interested in reading that article at the moment. I use x=0 because the output would essentially be the estimates of the means for the groups when x=0, similar to what is reported in the complete pooling model. $\endgroup$ – Demetri Pananos Jul 9 at 4:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.