#INTRODUCCIÓN
library(faraway)## Warning: package 'faraway' was built under R version 4.2.3Aunque el modelo lineal no es necesariamente inadecuado, la inclusión del bloque como un factor fijo puede también utilizarse. En este sentido el modelo mixto equivalente sería mucho más parsimonioso en término del número de parámetros.
# nos interesa analizar el efecto del factor (dieta) sobre la respuesta (ganancia de peso) una vez que se haya tenido en cuenta la variabilidad atribuible al bloque (camada).
library(nlme)
lme.rabbit1 <- lme(gain~ treat, random=~1|block, data=rabbit)
anova(lme.rabbit1)## numDF denDF F-value p-value
## (Intercept) 1 15 590.5297 <.0001
## treat 5 15 3.2818 0.0336El tratamiento es significativo, pues su p valor es 0.03 < 0.05.
summary(lme.rabbit1)## Linear mixed-effects model fit by REML
## Data: rabbit
## AIC BIC logLik
## 166.3569 175.7813 -75.17844
##
## Random effects:
## Formula: ~1 | block
## (Intercept) Residual
## StdDev: 4.657853 3.175519
##
## Fixed effects: gain ~ treat
## Value Std.Error DF t-value p-value
## (Intercept) 39.53540 2.130338 15 18.558278 0.0000
## treatb -2.50718 2.208700 15 -1.135139 0.2741
## treatc -0.18407 2.208700 15 -0.083340 0.9347
## treatd -0.88516 2.208700 15 -0.400759 0.6942
## treate -5.64602 2.208700 15 -2.556263 0.0219
## treatf 2.81003 2.208700 15 1.272254 0.2227
## Correlation:
## (Intr) treatb treatc treatd treate
## treatb -0.518
## treatc -0.518 0.500
## treatd -0.518 0.500 0.500
## treate -0.518 0.500 0.500 0.500
## treatf -0.518 0.500 0.500 0.500 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.3794203 -0.5213453 -0.1701517 0.6456859 1.4148990
##
## Number of Observations: 30
## Number of Groups: 10El intercepto y el tratamiento e son estadísticamente significativos al 5%.
#Otros modelos mixtos aplicables
#Comparamos los modelos.
#un modelo sin efecto aleatorio
lme.rabbit0 <- gls(gain~ treat, data=rabbit) #Se utilza gls para poder compararlo, el resultados es similar a un modelo lineal
#un modelo con efecto aleatorio
lme.rabbit1 <- lme(gain~ treat, data=rabbit, random=~1|block)
#En ambos fue utilizando un estimador REML (bajo la función gls, que permite comparar los parámetros de modelos. Si se hubiera hecho con lm, no es posible compararlos)
anova(lme.rabbit0, lme.rabbit1)## Model df AIC BIC logLik Test L.Ratio p-value
## lme.rabbit0 1 7 174.2621 182.5085 -80.13107
## lme.rabbit1 2 8 166.3569 175.7813 -75.17844 1 vs 2 9.905255 0.0016El contraste restringido vs no restringido muestra un p valor de 0.0016 < 0.05. Tenemos evidencia empírica suficiente para rechazar H0 al 5%. Existen diferencias entre el modelo general y el modelo de bloques aleatorizados. Se deben incluir, por tanto, estos bloques aleatorizados en el modelo.
#mejor opción comparar modelos anidados como en el caso de los efectos
#aleatorios.
lme.rabbit2 <- lme(gain~ 1, data=rabbit, random=~1|block, method="ML")
lme.rabbit3 <- lme(gain~ treat, data=rabbit, random=~1|block, method="ML")
anova(lme.rabbit2, lme.rabbit3)## Model df AIC BIC logLik Test L.Ratio p-value
## lme.rabbit2 1 3 188.8307 193.0343 -91.41536
## lme.rabbit3 2 8 183.7730 194.9825 -83.88648 1 vs 2 15.05776 0.0101Comparando los dos modelos de bloques aleatorios, se puede indicar que el tratamiento es conjuntamente significativo, es decir, produce medias diferentes cada uno de ellos.
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
)## # A tibble: 58 × 3
## package loadedversion source
## <chr> <chr> <chr>
## 1 boot 1.3-28 CRAN (R 4.2.2)
## 2 bslib 0.4.2 CRAN (R 4.2.2)
## 3 cachem 1.0.6 CRAN (R 4.2.2)
## 4 callr 3.7.3 CRAN (R 4.2.3)
## 5 cli 3.6.0 CRAN (R 4.2.2)
## 6 crayon 1.5.2 CRAN (R 4.2.3)
## 7 devtools 2.4.5 CRAN (R 4.2.3)
## 8 digest 0.6.31 CRAN (R 4.2.2)
## 9 ellipsis 0.3.2 CRAN (R 4.2.2)
## 10 evaluate 0.20 CRAN (R 4.2.2)
## # ℹ 48 more rows