En primer lugar, aquí están mis modelos:
lm1 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + (1|x10/x11), data=df, na.action=na.exclude)
lm2 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x12 + (1|x10/x11), data=df, na.action=na.exclude)
lm3 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12:x6 + x12:x7 + x12:x8 + x12:x9 + (1|x10/x11), data=df, na.action=na.exclude)
lm4 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12*x6 + x12*x7 + x12*x8 + x12*x9 + (1|x10/x11), data=df, na.action=na.exclude)
lm5 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12:x6 + x12:x9:x7 + x12:x9:x8 + (1|x10/x11), data=df, na.action=na.exclude)
lm6 <- lmer(log(y) ~ x1 + x2 + x3 + x4 + x5 + x12*x6 + x12*x9*x7 + x12*x9*x8 + (1|x10/x11), data=df, na.action=na.exclude)El conjunto de datos es bastante grande, así que está alojado aquí (csv guardado como txt): https://1drv.ms/t/s!AvjuPq_dbORpaySu39F0vPXbv5Q?e=ah7rsi
Me cuesta entender por qué la adición de términos de efecto principal en un modelo hace que se eliminen los niveles múltiples de los términos de interacción, como se observa en la información de resumen del modelo.
Soy consciente de que x12 y x9 son nominales (es decir, 2 niveles y 9 niveles respectivamente), y que algunas de las interacciones que implican a ambas variables son deficientes en cuanto al rango. Este no es el origen del problema, se están eliminando de la información de resumen las interacciones que no son deficientes en el rango.
lm4 y lm6 son variantes de lm3 y lm5, respectivamente, en las que se han incluido términos de efecto principal junto con interacciones (es decir, lm3 y lm5 no contienen términos de efecto principal).
El problema es que estoy encontrando que algunos de los niveles de los términos de interacción, que son clave en mi aplicación, se eliminan de la información de resumen donde se incluyen los efectos principales (es decir, en lm4 y lm6). Este problema es especialmente grave en lm4 (véase más adelante). No puedo entender por qué la inclusión de los términos del efecto principal justificaría la eliminación de los niveles de interacción de la información de resumen.
NOTA: Utilizo el paquete lmerTest para reforzar la información resumida en mi análisis.
> summary(lm4)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(y) ~ x1 + x2 + x3 + x4 + x5 + x12*x6 + x12*x7 + x12*x8 + x12*x9 + (1|x10/x11)
Data: df
REML criterion at convergence: 229807.1
Scaled residuals:
Min 1Q Median 3Q Max
-6.4111 -0.4604 0.0917 0.5612 5.0303
Random effects:
Groups Name Variance Std.Dev.
x10:x11 (Intercept) 0.2678 0.5175
x10 (Intercept) 0.2735 0.5230
Residual 0.2036 0.4512
Number of obs: 132083, groups: PTR:PLT, 39282; PLT, 526
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -5.772e-01 2.775e-01 8.191e+02 -2.080 0.03785 *
x1 1.089e-01 1.125e-03 7.422e+04 96.801 < 2e-16 ***
x2 -1.881e-02 3.498e-04 8.952e+04 -53.775 < 2e-16 ***
x3 1.498e-02 1.624e-04 9.646e+04 92.204 < 2e-16 ***
x4 -1.305e-02 1.144e-03 4.606e+03 -11.407 < 2e-16 ***
x5 1.888e-05 2.370e-05 4.689e+02 0.797 0.42604
x12NF -4.808e-01 1.789e-01 3.366e+03 -2.687 0.00724 **
x6 7.401e-03 7.154e-03 4.768e+02 1.035 0.30140
x7 1.139e-04 4.641e-05 9.878e+04 2.455 0.01411 *
x8 -8.096e-04 1.180e-04 6.469e+04 -6.858 7.03e-12 ***
x9BW -2.252e-01 4.507e-02 4.245e+04 -4.997 5.84e-07 ***
x9BY 2.444e-01 1.479e-01 3.051e+04 1.653 0.09838 .
x9LT -2.574e-01 1.174e-01 4.461e+04 -2.192 0.02837 *
x9MR -4.036e-01 1.536e-01 3.875e+04 -2.627 0.00862 **
x9PW 1.225e-01 1.421e-01 3.486e+04 0.862 0.38874
x9RT -1.846e-01 1.064e-01 3.850e+04 -1.735 0.08269 .
x9SB -2.126e-01 2.381e-02 4.841e+04 -8.930 < 2e-16 ***
x9SW -5.054e-01 9.772e-02 5.010e+04 -5.172 2.33e-07 ***
x12NF:x6 2.415e-03 8.312e-03 4.736e+02 0.291 0.77150
x12NF:x7 -2.817e-04 4.660e-05 9.853e+04 -6.045 1.50e-09 ***
x12NF:x8 8.554e-04 1.191e-04 6.624e+04 7.183 6.91e-13 ***
x12NF:x9BW 9.307e-03 4.866e-02 4.112e+04 0.191 0.84834
x12NF:x9LT 2.630e-01 1.289e-01 4.242e+04 2.041 0.04128 *
x12NF:x9MR 1.404e-01 1.639e-01 3.842e+04 0.856 0.39179
x12NF:x9RT 7.611e-02 1.204e-01 3.696e+04 0.632 0.52719
x12NF:x9SB 3.738e-01 2.644e-02 4.469e+04 14.136 < 2e-16 ***
x12NF:x9SW 5.014e-01 1.027e-01 4.789e+04 4.881 1.06e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation matrix not shown by default, as p = 27 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
fit warnings:
fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
