Modelos mixtos
Saydi Torres
2023-06-26
# Librerias
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.3.1## Warning: package 'readr' was built under R version 4.3.1## Warning: package 'forcats' was built under R version 4.3.1## Warning: package 'lubridate' was built under R version 4.3.1## ── 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)## Warning: package 'broom.mixed' was built under R version 4.3.1library(lme4)## Warning: package 'lme4' was built under R version 4.3.1## Loading required package: Matrix## Warning: package 'Matrix' was built under R version 4.3.1##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack1: MODELOS LINEALES MIXTOS
Tratamiento experimental, algunas plantas están refrigeradas y otras no refrigeradas con el fin de determinar cuanto carbono captan.
data("CO2")ggplot(CO2,
aes(x = conc, y = uptake, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment))
-A medida que la concentración de CO2 aumenta, se puede notar una
considerable variación. -Cada planta muestra un comportamiento único, y
existe variabilidad tanto en las plantas refrigeradas como en las no
refrigeradas. -Hay una variable asociada a la variabilidad que es
predecible grafiamnete -Existe una variable de tipo aleatoria.
# Para modelo lineal simple
fit1 <- lm(uptake ~ conc + Type + Treatment, data = CO2)
# Para agregar el Factor aleatorio
fit2 <- lmer(uptake ~ conc + Type + Treatment + (1 | Plant), data = CO2)# Para modelo lineal Simple
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.89e-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 2.08e-31
## 2 conc 0.0177 0.00230 7.72 2.87e-11
## 3 TypeMississippi -12.7 1.35 -9.37 1.67e-14
## 4 Treatmentchilled -6.86 1.35 -5.08 2.46e- 6# Para modelo lineal MIXTO
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 NA-Cada una de las plantas tendrá un intercepto distinto. -Se observa que nuestra variable es aleatoria, por lo que cada una de las plantas va a tener un intercepto distinto
Fit1 <- lm(uptake ~ I(log(conc)) + Type:Treatment, data = CO2)
Fit2 <- lmer(uptake ~ I(log(conc)) + Type:Treatment + (1|Plant), data = CO2)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientbroom::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.01e-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.6 4.08 -8.23 3.09e-12
## 2 I(log(conc)) 8.48 0.678 12.5 2.03e-20
## 3 TypeQuebec:Treatmentnonchilled 19.5 1.44 13.6 2.61e-22
## 4 TypeMississippi:Treatmentnonchilled 10.1 1.44 7.04 6.29e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 9.94e-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.6 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 NAGLM para efectos fijos
Se implementan dos factores Tiempo y Dieta
data("ChickWeight")# Vista Gráfica
ggplot(ChickWeight,
aes(x = Time, y = weight)) +
geom_point(aes(color = Diet)) +
geom_path(aes(color = Diet, group = Chick))
-Se observa que los pollos con la misma dieta presentan una gran
variabilidad en sus pesos. -Se puede inferir que a medida que pasa el
tiempo, los pollos tienden a ganar peso, y este aumento puede depender
del tipo de dieta aplicada. Además, se sugiere que la dieta 1 muestra
una menor efectividad en comparación con las otras dietas. -Debido a que
los valores de peso son siempre enteros positivos, existe una indicación
de que un modelo GLM de tipo Poisson, que tiene en cuenta la interacción
entre la dieta y el tiempo, puede ser apropiado para analizar estos
datos.
# El peso es el resultado de la interacción entre el tiempo y la dieta
fit1_poisson <- glm(weight ~ Diet:Time, data = ChickWeight, family = poisson)
#Agregar una familia
# Para agregar el Factor aleatorio Chick
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::glance(fit1_poisson)## # A tibble: 1 × 8
## null.deviance df.null logLik AIC BIC deviance df.residual nobs
## <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
## 1 22444. 577 -3890. 7791. 7812. 4038. 573 578broom::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::glance(fit2_poisson)## # A tibble: 1 × 7
## nobs sigma logLik AIC BIC deviance df.residual
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 578 1 -2645. 5302. 5328. 1294. 572broom.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 NA-A medida que transcurre el tiempo, se observa un incremento en el peso de los pollos. Además, se evidencia que la interacción entre el tiempo y la dieta es distinta para cada tipo de dieta. Estos hallazgos respaldan lo mencionado en el análisis descriptivo previo, confirmando que el peso está influenciado tanto por el factor tiempo como por la elección de la dieta.
-Se observa que la dieta 3 es la más efectiva en términos de aumento de peso en los pollos, ya que muestra una pendiente más pronunciada (estimada en 0.086) en comparación con las dietas 1, 2 y 4. -Hay variabilidad en los valores de intercepto al considerar la identidad de los pollos, lo que indica diferencias individuales en los pesos iniciales. -El aumento de peso en los pollos es exponencial a medida que pasa el tiempo, lo que implica que ganan más peso a medida que avanza el tiempo. Sin embargo, este aumento no es uniforme en todas las dietas utilizadas, mostrando diferencias en el ritmo de crecimiento.
2: Ruido para el modelo lineal mixto
Tratamiento experimental, algunas plantas están refrigeradas y otras no refrigeradas con el fin de determinar cuanto carbono captan.
ruido <- runif(84, 0.02, 0.04)
ruido## [1] 0.03678335 0.03340179 0.03607327 0.02011161 0.02521773 0.02861101
## [7] 0.02111284 0.03458931 0.02014120 0.03202492 0.02270467 0.02394491
## [13] 0.02325352 0.03025648 0.02037693 0.03860263 0.03100239 0.02460322
## [19] 0.03665638 0.03139130 0.03921750 0.02633830 0.02115290 0.03241917
## [25] 0.02407208 0.03463015 0.02439820 0.03817806 0.02675522 0.03975435
## [31] 0.03601620 0.03640423 0.02650162 0.02934394 0.03897278 0.03794943
## [37] 0.03754709 0.03694537 0.02665275 0.03661110 0.03989632 0.02785014
## [43] 0.02172464 0.02203790 0.02626405 0.02906599 0.03623746 0.02036326
## [49] 0.02692563 0.03405705 0.03393991 0.03514571 0.02814602 0.03432223
## [55] 0.03527813 0.03020741 0.02044152 0.02260959 0.02945558 0.02026133
## [61] 0.02222300 0.03833648 0.02819002 0.03174305 0.03099342 0.02385891
## [67] 0.02061509 0.02043090 0.02902841 0.03209216 0.03373249 0.03023525
## [73] 0.02961134 0.02058017 0.03841821 0.02037188 0.03867773 0.02332559
## [79] 0.03599576 0.03796744 0.02949132 0.02229967 0.03781590 0.02831596uptake2 <- CO2$uptake + ruido
(uptake2)## [1] 16.036783 30.433402 34.836073 37.220112 35.325218 39.228611 39.721113
## [8] 13.634589 27.320141 37.132025 41.822705 40.623945 41.423254 44.330256
## [15] 16.220377 32.438603 40.331002 42.124603 42.936656 43.931391 45.539217
## [22] 14.226338 24.121153 30.332419 34.624072 32.534630 35.424398 38.738178
## [29] 9.326755 27.339754 35.036016 38.836404 38.626502 37.529344 42.438973
## [36] 15.137949 21.037547 38.136945 34.026653 38.936611 39.639896 41.427850
## [43] 10.621725 19.222038 26.226264 30.029066 30.936237 32.420363 35.526926
## [50] 12.034057 22.033940 30.635146 31.828146 32.434322 31.135278 31.530207
## [57] 11.320442 19.422610 25.829456 27.920261 28.522223 28.138336 27.828190
## [64] 10.531743 14.930993 18.123859 18.920615 19.520431 22.229028 21.932092
## [71] 7.733732 11.430235 12.329611 13.020580 12.538418 13.720372 14.438678
## [78] 10.623326 18.035996 17.937967 17.929491 17.922300 18.937816 19.928316# Vista Gráfica
ggplot(CO2,
aes(x = conc, y = uptake2, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment))
-El análisis revela que la adición de ruido en los datos produce cambios
notables en la distribución de la captación de CO2 por parte de las
plantas. A pesar de esto, persiste la evidencia de que las plantas
originarias de Quebec tienen una mayor capacidad de captación en
comparación con las plantas de Mississippi. -Se observa claramente que
el tratamiento utilizado desempeña un papel importante en la capacidad
de captación, ya que en la mayoría de los casos, las plantas sometidas a
refrigeración presentaron una menor capacidad de captación. -La gráfica
también muestra que las diferencias en la captación de CO2 no son tan
notorias al agregar ruido a los datos en comparación con los resultados
anteriores.
-Se observa una alta variabilidad asociada a cada una de las plantas utilizadas, lo que indica que es necesario utilizar un modelo lineal mixto y aleatorizar la variable planta (Plant) para evitar la generación de efectos no deseados en la variable respuesta.
Fit3 <- lm(uptake2 ~ conc + Type + Treatment, data = CO2)
Fit4 <- lmer(uptake2 ~ conc + Type + Treatment + (1 | Plant), data = CO2)# MODELO LINEAL MIXTO
broom::glance(Fit3)## # 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.90e-20 3 -270. 551. 563.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>broom::tidy(Fit3)## # 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.96e-31
## 2 conc 0.0177 0.00230 7.72 2.87e-11
## 3 TypeMississippi -12.7 1.35 -9.37 1.66e-14
## 4 Treatmentchilled -6.86 1.35 -5.07 2.48e- 6broom.mixed::glance(Fit4)## # 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(Fit4)## # 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 NAFit3 <- lm(uptake2 ~ I(log(conc)) + Type:Treatment, data = CO2)
Fit4 <- lmer(uptake2 ~ I(log(conc)) + Type:Treatment + (1|Plant), data = CO2)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientbroom.mixed::glance(Fit3)## # 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.04e-29 4 -246. 504. 519.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>broom::tidy(Fit3)## # 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.20e-12
## 2 I(log(conc)) 8.48 0.678 12.5 2.03e-20
## 3 TypeQuebec:Treatmentnonchilled 19.5 1.44 13.6 2.63e-22
## 4 TypeMississippi:Treatmentnonchilled 10.1 1.44 7.04 6.34e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 9.92e-18
## 6 TypeMississippi:Treatmentchilled NA NA NA NAbroom.mixed::glance(Fit4)## # 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(Fit4)## # 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.84 10.6
## 4 fixed <NA> TypeMississippi:Treatmentnonch… 10.1 1.84 5.52
## 5 fixed <NA> TypeQuebec:Treatmentchilled 15.9 1.84 8.68
## 6 ran_pars Plant sd__(Intercept) 1.47 NA NA
## 7 ran_pars Residual sd__Observation 4.50 NA NAConclusiones:
-Al analizar los resultados obtenidos, se puede observar que la variable individuo (Planta) fue interceptada para evitar alterar los resultados de la variable respuesta. Esta interceptación se refleja en los valores estimados de los coeficientes en la columna “estimate”, los cuales son utilizados para aplicar el modelo en contextos diferentes al experimento actual, ya que no están influenciados por factores específicos como la planta utilizada en el experimento.
-En cuanto a los tratamientos, se observa que las plantas sometidas a refrigeración presentan una menor captación de carbono en comparación con aquellas que no fueron refrigeradas. Además, se confirma que las plantas de Mississipi, en general, muestran una menor captación de carbono en comparación con las plantas de Quebec.
-Es importante destacar que al interceptar el factor planta y aleatorizarlo, se evita considerar el coeficiente de variación que podría haber causado efectos no deseados en la variable respuesta.
-Los resultados obtenidos en la tabla también respaldan las conclusiones del análisis gráfico. Estos resultados explican por qué las plantas de Mississipi y aquellas sometidas al tratamiento de refrigeración presentan una menor capacidad de captación de carbono.
Conclusiones finales
-Los modelos de efectos mixtos son herramientas estadísticas valiosas en el análisis de datos, ya que permiten controlar variables irrelevantes y evitar que afecten la variable respuesta. Al considerar la variabilidad entre individuos y agregar un factor aleatorio en el modelo lineal mixto, es necesario obtener múltiples medidas para cada individuo evaluado. Esto garantiza un análisis más preciso y una mejor comprensión de los factores que influyen en los resultados en el campo de la agronomía.
-En la agricultura, donde múltiples variables pueden afectar el rendimiento de los cultivos, los modelos mixtos son especialmente útiles. Permiten separar las variaciones causadas por variables no deseadas, capturando así la variación no explicada y estableciendo una relación adecuada entre la variable respuesta y los tratamientos de interés. Al diferenciar las variables relevantes y obtener resultados específicos para cada variable de interés, los modelos lineales mixtos facilitan la toma de decisiones informadas en el ámbito agronómico.