| title: “Entrega final modelos mixtos” |
| author: “Luisa Perea” |
| date: “2023-06-25” |
| output: html_document |
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.3.1## Warning: package 'ggplot2' 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) ### Resumenes del tipo del modelo## 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, unpacklibrary(ggplot2)Base de Datos
Modelo lineal mixto Proporciona la felxibilidad necesario para modelar nos solo las medias no las varianzas y covarianzas de los datos.
es un experimentos con plantas, con diferentes dosis de carbono y cuando CO2 se va a captar esa es la variable respuesta
a diferentes concentraciones, cual capta mayor Co2 se usa ruido para tener una variablidad en los datos, ademas de que se pondra distribucion unifrome.
set.seed(123)
data("CO2")
uptake2 = CO2$uptake + runif (84, 0.02, 0.04)
CO2$uptake2 = uptake2###Analisis Descriptivo la variable planta se supone son escogidas de manera aleatoria, pero esto no influye sobre la capatacion de CO2 es una variable aletoria. Se va a usar el modelo lm4, pero primero se usara el modelo lineal
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()
Modelo linel simple Se usar uptake+ concentracion+ tratamiento = lm lmer
= factor aleatorio, esto va ahcer que cada una de las plantas tenga un
intercepto diferente.
Fit1 <- lm(uptake2 ~ conc + Type + Treatment, data = CO2 )Fit2 <- lmer(uptake2 ~ conc + Type + Treatment +(1 | Plant), data = CO2 )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- 6El broom= en el tidy cada uno de los factores cual va a hacer el estiamdor
broom.mixed::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>Modelo mixto Se usa lo siguiente: No se tendra R2 sino otros valores, y valores extra
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 NASe va apoder ver las variables aleatorias para cada planta e intercepto para efectos fijos y aleatorios
broom.mixed::tidy(Fit2) %>% View()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 / coefficientConclusion: tipo con el tratamiento tiene una interaccion, hay diferencias entre los tartamientos segun el ggplot.
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(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.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()broom.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 NASe pueden agregar familiar con lmer4: se escoge otra base de datos de peso de pollos, se va a tener tambien la identidada d eun pollo, ademas de peso, tiempo y dieta, como se muestra a continuacion:
data("ChickWeight")
View(ChickWeight)
ChickWeight$Chick%>% length()## [1] 578ChickWeight$Chick %>% unique()%>% length()## [1] 50Es un gl poisson, el peso es explicado por la Dieta en interaccion con el tiempo. La base de datos es poisson porque el peso siempre es positivo.
Fit1_Poisson <- glm(weight ~ Diet:Time, data = ChickWeight, family = poisson)El aumento del peso en el tiempo, cada dieta tendra diferentes pendientes
ggplot(ChickWeight, aes(x = Time, y = weight))+ geom_point(aes(color = Diet)) + geom_path(aes(color = Diet, group = Chick))
Interpretacion ggplot: se va a tener que el grupo se va a dar por el
pollo, cada pollo a pesar de tener dietas similares desde el punto de
partida, vana tener comportamientos distintos.
Fit2_Poisson <- glm(weight ~ Diet:Time, data = ChickWeight, family = poisson)Al tener medidas repetidas se puede hacer un Fit2, se va a hacer un glmer + los errores aleatorios por individuo
Fit1_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: 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 NAbroom.mixed::tidy(Fit2_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. 0En el modelo glm normal se tiene un intercepto de 3.85 y cuando se consideran los individuos este estimado o intercepto baja a 3.83, esto quiere decir que la identidad de los pollos capta alguna variabilidad, lo cual cauda que no sean iguales.
View(broom.mixed::tidy(Fit2_Poisson))Se observa que la dieta con la pendiente mas alta es la dieta 3 con una ganacia de 0.08624478g/ dia
View(broom::tidy(Fit1_Poisson))Se observa que
library(sjPlot)## Warning: package 'sjPlot' was built under R version 4.3.1## #refugeeswelcomesjPlot::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
Nota: El tomar en cuenta la variabilidad de un individuo vs otro va a agregar confiabilidad, cuando se toma en cuenta la identidad se debe tener en cuenta mas de una medida por un individuo, no solo es util para tener un modelo aleatorio.
En la agricultura se puede utilizar un modelo mixto cuando se estan haciendo trabajos de fitomejoramiento, en donde se cominan factores fijos, que son aquellos que el investigador puede controlar, como tambien factores aletorios que son los que estan fuera del control.
Los modelos mixtos basados en prediccion, se han usado en plantas como soya, maiz, caña en donde el investigador tiene como objetivo seleccionar los individuos y familas que presenten las caracteristicas deseadas para ser recomendadas en futuras investigaciones. Los efectos aleatorios en este caso serian la capacidad especifica de combinacion y el efecto de los bloques
REFERENCIAS Bandera-Fernández, Evelyn, & Pérez-Pelea, Leneidy. (2018). Los modelos lineales generalizados mixtos. Su aplicación en el mejoramiento de plantas. Cultivos Tropicales, 39(1), 127-133. Recuperado en 26 de junio de 2023, de http://scielo.sld.cu/scielo.php?script=sci_arttext&pid=S0258-59362018000100019&lng=es&tlng=es.