library(lme4)## Loading required package: Matrixlibrary(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(magrittr)PARTE 1 - VIDEO: MODELOS LINEALES MIXTOS
Plantas a distintas concentraciones de dióxido de carbono, en el que algunas plantas están refrigeradas y no refrigeradas (tratamiento experimental), con el fin de determinar cuanto carbono captan.
data("CO2") # Base de DatosAnálisis Descriptivo
ggplot (CO2, aes(x = conc, y = uptake, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment)) +
theme_bw()
* Se puede observar alta variación, conforme aumenta la concentración de CO2.
* Cada planta tiene un comportamiento diferente, donde siempre hay una variabilidad tanto en las plantas refrigeradas como las no refrigeradas.
* Variabilidad asociada a una variable (que no esperemos) y que es de manera predecible y por otro lado, una variable que no se querra utilizar para poder hacer predicciones a futuro.
* VARIABLE DE TIPO ALEATORIA.
Modelo Lineal Mixto
# Si nos enfrentaramos ante un modelo lineal simple:
Fit1 <- lm(uptake ~ conc + Type + Treatment, data = CO2)
# Para agregar un factor aleatorio
Fit2 <- lmer(uptake ~ 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.89e-20 3 -270. 551. 563.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>broom::tidy(Fit1) # Aplicado al modelo lineal simple## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 29.3 1.54 19.0 2.08e-31
## 2 conc 0.0177 0.00230 7.72 2.87e-11
## 3 TypeMississippi -12.7 1.35 -9.37 1.67e-14
## 4 Treatmentchilled -6.86 1.35 -5.08 2.46e- 6broom.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) #Aplicado al modelo lineal mixto## # 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 NA- Cada una de las plantas tendrá un intercepto distinto.
- Se observa que nuestra variable es aleatoria, por lo que se obtiene un estimador para cada planta y para cada observación.
Fit1 <- lm(uptake ~ I(log(conc)) + Type:Treatment, data = CO2)
Fit2 <- lmer(uptake ~ I(log(conc)) + Type:Treatment + (1|Plant), data = CO2)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientbroom::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.01e-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.6 4.08 -8.23 3.09e-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.29e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 9.94e-18
## 6 TypeMississippi:Treatmentchilled NA NA NA NAbroom.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.6 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 NAGLM para efectos fijos (Se implementan dos factores Tiempo y Dieta)
Análisis Descriptivo
data("ChickWeight")
ggplot(ChickWeight,aes(x = Time, y = weight)) + geom_point(aes(color =Diet)) + geom_path(aes(color = Diet, group = Chick))
* Gráficamente se podría determinar que los pollos ganan peso a medida que pasa el tiempo, dependiendo el tipo de dieta que sea aplicada.
* Posiblemente la dieta 1 es la que menor efectividad presenta con respecto a las otras dietas.
* El peso siempre tiene valores están dados en números enteros positivos, posiblemente adecuado para uso de GLM tipo Poisson (Interacción de dieta con el tiempo).
Fit1_Poisson <- glm(weight ~ Diet:Time, data =ChickWeight, family = poisson) # Se agrega familia
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 NA* Se evidencia aumento del peso conforme pasa el tiempo, así como la interacción existente en la que cada dieta es diferente con el factor tiempo. Confirmando lo mencionado en el análisis descriptivo.
sjPlot::plot_model(Fit1_Poisson, type = "eff") #Tipo efecto## Package `effects` is not available, but needed for `ggeffect()`. Either
## install package `effects`, or use `ggpredict()`. Calling `ggpredict()`
## now.## $Diet
##
## $Time
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* La dieta más eficiente con respecto al aumento de peso de los pollos es la dieta 3, teniendo en cuenta que presenta una pendiente mayor en comparación con las dietas 1, 2 y 4 con una estimación de 0.086.
* Se evidencia variación en los valores obtenidos en el intercepto al considerar la identidad de los pollos, notando variabilidad.
* El aumento del peso en los pollos es exponencial con respecto al tiempo, por lo que se puede decir que conforme pasa el tiempo los pollos ganan más peso, teniendo en cuenta, que el aumento no es igual en todas las dietas empleadas.
PARTE 2: MODELO LINEAL MIXTO CON RUIDO
noise <- runif(84, 0.02, 0.04)
noise## [1] 0.03468468 0.02755779 0.02687665 0.02708744 0.03930659 0.03706695
## [7] 0.02258138 0.02042335 0.03723657 0.02761657 0.02038949 0.03188724
## [13] 0.03719543 0.03008547 0.03731303 0.03006879 0.02438412 0.02160773
## [19] 0.02608369 0.02766264 0.02140072 0.03592416 0.02681550 0.03822169
## [25] 0.02668717 0.02825742 0.03352421 0.02834914 0.03902089 0.02349719
## [31] 0.03652049 0.03120365 0.02078254 0.03380533 0.03636345 0.02048352
## [37] 0.02094596 0.02269041 0.03616756 0.02837106 0.03897849 0.02097912
## [43] 0.03184269 0.03048980 0.03904742 0.02194354 0.02214648 0.03857901
## [49] 0.02683370 0.02085537 0.03440495 0.03333030 0.02211637 0.03488640
## [55] 0.03589830 0.02977545 0.02360882 0.03830937 0.02923339 0.02096916
## [61] 0.03821385 0.03903168 0.03128653 0.02367838 0.02255378 0.03523449
## [67] 0.02096292 0.03098422 0.02141483 0.02269463 0.02869803 0.03336503
## [73] 0.02932369 0.03596563 0.02268962 0.02185085 0.02010069 0.03339927
## [79] 0.03607232 0.02006592 0.03285910 0.02457775 0.03883053 0.02822755# Para trabajar con variable Uptake2
uptake2 <- CO2$uptake + noise
(uptake2)## [1] 16.034685 30.427558 34.826877 37.227087 35.339307 39.237067 39.722581
## [8] 13.620423 27.337237 37.127617 41.820389 40.631887 41.437195 44.330085
## [15] 16.237313 32.430069 40.324384 42.121608 42.926084 43.927663 45.521401
## [22] 14.235924 24.126816 30.338222 34.626687 32.528257 35.433524 38.728349
## [29] 9.339021 27.323497 35.036520 38.831204 38.620783 37.533805 42.436363
## [36] 15.120484 21.020946 38.122690 34.036168 38.928371 39.638978 41.420979
## [43] 10.631843 19.230490 26.239047 30.021944 30.922146 32.438579 35.526834
## [50] 12.020855 22.034405 30.633330 31.822116 32.434886 31.135898 31.529775
## [57] 11.323609 19.438309 25.829233 27.920969 28.538214 28.139032 27.831287
## [64] 10.523678 14.922554 18.135234 18.920963 19.530984 22.221415 21.922695
## [71] 7.728698 11.433365 12.329324 13.035966 12.522690 13.721851 14.420101
## [78] 10.633399 18.036072 17.920066 17.932859 17.924578 18.938831 19.928228Análisis Descriptivo
ggplot(CO2,
aes(x = conc, y = uptake2, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment))
###### * Gráficamente, al agregar ruido a los datos las diferencias no
son tan notorias con respecto a la captación de CO2 demostrada
anteriormente (Uptake).Sin embargo, es evidente que el tratamiento
empleado si influye en la captación de carbono, ya que, las plantas no
refrigeradas presentan mayor captación en comparación con las sometidas
a un tratamiento de refrigeración, tanto para las plantas de Quebec como
para las plantas de Mississipi.
* Como es posible observar, se presenta alta variabilidad asociada a cada una de las plantas utilizadas, por lo tanto, para emplear un modelo lineal mixto, es necesario aleatorizar la variable planta (Plant) para que no se generen efectos no deseados sobre la variable respuesta.
Fit3 <- lm(uptake2 ~ conc + Type + Treatment, data = CO2)
Fit4 <- lmer(uptake2 ~ conc + Type + Treatment + (1 | Plant), data = CO2)Modelo Lineal Mixto
broom::glance(Fit3)## # 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.89e-20 3 -270. 551. 563.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>broom::tidy(Fit3)## # 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.67e-14
## 4 Treatmentchilled -6.86 1.35 -5.08 2.45e- 6broom.mixed::glance(Fit4)## # 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(Fit4)## # 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 NAFit3 <- lm(uptake2 ~ I(log(conc)) + Type:Treatment, data = CO2)
Fit4 <- lmer(uptake2 ~ I(log(conc)) + Type:Treatment + (1|Plant), data = CO2)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientbroom::glance(Fit3)## # 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 6.95e-29 4 -246. 504. 519.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>broom::tidy(Fit3) ## # 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.18e-12
## 2 I(log(conc)) 8.48 0.678 12.5 2.02e-20
## 3 TypeQuebec:Treatmentnonchilled 19.5 1.44 13.6 2.59e-22
## 4 TypeMississippi:Treatmentnonchilled 10.1 1.44 7.04 6.22e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 9.83e-18
## 6 TypeMississippi:Treatmentchilled NA NA NA NAbroom.mixed::glance(Fit4)## # 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(Fit4)## # 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.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 NA* La variable individuo (Planta) fue interceptada, tal como se puede apreciar en los resultados obtenidos, con el fin de no alterar los resultados de la variable respuesta.
* Las plantas sometidas a un tratamiento de refrigeración presentan menor captación de carbono, con respecto a las plantas con tratamiento de no refrigeración.
* Las plantas de Mississipi, en general, presentan menor captación de carbono, en comparación con las plantas de Quebec, obteniendo valores estimados de 10.14 vs 19.52, respectivamente.
GLM para efectos fijos
data("ChickWeight")
ChickWeight$Chick %>% length() #Total de pollos## [1] 578ChickWeight$Chick %>% unique ()%>% length #Tipos de pollos## [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 1* Al obtener valores enteros positivos, se puede aplicar GML tipo Poisson.
ggplot(ChickWeight,aes(x = Time, y = weight)) + geom_point(aes(color =Diet)) + geom_path(aes(color = Diet, group = Chick))
###### * El comportamiento de los pollos es diferente a pesar de estar
sometidos a dietas semejantes. Por lo que es posible aplicar un GML.
Fit3_Poisson <- glm(weight ~ Diet:Time, data =ChickWeight, family = poisson)
Fit4_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(Fit3_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(Fit4_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