Setup.

El modelo lineal mixto se utiliza para hacer que las medias, varianzas y covarianzas sean más “flexibles”.

#install.packages

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 errors
library(broom.mixed)
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

Primer modelo.

data("CO2")

uptake2 = CO2$uptake + runif(84, 0.04, 0.09)
uptake2
##  [1] 16.057276 30.451652 34.859843 37.282537 35.384733 39.265259 39.759311
##  [8] 13.641123 27.380591 37.161402 41.866319 40.688173 41.485973 44.378798
## [15] 16.287411 32.489718 40.385843 42.169039 42.962640 43.979473 45.551488
## [22] 14.276920 24.185653 30.370135 34.668188 32.572756 35.460099 38.771378
## [29]  9.361548 27.343475 35.065480 38.844475 38.672115 37.575503 42.487042
## [36] 15.156786 21.071451 38.159608 34.070736 38.986360 39.653808 41.456118
## [43] 10.680049 19.274765 26.286861 30.070364 30.963291 32.448584 35.560308
## [50] 12.056549 22.078177 30.641623 31.844167 32.488167 31.180225 31.540280
## [57] 11.350736 19.464444 25.886039 27.960491 28.552228 28.167547 27.884246
## [64] 10.563452 14.955251 18.188820 18.945564 19.579244 22.260716 21.982617
## [71]  7.785831 11.448243 12.351045 13.055331 12.569392 13.785241 14.442602
## [78] 10.653848 18.055111 17.955749 17.969659 17.979829 18.968227 19.944886

Análisis Descriptivo.

ggplot(CO2, aes(x = conc, y = uptake2)) + geom_point()

data("CO2")

ggplot(CO2, aes(x = conc, y = uptake, color = Type)) + geom_point()

  • Si medimos cual es la captaciòn de dioxido de carbono de distintas plantas, en distintas concentraciones, podremos observar, en este caso, una amplia variacion, es decir, al principio hay menos, pero al final hay mayor variabilidad.

  • En general, de la grafica se puede concluir que las plantas de Quebes tienen un mayor consumo en comparacion con las de Mississippi.

data("CO2")

ggplot(CO2, aes(x = conc, y = uptake, color = Type)) + geom_point(aes(shape = Treatment)) + theme_bw()

Variables:

  • Plantas (variable aleatoria).

  • Tratamiento.

  • Concentracion.

  • Consumo.

data("CO2")

ggplot(CO2, aes(x = conc, y = uptake, color = Type)) + 
  geom_point(aes(shape = Treatment)) +
  geom_path(aes(group = Plant, lty = Treatment)) +
  theme_bw()

  • Si vemos a cada planta aislada, podemos apreciar comportamientos o tendencias inesperadas, si bien las plantas se escogieron de manera aleatoria, son del mismo tipo y por eso, esperariamos un comportamiento mas similar entre las repeticiones, por lo que podemos intuir que este efecto es causado por alguna variable que determinaria el resultado espeficico de cada una.

  • Para lidiar con esos modelos aleatorios, se utiliza el siguiente paquete.

# Lo que usariamos originalmente si tuvieramos un modelo lineal simple.

Fit1 <- lm(uptake ~ I(log(conc)) + Type:Treatment, data = CO2)

# Lo que usamos dadas las circunstancias en las que estamos con el factor aleatorio.

Fit2 <- lmer(uptake ~ I(log(conc)) + Type:Treatment + (1 | Plant), data = CO2)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

  • Se coloca la expresión “I(log())”, puesto a que la gráfica de los resultados toma una función de tipo logaritmica.

  • Como se notan diferencias entre el tipo y los tratamientos, se utiliza la expresión “Type:Treatment”, en lugar de sumarlos.

  • Lo que hace esta modificacion en particular es que cada una de las plantas va a tener un intercepto disitnto.

Para un modelo lineal.

# Estimadores.

broom::glance(Fit1)
broom::tidy(Fit1)

Para un modelo lineal mixto.

broom.mixed::glance(Fit2)
broom.mixed::tidy(Fit2)
broom.mixed::tidy(Fit2)

Segundo Modelo.

data("ChickWeight")

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?
tidy(Fit1_Poisson)

Analisis descriptivo.

ggplot(ChickWeight, aes(x = Time, y = weight)) + geom_point(aes(color = Diet))

ggplot(ChickWeight, aes(x = Time, y = weight)) +
  geom_point(aes(color = Diet)) +
  geom_path(aes(color = Diet, group = Chick))

Para un modelo lineal.

broom::tidy(Fit1_Poisson)

Para un modelo lineal mixto.

broom.mixed::tidy(Fit2_Poisson)  
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.

Referencia.