Ultima entega - Diseño de experimentos
Paula Andrea Sanchez Camargo
2023-06-24
#Modelos lineales mixtos
library(tidyverse)## Warning: package 'tidyverse' 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
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack#CO2 mas el ruido en uptake
data("CO2")
uptake2 = CO2$uptake + runif(84,0.02,0.04)
view(uptake2)#Analisis descriptivo
ggplot(CO2, aes(x = conc, y = uptake2)) + geom_point()
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)) +
geom_point(aes(color = Treatment))+
facet_wrap(~Type)
* Mediante la grafica se puede ver que hay una variacion relativamente
baja a menores concentraciones, sin embargo cuando esta concentracion es
alta, la variavilidad tambien lo es, por otro lado tambien se evidencia
que las plantas de Quebec captaron mayor C02 en diferentes
concentraciones, si bien las plantas de Mississippi tienen un
comportamiento parecido, las de Quebec fueron capaces de captar mas CO2,
adicionalmente el tratamiento de nochilled tiene mayor efecto en la
captura de CO2 en las plantas provenientes de Quebec y Mississipi.
Para poder utilizar el modelo con plantas que no estén dentro de este experimento y adaptar esa variable aleatoria se hace uso de lme4
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 / coefficientCon lme4 el modelo puede ser aplicado a plantas que no fueron tenidas en cuenta para el experimento, por lo que cada planta tendra un intercepto distinto
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.09e-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.24e-12
## 2 I(log(conc)) 8.48 0.678 12.5 2.05e-20
## 3 TypeQuebec:Treatmentnonchilled 19.5 1.44 13.6 2.63e-22
## 4 TypeMississippi:Treatmentnonchilled 10.1 1.44 7.04 6.28e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 1.00e-17
## 6 TypeMississippi:Treatmentchilled NA NA NA NACon la función Broom::Tidy se optiene el estimado de los factores
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 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.33
## 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.53
## 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()- Mediante tidy se observan los valores estimados para las interacciones entre las plantas y el enfriamento, donde la que tiene una estaimación mayor es TypeQuebec:Treatmentnonchilled con 19.516, seguida por TypeQuebec:Treatmentchilled con 15.936, teniendo con esto la confirmacion de lo visto en el analisis descriptivo donde con cualquiera de los dos tratamientos las plantas de Quebec tienen mayor captura de CO2.
#GLMER * Con esta función se agrega una familia, en este caso de tipo poisson ya que no hay datos decimales y son valores positivos
- La base de datos ChickWeight (peso de pollos) es similar a la anterior, ya que proporciona peso, tiempo y dieta, la identidad del pollo los cuales son 50
data("ChickWeight")ggplot(ChickWeight, aes(x = Time, y = weight)) +
geom_point(aes(color = Diet)) +
geom_path(aes(color = Diet, group = Chick))
*En la gráfica se observa que aunque los pollos tengan dietas similares
van a tener comportamientos distintos aun habiendo partido de un mismo
punto, donde todos con el transcurso del tiempo aumentan de peso
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?*El peso es explicado por la dieta en interacción con el tiempo
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 NA*tidy nos muestra un aumento de peso en el tiempo y las interacciones de cada dieta con el tiempo donde cada dieta es diferente
*La dieta con un mayor aumento en el peso es la 3 con una ganancia diaria de 0,086 y seguida por la dieta 4 con 0,076 de ganancia diaria
*Para poder realizar este tipo de modelo es necesario tener más de una medida por cada individuo
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
Conclusión
Este tipo de modelos es de gran relevancia e importancia ya que permite analizar datos de mediciones temporales a individuos modelando tanto variables fijas como aleatorias, por ejemplo en la agricultura es clave ya que estas variables se pueden estimar ya sea el rendimiento o respuesta teniendo en cuenta el efecto aleatorio entre los cultivos y dentro de los mismos, pues durante el desarrollo de los estos hay muchas variables que inciden en el rendimiento de la planta entonces se pueden identificar las variables que son necesarias para comprender el efecto del tratamiento.