MODELOS MIXTOS
Dario_Gonzalez
2023-06-25
#install.packages("tidyverse")
#install.packages("broom.mixed")
#install.packages("lme4")
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 --
## v dplyr 1.1.2 v readr 2.1.4
## v forcats 1.0.0 v stringr 1.5.0
## v ggplot2 3.4.2 v tibble 3.2.1
## v lubridate 1.9.2 v tidyr 1.3.0
## v purrr 1.0.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## i 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#Datos co2
data("CO2")
uptake2 = CO2$uptake + runif(84,0.02,0.04)#Analisis descriptivo
ggplot(CO2, aes(x = conc, y = uptake2, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment)) +
theme_bw()
De acuerdo al grafico, podemos decir que existe un bajo consumo y una baja variacion a una concentracion baja. Se observa que las plantas de Quebec captan mayor co2 que las de Misisipi, aunque en misisipi se observa una mayor diferencia en los tratamientos de enfriado y no enfriado.
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 / coefficientEn este caso, se usa el “lmer” para poder utilizar este modelo en plantas que no estan dentro de este experimento (variable aleatoria), en donde cada una tendra un intercepto diferente.
#Modelo lineal simple
broom::glance(Fit1)## # A tibble: 1 x 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.00e-29 4 -246. 504. 519.
## # i 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>broom::tidy(Fit1)## # A tibble: 6 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -33.5 4.08 -8.22 3.19e-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.30e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 9.91e-18
## 6 TypeMississippi:Treatmentchilled NA NA NA NACon “Tidy” obtenemos el estimado de cada factor.
#MOdelo lineal mixto
broom.mixed::glance(Fit2)## # A tibble: 1 x 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 x 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.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()*En “tidy” podemos ver los valores estimados para la interacción entre las plantas y el enfriamento. la de mayor estimacion es en el tratamiento no enfriado en Quebec seguida por el tratamiento enfiado en Quebec, confirmando lo interpretado en el analisis gráfico en donde con cualquiera de los dos tratamientos las plantas de Quebec tienen mayor consumo de CO2.
#GLMER Con esta función se agrega la familia poisson porque no se tienen datos decimales y son valores positivos. La base de datos “peso de pollos” es muy similar a la anterior, se tiene el peso, tiempo, dieta y la identidad del pollo.
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 los pollos tienen dietas similares, pero
van a tener comportamientos distinto aun partiendo de un mismo
punto.
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 x 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 x 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 NAcon la funcion “tidy” nos muestra un aumento de peso en el tiempo y una interración 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.
*Para poder realizar este tipo de modelo es necesario que se tenga más de una medida por cada individuo
Conclusión
Este modelo es muy importante en la agricultura, como tambien para otros campos, porque nos permite llevar a campo experimentos que son evaluados con factores requeridos por un cultivo, teniendo en cuenta la variable aleatoria