Modelos mixtos VyN
Valentina Melo y Natalia Menjura
2023-06-26
PRIMER EJERCICIO - NIVEL DE CAPTACIÓN DE CO2 EN PLANTAS
Se preparan las librerías a emplear en el modelo.
- Las más importantes para modelos lineales mixtos son el paquete “lme4” y “broom.mixed” que es para tener resumenes del modelo.
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(broom.mixed)
library(lme4)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpackAñadiendo la base de datos
- Se trabajará con la base de datos de R de “CO2” que contiene datos de un experimento sobre plantas que son de Quebec o de Mississippi que tienen distintos tipos de tratamiento en zonas enfriadas y no enfriadas. Se evalúa las plantas a distintas concentraciones de dioxido de carbono y miden el nivel de captación de CO2.
data("CO2")
View(CO2)
head(CO2)## Grouped Data: uptake ~ conc | Plant
## Plant Type Treatment conc uptake
## 1 Qn1 Quebec nonchilled 95 16.0
## 2 Qn1 Quebec nonchilled 175 30.4
## 3 Qn1 Quebec nonchilled 250 34.8
## 4 Qn1 Quebec nonchilled 350 37.2
## 5 Qn1 Quebec nonchilled 500 35.3
## 6 Qn1 Quebec nonchilled 675 39.2Generando ruido en la base de datos del nivel de captación de CO2
- Añadimos una variación en los datos de captación de CO2. La nueva variable la llamaremos entonces uptake2 (variable de captación de CO2 de la base de datos (uptake) más variación (ruido)).
set.seed(123)
ruido = runif(84, 0.02, 0.04)
uptake2 = CO2$uptake+ruido
data.frame(uptake2)## uptake2
## 1 16.025752
## 2 30.435766
## 3 34.828180
## 4 37.237660
## 5 35.338809
## 6 39.220911
## 7 39.730562
## 8 13.637848
## 9 27.331029
## 10 37.129132
## 11 41.839137
## 12 40.629067
## 13 41.433551
## 14 44.331453
## 15 16.222058
## 16 32.437996
## 17 40.324922
## 18 42.120841
## 19 42.926558
## 20 43.939090
## 21 45.537791
## 22 14.233856
## 23 24.132810
## 24 30.339885
## 25 34.633114
## 26 32.534171
## 27 35.430881
## 28 38.731883
## 29 9.325783
## 30 27.322942
## 31 35.039260
## 32 38.838046
## 33 38.633814
## 34 37.535909
## 35 42.420492
## 36 15.129556
## 37 21.035169
## 38 38.124328
## 39 34.026364
## 40 38.924633
## 41 39.622856
## 42 41.428291
## 43 10.628274
## 44 19.227377
## 45 26.223049
## 46 30.022776
## 47 30.924661
## 48 32.429319
## 49 35.525319
## 50 12.037157
## 51 22.020917
## 52 30.628844
## 53 31.835978
## 54 32.422438
## 55 31.131219
## 56 31.524131
## 57 11.322551
## 58 19.435066
## 59 25.837901
## 60 27.927489
## 61 28.533302
## 62 28.121897
## 63 27.827679
## 64 10.525488
## 65 14.936293
## 66 18.128970
## 67 18.936201
## 68 19.536248
## 69 22.235887
## 70 21.928797
## 71 7.735090
## 72 11.432584
## 73 12.334204
## 74 13.020012
## 75 12.529506
## 76 13.724402
## 77 14.427596
## 78 10.632255
## 79 18.027036
## 80 17.922223
## 81 17.924872
## 82 17.933361
## 83 18.928353
## 84 19.935764Graficando el experimento
Observamos la captación de CO2 de las plantas a distintas concentraciones de forma gráfica. Para esto establecemos en el eje “X” la concentracion de CO2 y en el eje “Y” el nivel de captación. Para diferenciar los valores, añadimos color según el tipo y añadimos una forma según el tipo de tratamiento. Finalmente añadimos las lineas de proyección de los tratamientos, en donde, las diferenciamos por grupo de planta y por tratamiento. Esto con el fin de realizar estimaciones del tipo de comportamiento de otras plantas para el experimento.
En la gráfica se observa variación donde las plantas de Quebec presentan mayor capacidad de captación de CO2 respecto a las plantas de Mississippi.
ggplot(CO2, aes(x = conc, y = uptake2, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment)) +
theme_bw()
Solucionando el modelo de tipo aleatorio
- Se empleará el paquete lme4 para lidiar con el modelo aleatorio. Para esto primero se desarrolla un modelo normal, lineal simple, donde se mira la captación (uptake2) explicada con la concentración (conc) más el tipo y el tratamiento. Para agregar el factor aleatorio empleamos el modelo lmer y con esto, después del tratamiento, le añadimos el factor aleatorio de la planta (1|Plant, esto hace que cada una de las plantas tenga un intercepto distinto).
Fit1 <- lm(uptake2 ~ conc + Type + Treatment, data = CO2)
Fit2 <- lmer(uptake2 ~ conc + Type + Treatment + (1|Plant), data = CO2)Observando los modelos
- Se emplea el paquete “broom” que nos permite transformar los modelos en datos ordenados. En este caso usamos el “broom::glance” para extraer las metricas de la calidad del modelo y aplicarlas al mismo para obtener un data frame de una única fila. Por otra parte el “broom::tidy” nos permite construir un resumen con los hallazgos estadísticos del modelo. Con la función “tidy” el resumen nos muestra para cada uno de los factores cual va a ser el estimado, el error estándar, el estadístico y el p - valor.
broom::glance(Fit1)## # A tibble: 1 × 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.684 0.672 6.19 57.7 5.88e-20 3 -270. 551. 563.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>broom::tidy(Fit1)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 29.3 1.54 19.0 1.94e-31
## 2 conc 0.0177 0.00230 7.72 2.88e-11
## 3 TypeMississippi -12.7 1.35 -9.37 1.66e-14
## 4 Treatmentchilled -6.86 1.35 -5.08 2.47e- 6Empleando el modelo lineal mixto
Emplearemos la función de “broom.mixed” para generar los mismos resúmenes que para los de la función “broom” siendo en este caso para los del modelo lineal mixto (el generado con el modelo lmer).
A diferencia con el “broom::glance” del modelo lineal simple, para la función de “broom.mixed::glance” del modelo lineal mixto no se tiene el valor de R cuadrado.
Proyectando el modelo, se pueden observar las variables y se comprueba que para cada planta y para cada observación existe una estimación. Se tiene también que se divide entre efectos fijos y efectos aleatorios.
broom.mixed::glance(Fit2)## # A tibble: 1 × 7
## nobs sigma logLik AIC BIC REMLcrit df.residual
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 84 6.00 -272. 556. 571. 544. 78broom.mixed::tidy(Fit2)## # A tibble: 6 × 6
## effect group term estimate std.error statistic
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 fixed <NA> (Intercept) 29.3 1.72 17.0
## 2 fixed <NA> conc 0.0177 0.00223 7.96
## 3 fixed <NA> TypeMississippi -12.7 1.64 -7.72
## 4 fixed <NA> Treatmentchilled -6.86 1.64 -4.18
## 5 ran_pars Plant sd__(Intercept) 1.71 NA NA
## 6 ran_pars Residual sd__Observation 6.00 NA NAbroom.mixed::tidy(Fit2) %>% View()Añadiendo logaritmo al modelo y mirando el comportamiento del tipo con la interacción
- Añadiremos al modelo la función de logaritmo I(log), esto debido a que en la gráfica se observa el comportamiento de los datos de tipo logarítmica y se observa en la gráfica que existe interacción entre el tipo con el tratamiento por lo cual se puede asumir que existen diferencias.
Fit1 <- lm(uptake2 ~ I(log(conc)) + Type:Treatment, data = CO2)
Fit2 <- lmer(uptake2 ~ I(log(conc)) + Type:Treatment + (1 | Plant), data = CO2)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientResumiendo los modelos
- Se debe tener presente para diferenciar los modelos. Si bien las funciones se interpretan de la misma mandera, tener presente el tipo de modelo que se está desarrollando nos ayuda a realizar predicciones y estimaciones más exactas según el tipo de experimento que se está desarrollando.
broom::glance(Fit1)## # A tibble: 1 × 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.823 0.814 4.67 91.7 7.02e-29 4 -246. 504. 519.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>broom::tidy(Fit1)## # A tibble: 6 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -33.5 4.08 -8.22 3.21e-12
## 2 I(log(conc)) 8.48 0.678 12.5 2.04e-20
## 3 TypeQuebec:Treatmentnonchilled 19.5 1.44 13.6 2.61e-22
## 4 TypeMississippi:Treatmentnonchilled 10.1 1.44 7.04 6.35e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 9.95e-18
## 6 TypeMississippi:Treatmentchilled NA NA NA NAbroom.mixed::glance(Fit2)## # A tibble: 1 × 7
## nobs sigma logLik AIC BIC REMLcrit df.residual
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 84 4.50 -241. 496. 513. 482. 77broom.mixed::tidy(Fit2)## # A tibble: 7 × 6
## effect group term estimate std.error statistic
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 fixed <NA> (Intercept) -33.5 4.02 -8.34
## 2 fixed <NA> I(log(conc)) 8.48 0.654 13.0
## 3 fixed <NA> TypeQuebec:Treatmentnonchilled 19.5 1.83 10.6
## 4 fixed <NA> TypeMississippi:Treatmentnonch… 10.1 1.83 5.52
## 5 fixed <NA> TypeQuebec:Treatmentchilled 15.9 1.83 8.69
## 6 ran_pars Plant sd__(Intercept) 1.47 NA NA
## 7 ran_pars Residual sd__Observation 4.50 NA NAbroom.mixed::tidy(Fit2) %>% View()SEGUNDO EJERCICIO - PESO DE POLLOS
Agregando la base de datos y graficando el experimento
- En la base de datos se tienen las variables de tiempo, peso, dieta y la identidad del pollo. Para graficar el experimento empleamos la función “ggplot” (igualmente empleada en el experimento 1), donde se establecio en el eje “X” la variable de tiempo y en el eje “Y” la variable de peso. Para diferenciar, se establecio el color de cada punto en función de la dieta y se proyectaron las líneas de gráfico en función de la dieta y por grupo de identidad de pollo.
data("Chickweight")## Warning in data("Chickweight"): data set 'Chickweight' not foundview(ChickWeight)
ggplot(ChickWeight, aes(x = Time, y = weight)) +
geom_point(aes(color = Diet)) +
geom_path(aes(color = Diet, group = Chick))
Proyectando la base de datos del experimento
- Primero se extrajo de la base de datos el tolat de pollos del experimento y con base en esta se extrajo la cantidad de determinantes únicos, obteniendo que se cuenta con 50 pollos distintos. Se observa en la base de datos que todos los valores de peso son positivos y no se cuenta con valores decimales.
ChickWeight$Chick## [1] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3
## [26] 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 5 5
## [51] 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7
## [76] 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9
## [101] 9 9 9 9 9 9 9 10 10 10 10 10 10 10 10 10 10 10 10 11 11 11 11 11 11
## [126] 11 11 11 11 11 11 12 12 12 12 12 12 12 12 12 12 12 12 13 13 13 13 13 13 13
## [151] 13 13 13 13 13 14 14 14 14 14 14 14 14 14 14 14 14 15 15 15 15 15 15 15 15
## [176] 16 16 16 16 16 16 16 17 17 17 17 17 17 17 17 17 17 17 17 18 18 19 19 19 19
## [201] 19 19 19 19 19 19 19 19 20 20 20 20 20 20 20 20 20 20 20 20 21 21 21 21 21
## [226] 21 21 21 21 21 21 21 22 22 22 22 22 22 22 22 22 22 22 22 23 23 23 23 23 23
## [251] 23 23 23 23 23 23 24 24 24 24 24 24 24 24 24 24 24 24 25 25 25 25 25 25 25
## [276] 25 25 25 25 25 26 26 26 26 26 26 26 26 26 26 26 26 27 27 27 27 27 27 27 27
## [301] 27 27 27 27 28 28 28 28 28 28 28 28 28 28 28 28 29 29 29 29 29 29 29 29 29
## [326] 29 29 29 30 30 30 30 30 30 30 30 30 30 30 30 31 31 31 31 31 31 31 31 31 31
## [351] 31 31 32 32 32 32 32 32 32 32 32 32 32 32 33 33 33 33 33 33 33 33 33 33 33
## [376] 33 34 34 34 34 34 34 34 34 34 34 34 34 35 35 35 35 35 35 35 35 35 35 35 35
## [401] 36 36 36 36 36 36 36 36 36 36 36 36 37 37 37 37 37 37 37 37 37 37 37 37 38
## [426] 38 38 38 38 38 38 38 38 38 38 38 39 39 39 39 39 39 39 39 39 39 39 39 40 40
## [451] 40 40 40 40 40 40 40 40 40 40 41 41 41 41 41 41 41 41 41 41 41 41 42 42 42
## [476] 42 42 42 42 42 42 42 42 42 43 43 43 43 43 43 43 43 43 43 43 43 44 44 44 44
## [501] 44 44 44 44 44 44 45 45 45 45 45 45 45 45 45 45 45 45 46 46 46 46 46 46 46
## [526] 46 46 46 46 46 47 47 47 47 47 47 47 47 47 47 47 47 48 48 48 48 48 48 48 48
## [551] 48 48 48 48 49 49 49 49 49 49 49 49 49 49 49 49 50 50 50 50 50 50 50 50 50
## [576] 50 50 50
## 50 Levels: 18 < 16 < 15 < 13 < 9 < 20 < 10 < 8 < 17 < 19 < 4 < 6 < 11 < ... < 48ChickWeight$Chick %>% length()## [1] 578ChickWeight$Chick %>% unique() %>% length()## [1] 50view(ChickWeight)Modelo del experimento
Como no se tienen valores decimales en cuanto al peso de los pollos, se puede estimar que el modelo es un glm de tipo poisson. En este caso establecemos que el peso es explicado con la dieta en interacción con el tiempo.
Emplearemos la función “tidy” para extraer el resumen del modelo estadístico. Esta función nos muestra en general que se presenta el aumento del peso en el tiempo y va a existir una interacción respecto a la dieta.
Fit1_Poisson <- glm(weight ~ Diet:Time, data = ChickWeight, family = poisson)
tidy(Fit1_Poisson)## # A tibble: 5 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.86 0.00931 415. 0
## 2 Diet1:Time 0.0656 0.000722 90.8 0
## 3 Diet2:Time 0.0754 0.000768 98.1 0
## 4 Diet3:Time 0.0862 0.000729 118. 0
## 5 Diet4:Time 0.0823 0.000756 109. 0Modelo lineal mixto
Gracias a la gráfica, se puede observar que existen pollos que aunque cuentan con el mismo tipo de dieta, presentan comportamientos diferentes. Con base a esto, gracias a estas medidas repetidas, nos permite ajustar el modelo con la función “glmer” donde podemos mirar el comportamiento de los pesos de los pollos explicado con la interacción entre el peso y el tiempo, más una nueva variable que en este caso serán los errores aleatorios de los individuos (1|Chick).
Establecido el modelo, emplearemos la función “tidy” para comparar los resumenes de cada uno de los modelos, donde se observa que existen unas diferencias muy pequeñas en el intercepto, siendo estas menores cuando se consideran los individuos (relacionado al modelo glmer para modelos mixtos).
Los resultados para el modelo lineal mixto detalla que la identidad de los pollos recibe alguna variabilidad la cual impica que existen diferencias para el experimento.
Observando los resultados obtenidos para el modelo lineal mixto se evidencia que la pendiente más alta se encuentra con el tipo de dieta 3 seguido de la dieta 4, 2 y finalmente 1. Por otra parte, en el modelo lienal simple cambian los resultados siendo la dieta 2 la que obtuvo mejores resultados por lo cual se detalla el cambio en el modelo.
Fit2_Poisson <- glmer(weight ~ Diet:Time + (1|Chick), data = ChickWeight, family = poisson)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?broom::tidy(Fit1_Poisson)## # A tibble: 5 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.86 0.00931 415. 0
## 2 Diet1:Time 0.0656 0.000722 90.8 0
## 3 Diet2:Time 0.0754 0.000768 98.1 0
## 4 Diet3:Time 0.0862 0.000729 118. 0
## 5 Diet4:Time 0.0823 0.000756 109. 0broom.mixed::tidy(Fit2_Poisson)## # A tibble: 6 × 7
## effect group term estimate std.error statistic p.value
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 fixed <NA> (Intercept) 3.84 0.0315 122. 0
## 2 fixed <NA> Diet1:Time 0.0666 0.00105 63.7 0
## 3 fixed <NA> Diet2:Time 0.0749 0.00129 57.9 0
## 4 fixed <NA> Diet3:Time 0.0868 0.00123 70.5 0
## 5 fixed <NA> Diet4:Time 0.0763 0.00126 60.4 0
## 6 ran_pars Chick sd__(Intercept) 0.213 NA NA NAview(broom.mixed::tidy(Fit2_Poisson))
view(broom::tidy(Fit1_Poisson))Generando un plot_model de tipo efecto
- Con la función “sjPlot” se realiza el modelo del experimento para detallar los resultados obtenidos. En esta se observa gráficamente que los resultados en ganancia de peso respecto a dieta son más favorables para el tipo de dieta 3 al igual que se detalla la ganacia de peso de los pollos respecto al tiempo.
sjPlot::plot_model(Fit1_Poisson,type = "eff")## Package `effects` is not available, but needed for `ggeffect()`. Either
## install package `effects`, or use `ggpredict()`. Calling `ggpredict()`
## now.## $Diet
##
## $Time
RESUMIENDO
Tomar en cuenta la variabiblidad de un individuo frente a otro mejora el modelo para realizar predicciones del comportamiento de los datos.
Siempre se debe tener más de una medida por individuo. Si solo se tiene una medida por individuo, esto no es viable para generar un modelo aleatorio.
Cuando se tiene en cuenta el error de cada individio, permite ajustar el modelo para realizar un modelo de tipo mixto.
RESOLUCIÓN DE LAS PREGUNTAS
#1. ¿Por qué utiliza modelo mixto?
- Se utiliza un modelo mixto en este caso debido a que los datos presentan una estructura jerárquica o de medidas repetidas, donde se tienen múltiples observaciones para cada individuo (planta en este ejemplo). Al emplear un modelo mixto, se puede tener en cuenta tanto los efectos fijos (variables explicativas a nivel de grupo) como los efectos aleatorios (variabilidad de los individuos) en el análisis de los datos. Respecto al ejemplo del script, se analizan los datos de captación de CO2 en plantas de Quebec y Mississippi en diferentes tratamientos y concentraciones de CO2, se emplea un modelo mixto porque se tienen multiples mediciones para cada planta, al considerar el efecto aleatorio de la planta, se puede capturar la variabilidad entre las plantas y tener en cuenta las diferencias individuales en la captacion de CO2. Esto resulta ser bastante importante, ya que cada planta puede tener su propio nivel base de capatación de CO2 debido a diferentes factores como la genética o el ambiente. El uso de un modelo mixto tambien permite controlar la correlacion entre las observaciones dentro de cada planta, ya que se asume que las mediciones realizadas en la misma planta estan correlacionadas, esto es particularmente relevante cuando se trabaja con datos longitudinales o con medidas repetidas en el tiempo.
#2. ¿Qué ventajas tiene este tipo de modelos en agricultura?
- En muchos estudios agricolas, los diseños experimentales pueden ser desbalanceados debido a diversas razones, como la pérdida de datos, mediciones faltantes o diferentes tamaños de muestra en diferentes unidades experimentales. Los modelos mixtos pueden trabajar de manera eficiente frente a situaciones con datos desbalanceados, evitando la exclusión de observaciones valiosas y maximizando la utilización de los datos disponibles. Tambien se puede tener en cuenta que en la agricultura, es comun tener datos recopilados de manera jerárquica, donde se realizan mediciones repetidas en diferentes unidades experimentales, como parcelas, lotes o campos, en este caso los modelos mixtos permiten tener en cuenta esta estructura jerárquica al incorporar efectos aleatorios para obtener la variabilidad entre las unidades experimentales y las mediciones repetidas en el tiempo. Los modelos mixtos son flexibles y pueden utilizarse para analizar tanto variables de respuesta continua como categóricas. En la agricultura, esto puede ser especialmente útil, ya que se pueden modelar diversos resultados, como rendimiento de cultivos, contenido de nitrientes, presencia de enfermedades, etc.