Modelos mixtos
Carlos S
2023-06-25
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, unpacklibrary(ggplot2)
data("CO2")
uptake2 = CO2$uptake + runif(84, 0.02, 0.04)
view(uptake2)
ggplot(CO2, aes(x = conc, y = uptake2, color = Type)) + geom_point()
• Captación de CO2 de distintas plantas tiene alta variación. Se observa
plantas de Quebecen la parte superior y Missisipi en la inferior. •
Tratamiento experimental = enfriamiento
ggplot(CO2, aes(x = conc, y = uptake2, color = Type)) + geom_point(aes(shape = Treatment)) + theme_bw()
Se tiene una variable adicional = Planta (individuo), variable aleatoria, no predecible.
ggplot(CO2, aes(x = conc, y = uptake2)) +
geom_point(aes(color = Treatment))+
facet_wrap(~Type)
• Se tiene variable adicional = Planta (individuo), variable aleatoria,
no predecible.
ggplot(CO2, aes(x = conc, y = uptake2, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment)) +
theme_bw()
ggplot(CO2, aes(x = conc, y = uptake2, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment)) +
theme_bw() + facet_wrap(~Type)
• Cada planta segregante en cada tratamiento (enfriada y sin enfriar) =
variación relacionada con la planta. • Las pendientes son las mismas,
pero las intersecciones son diferentes. Análisis descriptivo. • Se ha
observado que las plantas en Quebec capturan o consumen concentraciones
exponencialmente más altas de CO2 a medida que aumentan las
concentraciones ambientales/en el aire. • Además, tanto en Quebec como
en Mississippi, la respuesta de las plantas no refrigeradas es buena. •
Cada planta/planta se comporta de manera diferente.
Fit1 <- lm(uptake2 ~ I(log(conc)) + Type:Treatment, data = CO2)
summary(Fit1)##
## Call:
## lm(formula = uptake2 ~ I(log(conc)) + Type:Treatment, data = CO2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.7140 -2.9014 0.5783 2.7617 8.8677
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -33.5265 4.0773 -8.223 3.21e-12 ***
## I(log(conc)) 8.4842 0.6784 12.506 < 2e-16 ***
## TypeQuebec:Treatmentnonchilled 19.5172 1.4404 13.549 < 2e-16 ***
## TypeMississippi:Treatmentnonchilled 10.1362 1.4404 7.037 6.35e-10 ***
## TypeQuebec:Treatmentchilled 15.9381 1.4404 11.065 < 2e-16 ***
## TypeMississippi:Treatmentchilled NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.668 on 79 degrees of freedom
## Multiple R-squared: 0.8227, Adjusted R-squared: 0.8137
## F-statistic: 91.64 on 4 and 79 DF, p-value: < 2.2e-16• Consumo de CO2 - Variable de respuesta a la concentración de CO2: Existe una relación entre las variables relacionadas con el consumo de CO2 de la planta que son importantes en el modelo. • Con base en la pendiente del intercepto, se rechaza la hipótesis nula (no igual a 0). • R-cuadrado ajustado = 81,3 % de los datos están relacionados (más del 70 % y los datos son positivos) • Hay sigularidad en los datso TypeMississippi:Treatmentchilled
ggplot(CO2, aes(x = conc, y = uptake2, color = Type)) +
geom_point() +
geom_smooth(method='lm', formula=y~x, se=FALSE, col='darkblue') +
theme_light()
#Normalidad de Residuos
shapiro.test(Fit1$residuals)##
## Shapiro-Wilk normality test
##
## data: Fit1$residuals
## W = 0.97699, p-value = 0.1374• El valor p es mayor que 5, los datos cumplen con el supuesto de normalidad • En este caso, el factor aleatorio son las plantas, por lo que se utiliza un modelo lineal mixto. • Con el paquete lme4, el modelo se puede exportar a plantas que no forman parte del experimento, por lo que cada planta tiene un punto de intersección diferente.
#Agregar el factor aleatorio: (Planta) - Para estimar los parámetros del modelos
Fit2 <- lmer(uptake2 ~ I(log(conc)) + Type:Treatment + (1|Plant), data = CO2)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientsummary(Fit2)## Linear mixed model fit by REML ['lmerMod']
## Formula: uptake2 ~ I(log(conc)) + Type:Treatment + (1 | Plant)
## Data: CO2
##
## REML criterion at convergence: 482.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6927 -0.5397 0.1038 0.6914 1.8913
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant (Intercept) 2.156 1.468
## Residual 20.258 4.501
## Number of obs: 84, groups: Plant, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -33.5265 4.0220 -8.336
## I(log(conc)) 8.4842 0.6542 12.969
## TypeQuebec:Treatmentnonchilled 19.5172 1.8349 10.637
## TypeMississippi:Treatmentnonchilled 10.1362 1.8349 5.524
## TypeQuebec:Treatmentchilled 15.9381 1.8349 8.686
##
## Correlation of Fixed Effects:
## (Intr) I(l()) TypQbc:Trtmntn TypM:T
## I(log(cnc)) -0.947
## TypQbc:Trtmntn -0.228 0.000
## TypMssssp:T -0.228 0.000 0.500
## TypQbc:Trtmntc -0.228 0.000 0.500 0.500
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientLa media de los residuos cambia en comparación con el modelo lineal simple.
# Valores efectos variables:
ranef(Fit2)## $Plant
## (Intercept)
## Qn1 -0.89954179
## Qn2 -0.07459541
## Qn3 0.97413721
## Qc1 -0.76063525
## Qc3 0.35510527
## Qc2 0.40552998
## Mn3 -0.78421160
## Mn2 0.59197757
## Mn1 0.19223404
## Mc2 -1.56718838
## Mc3 0.63390689
## Mc1 0.93328150
##
## with conditional variances for "Plant"# Valores efectos fijos:
fixef(Fit2)## (Intercept) I(log(conc))
## -33.526518 8.484187
## TypeQuebec:Treatmentnonchilled TypeMississippi:Treatmentnonchilled
## 19.517188 10.136157
## TypeQuebec:Treatmentchilled
## 15.938135# 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.823 0.814 4.67 91.6 7.07e-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.03e-20
## 3 TypeQuebec:Treatmentnonchilled 19.5 1.44 13.5 2.64e-22
## 4 TypeMississippi:Treatmentnonchilled 10.1 1.44 7.04 6.35e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 1.00e-17
## 6 TypeMississippi:Treatmentchilled NA NA NA NA• Tiddy nos muestra cual va a ser el estimador El modelo lineal mixto no se tienen R cuadrado:
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 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 NA• vemos las variables aleatorias con sus estimaciones para cada observación; divide en efectos fijos y aleatorios. El tratamiento de enfriamiento y el origen de Mississippi tienen menos efecto sobre la absorción de CO2. Hay diferencias entre las plantas individuales. Puede ver los valores muy bien estimados para las interacciones entre las plantas y el enfriamiento, siendo la estimación más alta TypeQuebec:Treatmentnonchilled y luego TypeQuebec:Treatmentchilled. Esto le permite confirmar el mayor secuestro de CO2 encontrado en los dos tratamientos de la planta de Quebec a través del análisis gráfico. ### RESUMEN: • Aumenta la precisión de la investigación de cultivos, porque tiene en cuenta la variabilidad de los individuos, permite mediciones repetidas en el tiempo; considera variables fijas y aleatorias (dentro y fuera de un conjunto o bloque). • Puede ser útil principalmente en programas de mejoramiento genético y condiciones óptimas de diferentes cultivos a través de diferentes períodos y/o etapas fenológicas.
EFECTOS FIJOS GLMER.
data("ChickWeight")
#Total pollos
ChickWeight$Chick %>% length()## [1] 578#Tipos de pollos
ChickWeight$Chick %>% unique ()%>% length## [1] 50head(ChickWeight)## Grouped Data: weight ~ Time | Chick
## weight Time Chick Diet
## 1 42 0 1 1
## 2 51 2 1 1
## 3 59 4 1 1
## 4 64 6 1 1
## 5 76 8 1 1
## 6 93 10 1 1NOTA: Como se tienen valores positivos en el peso y no se tiene decimales, por tal motivo, se podría analizar como un gml tipo poisson.
ggplot(ChickWeight, aes(x = Time, y = weight)) +
geom_point(aes(color = Diet)) +
geom_path(aes(color = Diet, group = Chick))
Hay pollos donde a pesar de que se tienen dietas simialres tienen comportamientos distintos. Esto indica que si se tiene estas medidades repetitas se puede hacer un glmer.
#gml tipo poisson
Fit1_Poisson <- glm(weight ~ Diet:Time, data = ChickWeight, family = poisson)
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 NASe observan diferencias entre los modelos, porque, en el modelo lineal da un estimate del intercepto de 3.85 y cuando se aleatorizan la variable individuo en el modelo mixtose observaun ligero cambio de ese intercepto a 3.83. Por otro lado, se observa que la identidad del pollo capta alguna variabilidad y esto genera que no sean iguales los modelos. Para concluir si se toma en cuenta la variablidad entre individuos esto mejora el modelo, asimismo, cuando se toma la indentidad del pollo es necesario que se tenga más de una medida por cada individuo para poder hacer el modelo aleatorio.
CONCLUSIÓN: Se usa el modelo lineal mixto porque se observa que la identidad de los individuos genera un cambio en los resultados, por esta razón, este modelo permite dividir entre variables fijas (controlados) y aleatorias (no controlados). En cambio, el modelo lineal solo toma en cuenta un efecto.
Importancia del Modelo lineal mixto en la agricultura:
El script utiliza un modelo lineal mixto porque considera la presencia de efectos aleatorios, en este caso, el efecto de las plantas individuales en la captación de CO2. El modelo mixto es adecuado cuando se tienen variables aleatorias o efectos no predecibles que influyen en la respuesta. En agricultura, los modelos lineales mixtos ofrecen varias ventajas importantes: 1. Consideración de la variabilidad individual: Los modelos mixtos tienen la capacidad de tener en cuenta la variabilidad individual que puede existir entre las plantas, animales o parcelas de una misma especie. Esto es especialmente relevante en estudios agrícolas, ya que las características de las plantas individuales pueden variar significativamente y afectar los resultados de las mediciones. 2. Medidas repetidas: Los modelos mixtos son útiles cuando se realizan mediciones repetidas en el tiempo o en diferentes condiciones. Permiten modelar la estructura de correlación entre las observaciones repetidas, teniendo en cuenta la dependencia entre ellas. En agricultura, esto es común en estudios de seguimiento de cultivos a lo largo del tiempo o en ensayos que involucran múltiples parcelas o tratamientos. 3. Manejo de efectos fijos y aleatorios: Los modelos mixtos pueden incorporar tanto efectos fijos como aleatorios. Los efectos fijos representan factores que se consideran constantes y de interés en el estudio, como los tratamientos experimentales o las variables ambientales. Los efectos aleatorios capturan la variabilidad no predecible asociada a las unidades de observación, como las plantas individuales. Al considerar ambos tipos de efectos, se obtiene una estimación más precisa de los efectos fijos y se tiene en cuenta la estructura de dependencia entre las observaciones. 4. Flexibilidad en el diseño experimental: Los modelos mixtos permiten manejar diseños experimentales complejos, como diseños completamente aleatorizados, bloques aleatorizados, parcelas divididas y diseños de medidas repetidas. Esto es relevante en agricultura, donde los estudios a menudo involucran diseños experimentales específicos para controlar fuentes de variabilidad y obtener conclusiones más sólidas.