Comparar las pendientes de la regresión lineal entre modelos no anidados con tamaños de conjunto de datos diferentes

Question

Me gustaría probar si las pendientes de dos modelos de regresión lineal difieren. Sin embargo, la advertencia es que una de las regresiones se ajusta a un subconjunto de los datos, y la otra se ajusta a todo el conjunto de datos. Los dos modelos utilizan las mismas variables predictoras y no son anidados. las pruebas de razón de verosimilitud requieren que los modelos comparados estén anidados y ajustados al mismo tamaño de conjunto de datos. Se ha sugerido el uso de AIC y BIC para comparar modelos no anidados, pero no me interesa cuál modelo se ajusta mejor, solo si las pendientes difieren significativamente entre sí.

Ejemplo (R):

library('lmtest')   
df <- data.frame(time=c(1,2,3,4,5,6,7,8,9,10), response=c(1,2,4,3,7,7,6,6,7,10))
lm1 <- lm(response ~ time, data=df)
lm2 <- lm(response ~ time, data=df[1:6,])
lrtest(lm1, lm2) # este enfoque no funciona porque los conjuntos de datos son de tamaños diferentes

¿Qué otro enfoque podría funcionar aquí? ¿Puedo simplemente mirar los intervalos de confianza de las pendientes y verificar si se superponen? ¡Gracias por tu ayuda!

Answer 1

Bienvenido aquí.

Una solución simple es definir una nueva variable binaria que tome el valor $1$ para el subconjunto de los datos y $0$ para todas las demás observaciones que no están incluidas en el subconjunto. Luego puedes interactuar la pendiente en la que estás interesado con la variable binaria para ver si encuentras diferencias estadísticamente significativas en la pendiente entre el subconjunto y el conjunto de datos completo:

df$sample = 0 # 'sample' toma el valor 0 para todo el conjunto de datos
df[1:6, "sample"] = 1 # 'sample' toma el valor 1 para el subconjunto
linear_regression_model_with_interaction = lm(response ~ time + time:sample, data = df)
summary(linear_regression_model_with_interaction)

Los resultados se ven así. En los datos de ejemplo que presentaste, encontramos diferencias estadísticas en la pendiente estimada en el nivel del 10 por ciento (ignorando la posible heterocedasticidad, etc., aquí, o si sería sabio incluir también un efecto principal para la variable 'sample'):

Llamada:
lm(formula = response ~ time + time:sample, data = df)

Residuos:
     Min       1Q   Mediana       3Q      Máx. 
-1.61401 -0.63063  0.03804  0.58723  1.30236 

Coeficientes:
            Estimado Error estándar Valor t Pr(>|t|)    
(Intercept)  -0.5893     0.9292  -0.634 0.546128    
time          0.9287     0.1230   7.548 0.000132 ***
time:sample   0.3721     0.1632   2.281 0.056591 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error estándar residual: 1.027 en 7 grados de libertad
R-cuadrado múltiple:  0.8915,    R-cuadrado ajustado:  0.8606 
Estadístico F: 28.77 con 2 y 7 DF,  valor p: 0.0004201

PD: La posible heterocedasticidad se puede tener en cuenta a través del error estándar robusto a la heterocedasticidad:

library('lmtest')
library('sandwich')
coeftest(linear_regression_model_with_interaction, vcov = vcovHC(linear_regression_model_with_interaction, type = "HC0"))

Pero ten en cuenta que los errores estándar robustos podrían no ser válidos con muestras pequeñas, este podría ser el motivo por el cual el valor p baja aquí (por supuesto, la heterocedasticidad puede sesgar tu estimación de la varianza en ambas direcciones pero usualmente el error estándar se vuelve más grande).

Prueba t de coeficientes:

             Estimado Error estándar Valor t  Pr(>|t|)    
(Intercepto) -0.589293   0.357539 -1.6482   0.14330    
time         0.928694   0.079517 11.6792 7.625e-06 ***
time:sample  0.372132   0.109051  3.4124   0.01125 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1