Ejemplo multinivel vs regresión clásica

Vamos a utilizar datos de actividad radiológica del gas radón en el estado de Minnesota. Los datos se han tomado del libro Data Analysis Using Regression and Multilevel/Hierarchical Models

Cargamos los datos y definimos variables

setwd("/home/jose/Dropbox/Multinivel")
srrs2 <- read.table("data/ARM_Data/radon/srrs2.dat", header = T, sep = ",")
mn <- srrs2$state == "MN"
radon <- srrs2$activity[mn]
log.radon <- log(ifelse(radon == 0, 0.1, radon))
floor <- srrs2$floor[mn]  # 0 for basement, 1 for first floor
n <- length(radon)
y <- log.radon
x <- floor

# Asignamos un número a cada condado
county.name <- as.vector(srrs2$county[mn])
uniq <- unique(county.name)
J <- length(uniq)
county <- rep(NA, J)
for (i in 1:J) {
    county[county.name == uniq[i]] <- i
}

Ajustamos un modelo lineal sin tener en cuenta los condados.

lm.pooled <- lm(y ~ x)
coef(lm.pooled)

## (Intercept)           x 
##      1.3267     -0.6134

Ahora ajustamos un modelo lineal teniendo en cuenta los condados

Quitamos el intercept y así no es necesario asignar un condado de referencia

lm.unpooled <- lm(y ~ x + factor(county) - 1)
coef(lm.unpooled)[1:5]

##               x factor(county)1 factor(county)2 factor(county)3 
##         -0.7205          0.8405          0.8748          1.5287 
## factor(county)4 
##          1.5527

Gráficamente

Primero seleccionamos unos pocos condados

x.jitter <- x + runif(n, -0.05, 0.05)
display8 <- c(36, 1, 35, 21, 14, 71, 61, 70)  # counties to be displayed
y.range <- range(y[!is.na(match(county, display8))])

par(mfrow = c(2, 4), mar = c(4, 4, 3, 1), oma = c(1, 1, 2, 1))
for (j in display8) {
    plot(x.jitter[county == j], y[county == j], xlim = c(-0.05, 1.05), ylim = y.range, 
        xlab = "floor", ylab = "log radon level", cex.lab = 1.2, cex.axis = 1.1, 
        pch = 20, mgp = c(2, 0.7, 0), xaxt = "n", yaxt = "n", cex.main = 1, 
        main = uniq[j])
    axis(1, c(0, 1), mgp = c(2, 0.7, 0), cex.axis = 1.1)
    axis(2, seq(-1, 3, 2), mgp = c(2, 0.7, 0), cex.axis = 1.1)
    curve(coef(lm.pooled)[1] + coef(lm.pooled)[2] * x, lwd = 2, col = "blue", 
        add = TRUE)
    curve(coef(lm.unpooled)[j + 1] + coef(lm.unpooled)[1] * x, col = "red", 
        add = TRUE)
}

plot of chunk unnamed-chunk-5

La línea azul es el ajuste general, sin tener en cuenta los condados, mientras que la línea roja es el ajuste cuando se ha tenido en cuenta los diferentes condados.

Análisis multinivel

Con el análisis multinivel se llega a una solución intermedia entre el modelo que no tiene en cuenta la variable county y el modelo que las mete como variables indicadoras.

Ajuste modelo multinivel con lmer

library(lme4)

## Loading required package: Matrix

## Loading required package: lattice

## Attaching package: 'lme4'

## The following object(s) are masked from 'package:stats':
## 
## AIC, BIC

M1 <- lmer(y ~ x + (1 | county))
summary(M1)

## Linear mixed model fit by REML 
## Formula: y ~ x + (1 | county) 
##   AIC  BIC logLik deviance REMLdev
##  2179 2199  -1086     2164    2171
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  county   (Intercept) 0.108    0.328   
##  Residual             0.571    0.756   
## Number of obs: 919, groups: county, 85
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   1.4616     0.0516   28.34
## x            -0.6930     0.0704   -9.84
## 
## Correlation of Fixed Effects:
##   (Intr)
## x -0.288

Coeficientes del modelo

# Para los 5 condados seleccionados
coef(M1)[[1]][display8, ]

##    (Intercept)      x
## 36      1.8681 -0.693
## 1       1.1915 -0.693
## 35      1.0870 -0.693
## 21      1.6311 -0.693
## 14      1.8380 -0.693
## 71      1.4829 -0.693
## 61      1.1995 -0.693
## 70      0.8899 -0.693

Separando las partes fijas y mixtas

# fixed and random effects
fixef(M1)

## (Intercept)           x 
##       1.462      -0.693

ranef(M1)[[1]][display8, ]

## [1]  0.40647 -0.27009 -0.37457  0.16952  0.37644  0.02132 -0.26207 -0.57170

Gráficamente

a.hat.M1 <- coef(M1)$county[, 1]  # 1st column is the intercept
b.hat.M1 <- coef(M1)$county[, 2]  # 2nd element is the slope

par(mfrow = c(2, 4))
for (j in display8) {
    plot(x.jitter[county == j], y[county == j], xlim = c(-0.05, 1.05), ylim = y.range, 
        xlab = "floor", ylab = "log radon level", main = uniq[j], cex.lab = 1.2, 
        cex.axis = 1.1, pch = 20, mgp = c(2, 0.7, 0), xaxt = "n", yaxt = "n", 
        cex.main = 1.1)
    axis(1, c(0, 1), mgp = c(2, 0.7, 0), cex.axis = 1)
    axis(2, c(-1, 1, 3), mgp = c(2, 0.7, 0), cex.axis = 1)
    curve(coef(lm.pooled)[1] + coef(lm.pooled)[2] * x, lwd = 2, col = "blue", 
        add = TRUE)
    curve(coef(lm.unpooled)[j + 1] + coef(lm.unpooled)[1] * x, col = "red", 
        add = TRUE)
    curve(a.hat.M1[j] + b.hat.M1[j] * x, lwd = 2, lty = 2, col = "darkgreen", 
        add = TRUE)
}

plot of chunk unnamed-chunk-9

Línea azul : Modelo lineal sin variable condados (Mismo modelo para todos los condados)
Línea roja: Modelo añadiendo condados al modelo anterior como un factor. (Modelo con distinto término constante para cada condado e igual pendiente)
Línea verde: Modelo multinivel

Se observa que en los condados con poco tamaño muestral, el modelo multinivel se parece más al modelo lineal, mientras que en los condados con mayor tamaño muestral se parece más al modelo lineal añadiendo el factor como variable.

Es decir, para los condados con poco tamaño muestral, tiene más en cuenta los datos globales que los datos del condado. Para condados con mayor tamaño muestral tiene más en cuenta los datos del condado que los generales.