Modelo lineales mixtos (variables aleatorias)
Oscar A. Gómez
25/6/2023
El objtivo del ejecicio es algunas variaciones que puede tener las interpretaciones de un modelo lineal a un modelo donde se tiene en cuenta una variable aleatoria
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.1.3## Warning: package 'ggplot2' was built under R version 4.1.3## Warning: package 'tibble' was built under R version 4.1.3## Warning: package 'tidyr' was built under R version 4.1.3## Warning: package 'readr' was built under R version 4.1.3## Warning: package 'purrr' was built under R version 4.1.3## Warning: package 'dplyr' was built under R version 4.1.3## Warning: package 'stringr' was built under R version 4.1.3## Warning: package 'forcats' was built under R version 4.1.3## Warning: package 'lubridate' was built under R version 4.1.3## -- 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.1.3library(lme4)## Warning: package 'lme4' was built under R version 4.1.3## Loading required package: Matrix## Warning: package 'Matrix' was built under R version 4.1.3##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(ggplot2)Para esto llamamos algunas librerias que nos ayudara usar una base de datos tidy, la cual nos permitira una discriminación de los datos por observaciones por cada individuo (siendo nuestro principal interés la planta) y así poder una comparación de un modelo lineal a un modelo lineal mixto.
Ahora bien, los datos con los que vamos trabajar estan enfocados a la captación de CO2 de una especie de pasto de dos sitios (Quebec y Mississippi) a diferentes niveles de concentración de CO2 (95,175, 250,350,500, 675 y 1000). En este ejercicio la mitad de las plantas recibieron un enfriamiento durante la noche previo a la realización del experimento y la otra mitad no se realizo dicho procedimiento.
data("CO2")
datos=CO2
datos2=datos %>% mutate(uptake2= uptake + runif(84,0.02,0.04))
datos2 ## Grouped Data: uptake ~ conc | Plant
## Plant Type Treatment conc uptake uptake2
## 1 Qn1 Quebec nonchilled 95 16.0 16.033301
## 2 Qn1 Quebec nonchilled 175 30.4 30.433364
## 3 Qn1 Quebec nonchilled 250 34.8 34.826407
## 4 Qn1 Quebec nonchilled 350 37.2 37.231181
## 5 Qn1 Quebec nonchilled 500 35.3 35.331208
## 6 Qn1 Quebec nonchilled 675 39.2 39.233870
## 7 Qn1 Quebec nonchilled 1000 39.7 39.724019
## 8 Qn2 Quebec nonchilled 95 13.6 13.632533
## 9 Qn2 Quebec nonchilled 175 27.3 27.337327
## 10 Qn2 Quebec nonchilled 250 37.1 37.122287
## 11 Qn2 Quebec nonchilled 350 41.8 41.830542
## 12 Qn2 Quebec nonchilled 500 40.6 40.621048
## 13 Qn2 Quebec nonchilled 675 41.4 41.428868
## 14 Qn2 Quebec nonchilled 1000 44.3 44.321618
## 15 Qn3 Quebec nonchilled 95 16.2 16.233135
## 16 Qn3 Quebec nonchilled 175 32.4 32.424431
## 17 Qn3 Quebec nonchilled 250 40.3 40.322911
## 18 Qn3 Quebec nonchilled 350 42.1 42.123033
## 19 Qn3 Quebec nonchilled 500 42.9 42.933405
## 20 Qn3 Quebec nonchilled 675 43.9 43.927918
## 21 Qn3 Quebec nonchilled 1000 45.5 45.539384
## 22 Qc1 Quebec chilled 95 14.2 14.222619
## 23 Qc1 Quebec chilled 175 24.1 24.139661
## 24 Qc1 Quebec chilled 250 30.3 30.339143
## 25 Qc1 Quebec chilled 350 34.6 34.635674
## 26 Qc1 Quebec chilled 500 32.5 32.522596
## 27 Qc1 Quebec chilled 675 35.4 35.436848
## 28 Qc1 Quebec chilled 1000 38.7 38.729574
## 29 Qc2 Quebec chilled 95 9.3 9.336692
## 30 Qc2 Quebec chilled 175 27.3 27.329999
## 31 Qc2 Quebec chilled 250 35.0 35.038788
## 32 Qc2 Quebec chilled 350 38.8 38.823389
## 33 Qc2 Quebec chilled 500 38.6 38.634258
## 34 Qc2 Quebec chilled 675 37.5 37.525838
## 35 Qc2 Quebec chilled 1000 42.4 42.439918
## 36 Qc3 Quebec chilled 95 15.1 15.127012
## 37 Qc3 Quebec chilled 175 21.0 21.035854
## 38 Qc3 Quebec chilled 250 38.1 38.120143
## 39 Qc3 Quebec chilled 350 34.0 34.026317
## 40 Qc3 Quebec chilled 500 38.9 38.932325
## 41 Qc3 Quebec chilled 675 39.6 39.638946
## 42 Qc3 Quebec chilled 1000 41.4 41.423030
## 43 Mn1 Mississippi nonchilled 95 10.6 10.632470
## 44 Mn1 Mississippi nonchilled 175 19.2 19.224016
## 45 Mn1 Mississippi nonchilled 250 26.2 26.228520
## 46 Mn1 Mississippi nonchilled 350 30.0 30.034147
## 47 Mn1 Mississippi nonchilled 500 30.9 30.932403
## 48 Mn1 Mississippi nonchilled 675 32.4 32.422861
## 49 Mn1 Mississippi nonchilled 1000 35.5 35.536179
## 50 Mn2 Mississippi nonchilled 95 12.0 12.024597
## 51 Mn2 Mississippi nonchilled 175 22.0 22.028053
## 52 Mn2 Mississippi nonchilled 250 30.6 30.636523
## 53 Mn2 Mississippi nonchilled 350 31.8 31.837113
## 54 Mn2 Mississippi nonchilled 500 32.4 32.423034
## 55 Mn2 Mississippi nonchilled 675 31.1 31.121907
## 56 Mn2 Mississippi nonchilled 1000 31.5 31.536667
## 57 Mn3 Mississippi nonchilled 95 11.3 11.335626
## 58 Mn3 Mississippi nonchilled 175 19.4 19.427659
## 59 Mn3 Mississippi nonchilled 250 25.8 25.839482
## 60 Mn3 Mississippi nonchilled 350 27.9 27.929525
## 61 Mn3 Mississippi nonchilled 500 28.5 28.539691
## 62 Mn3 Mississippi nonchilled 675 28.1 28.127650
## 63 Mn3 Mississippi nonchilled 1000 27.8 27.836101
## 64 Mc1 Mississippi chilled 95 10.5 10.528489
## 65 Mc1 Mississippi chilled 175 14.9 14.932507
## 66 Mc1 Mississippi chilled 250 18.1 18.139590
## 67 Mc1 Mississippi chilled 350 18.9 18.935747
## 68 Mc1 Mississippi chilled 500 19.5 19.524247
## 69 Mc1 Mississippi chilled 675 22.2 22.232228
## 70 Mc1 Mississippi chilled 1000 21.9 21.935687
## 71 Mc2 Mississippi chilled 95 7.7 7.728417
## 72 Mc2 Mississippi chilled 175 11.4 11.429615
## 73 Mc2 Mississippi chilled 250 12.3 12.328193
## 74 Mc2 Mississippi chilled 350 13.0 13.028806
## 75 Mc2 Mississippi chilled 500 12.5 12.531404
## 76 Mc2 Mississippi chilled 675 13.7 13.737137
## 77 Mc2 Mississippi chilled 1000 14.4 14.429513
## 78 Mc3 Mississippi chilled 95 10.6 10.620625
## 79 Mc3 Mississippi chilled 175 18.0 18.024320
## 80 Mc3 Mississippi chilled 250 17.9 17.930472
## 81 Mc3 Mississippi chilled 350 17.9 17.939499
## 82 Mc3 Mississippi chilled 500 17.9 17.925626
## 83 Mc3 Mississippi chilled 675 18.9 18.933397
## 84 Mc3 Mississippi chilled 1000 19.9 19.930346ggplot(datos2, aes(x=conc, y=uptake2, color=Type))+
geom_point(aes(shape=Treatment))+
geom_path(aes(group = Plant, lty=Treatment))+
theme_bw()
La grafica reflaja el compotamiento de la captación de CO2 de las plantas de pastos discriminada por los sitios Quebec y Missippi, además por plantas a las que se le realizo enfriamiento o no. Nos referiremos como tratamiento si se realizo enfriamiento o no y tipo a los lugares de procedencia de las plantas.
La grafica nos da cuenta que las plantas captan más CO2 a medida que se aumenta los niveles de concentración del mismo, primero de forma exponencial y despues se estabiliza aunmentado paulatinamente. Captan más CO2 las plantas de Quebec con respecto a Missippi; del mismo modo, podemos afirmar que las plantas a las que no se aplico el enfriamiento se captó más CO2. Entonces así se puede afirmar que las plantas procedentes de quebec y que no se les aplico enfriamiento obtuvieron una captación mayor de CO2, mientras que las plantas que proceden de Mississipi y que recibieron el enfriamiento captaron menos CO2.
Para la comparación de los modelos vamos a aplicarle a los datos las funciones glace y tidy del paquete bloom que nos ayudara primero ordenar los datos del tal forma pudieramos obtener metricas de calidad de los modelos, además de tener un resumen estadistico con los parametros de estimación, desviación estandar, estadistico y el p-value.
Se va utilizar bloom y lm para el caso de modelo lineal y bloom.mixed y lmer para el modelo mixto donde se contempla la variable aleatoria de planta.
fit1= lm(uptake2~I(log(conc))+Type:Treatment, data = datos2)
fit2= lmer(uptake2~I(log(conc))+Type:Treatment + (1| Plant), data = datos2)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientbroom::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 6.91e-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.23 3.16e-12
## 2 I(log(conc)) 8.48 0.678 12.5 2.00e-20
## 3 TypeQuebec:Treatmentnonchilled 19.5 1.44 13.6 2.59e-22
## 4 TypeMississippi:Treatmentnonchilled 10.1 1.44 7.04 6.24e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 9.82e-18
## 6 TypeMississippi:Treatmentchilled NA NA NA NAbroom.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) ## # 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 NAPodemos ver que hay diferencias entre ambos modelos si utlizamos uno o el otro. Es de resaltar que analizando el modelo mixto la variable aleatoria nos proporciona información adicional en el experimento, teniendo una estimación propia cada planta considerandose así una división entre efectos fijos y efectos aleatorios.
Como conclusiones de la comparación de los modelos en primer lugar, se observa una variación de la interpretación de los datos si se considera un modelos lineal simbre o un modelo lineal mixto. En segundo lugar, al considerarse un modelo mixto, se emplea una variable aleatoria que nos dara información adicional al otro modelo, siendo por un lado la discriminación por observaciones o individuos (en este caso la planta) llevandonos a una interpretación que tiene en cuenta las variación de la captación de CO2 por planta y no el promedio del grupo en general. Además, puede aplicarse el experimento en otros conceptos al contemplar una variación entre diferentes sujetos.Por ultimo, en experimentos relacionados con la agricultura resulta ventajoso utilizar un modelo mixto donde contemplen la variación por planta de los tratamientos ya que entendemos que las plantas se pueden comportar de manera diferente y entre más entendamos y registremos esta diferencia podremos entender la gama de variaciones que se puedan presentar en campo que puede como respuesta la aplicación de los tratamientos en determinados sujetos. De esta manera, se puede llegar a interpretaciones más contundentes referentes a diferentes problematicas de interes agronomico.
datos3=ChickWeight
datos3%>%length()## [1] 4datos3%>%unique()%>%length()## [1] 4ggplot(datos3, aes(x=Time, y=weight))+
geom_point(aes(color= Diet))+
geom_path(aes(color= Diet, group= Chick))
fit1_Poisson = glm(weight ~ Diet:Time, data = datos3, family = poisson)
fit2_Poisson = glmer(weight ~ Diet:Time + (1|Chick), data = datos3, 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)## # 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::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 NA