Inferencia no válida cuando las observaciones no son independientes

Aprendí en las estadísticas elementales que, con un modelo lineal general, para que las inferencias sean válidas, las observaciones deben ser independientes. Cuando se produce la agrupación, la independencia ya no puede conducir a una inferencia no válida a menos que esto se tenga en cuenta. Una forma de dar cuenta de dicha agrupación es mediante el uso de modelos mixtos. Me gustaría encontrar un conjunto de datos de ejemplo, simulado o no, que lo demuestre claramente. Intenté usar uno de los conjuntos de datos de muestra en el sitio de UCLA para analizar datos agrupados

> require(foreign)
> require(lme4)
> dt <- read.dta("http://www.ats.ucla.edu/stat/stata/seminars/svy_stata_intro/srs.dta")

> m1 <- lm(api00~growth+emer+yr_rnd, data=dt)
> summary(m1)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 740.3981    11.5522  64.092   <2e-16 ***
growth       -0.1027     0.2112  -0.486   0.6271    
emer         -5.4449     0.5395 -10.092   <2e-16 ***
yr_rnd      -51.0757    19.9136  -2.565   0.0108 * 


> m2 <- lmer(api00~growth+emer+yr_rnd+(1|dnum), data=dt)
> summary(m2)

Fixed effects:
             Estimate Std. Error t value
(Intercept) 748.21841   12.00168   62.34
growth       -0.09791    0.20285   -0.48
emer         -5.64135    0.56470   -9.99
yr_rnd      -39.62702   18.53256   -2.14

A menos que me falte algo, estos resultados son lo suficientemente similares como para no pensar que el resultado lm()no es válido. He mirado algunos otros ejemplos (por ejemplo, 5.2 del Centro de la Universidad de Bristol para el modelado multinivel ) y descubrí que los errores estándar tampoco son terriblemente diferentes (no estoy interesado en los efectos aleatorios del modelo mixto, pero vale la pena señalar que el ICC de la salida del modelo mixto es 0.42).

Por lo tanto, mis preguntas son 1) bajo qué condiciones los errores estándar serán notablemente diferentes cuando ocurra la agrupación, y 2) alguien puede proporcionar un ejemplo de dicho conjunto de datos (simulado o no).

En primer lugar, tiene razón, este conjunto de datos quizás no sea el mejor para comprender el modelo mixto. Pero veamos primero por qué

require(foreign)
dt <- read.dta("http://www.ats.ucla.edu/stat/stata/seminars/svy_stata_intro/srs.dta")

length(dt$dnum)          # 310
length(unique(dt$dnum))  # 187 
sum(table(dt$dnum)==1)   # 132

Verá que tiene 310 observaciones y 187 grupos, de los cuales 132 tienen solo una observación. Esto no significa que no debamos usar el modelado multinivel, sino que no obtendremos resultados muy diferentes como usted dijo.

Motivación de modelado multinivel

La motivación para usar el modelado multinivel comienza desde el diseño en sí, y no solo desde los resultados del análisis realizado. Por supuesto, el ejemplo más común es tomar múltiples observaciones de personas, pero para hacer las cosas más extremas para dar una comprensión más fácil de la situación, piense en preguntar a las personas de diferentes países del mundo acerca de sus ingresos. Los mejores ejemplos son aquellos que tienen mucha heterogeneidad, ya que tomar grupos que son homogéneos en el resultado del examen, por supuesto, no hará mucha diferencia.

Ejemplo

Entonces, simulemos algunos datos para aclarar las cosas, la simulación funciona mejor ya que con los datos de la vida real no es tan obvio. Imagina que tomas$10$ países y preguntas $100$individuos de cada país sobre sus ingresos yy algo más xque tiene un efecto positivo en los ingresos con coeficiente$0.5$.

set.seed(1)
I <- 100
J <- 10
n <- I*J
i <- rep(1:I, each=J)
j <- rep(1:J,I)
x <- rnorm(n,mean=0, sd=1)
beta0  <- 1000
beta1  <- 0.5
sigma2 <- 1
tau2   <- 200
u <- rep(rnorm(I,mean=0,sd=sqrt(tau2)),each=J)
y <- beta0 + beta1*x + u + rnorm(n,mean=0, sd=sqrt(sigma2))

Entonces, ejecutando un modelo lineal obtienes

> summary(lm(y~x))

Coefficients:
            Estimate Std. Error  t value Pr(>|t|)    
(Intercept) 999.8255     0.4609 2169.230   <2e-16 ***
x             0.5728     0.4456    1.286    0.199    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 14.57 on 998 degrees of freedom
Multiple R-squared:  0.001653,  Adjusted R-squared:  0.0006528 
F-statistic: 1.653 on 1 and 998 DF,  p-value: 0.1989

y concluyes que xno tiene efecto estadístico en y. Vea qué tan grande es el error estándar. Pero ejecutar un modelo de intercepción aleatoria

> summary(lmer(y~x + (1|i)))

Random effects:
 Groups   Name        Variance Std.Dev.
 i        (Intercept) 213.062  14.597  
 Residual               1.066   1.032  
Number of obs: 1000, groups:  i, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept) 999.8247     1.4600   684.8
x             0.4997     0.0327    15.3

verá cuánto ha cambiado el error estándar de la estimación. Al observar la parte del efecto aleatorio, vemos cómo se ha descompuesto la variabilidad: la mayor parte de la variabilidad en el ingreso es entre países, y dentro de los países las personas tienen ingresos más similares. En palabras simples, lo que sucedió aquí es que no se tiene en cuenta la agrupación del efecto de x"perderse" (si podemos usar este tipo de término), sino que al descomponer la variabilidad se encuentra lo que realmente se debe obtener.