Trabajo Final. Modelos Mixtos
Maria Jose Duarte H y Sergio Alejandro Castellanos
2023-06-22
MODELOS MIXTOS
Se utiliza para cuantificar la variación que genera una variable aleatoria en la variable respuesta.
Se tiene en cuenta asociación entre observaciones que se toman al mismo sujeto.
Se tienen crecimientos lineales en los que cada sujeto tiene su propio intercepto.
Se utiliza la técnica de los minimos cuadrados para hacer la regresión lineal y estimar los parametros (intercepto y pendiente)
Desarrollo
Se hace uso de una base de datos predeterminada a la cual se le agrega un “ruido” utilizando la función runif. Esta base de datos se caracteriza por tener plantas en dos lugares, sometidas a distintas concentraciones de CO2 y a tratamientos como el enfriamiento o no enfriamiento y se analiza la capatación de CO2.
#LIibrerias necesarias
library(Matrix)
library(lme4) #Esta se utiliza para modelos aleatorios.
#Base de datos a usar:
datos = data("CO2")
print(CO2)## Grouped Data: uptake ~ conc | Plant
## Plant Type Treatment conc uptake
## 1 Qn1 Quebec nonchilled 95 16.0
## 2 Qn1 Quebec nonchilled 175 30.4
## 3 Qn1 Quebec nonchilled 250 34.8
## 4 Qn1 Quebec nonchilled 350 37.2
## 5 Qn1 Quebec nonchilled 500 35.3
## 6 Qn1 Quebec nonchilled 675 39.2
## 7 Qn1 Quebec nonchilled 1000 39.7
## 8 Qn2 Quebec nonchilled 95 13.6
## 9 Qn2 Quebec nonchilled 175 27.3
## 10 Qn2 Quebec nonchilled 250 37.1
## 11 Qn2 Quebec nonchilled 350 41.8
## 12 Qn2 Quebec nonchilled 500 40.6
## 13 Qn2 Quebec nonchilled 675 41.4
## 14 Qn2 Quebec nonchilled 1000 44.3
## 15 Qn3 Quebec nonchilled 95 16.2
## 16 Qn3 Quebec nonchilled 175 32.4
## 17 Qn3 Quebec nonchilled 250 40.3
## 18 Qn3 Quebec nonchilled 350 42.1
## 19 Qn3 Quebec nonchilled 500 42.9
## 20 Qn3 Quebec nonchilled 675 43.9
## 21 Qn3 Quebec nonchilled 1000 45.5
## 22 Qc1 Quebec chilled 95 14.2
## 23 Qc1 Quebec chilled 175 24.1
## 24 Qc1 Quebec chilled 250 30.3
## 25 Qc1 Quebec chilled 350 34.6
## 26 Qc1 Quebec chilled 500 32.5
## 27 Qc1 Quebec chilled 675 35.4
## 28 Qc1 Quebec chilled 1000 38.7
## 29 Qc2 Quebec chilled 95 9.3
## 30 Qc2 Quebec chilled 175 27.3
## 31 Qc2 Quebec chilled 250 35.0
## 32 Qc2 Quebec chilled 350 38.8
## 33 Qc2 Quebec chilled 500 38.6
## 34 Qc2 Quebec chilled 675 37.5
## 35 Qc2 Quebec chilled 1000 42.4
## 36 Qc3 Quebec chilled 95 15.1
## 37 Qc3 Quebec chilled 175 21.0
## 38 Qc3 Quebec chilled 250 38.1
## 39 Qc3 Quebec chilled 350 34.0
## 40 Qc3 Quebec chilled 500 38.9
## 41 Qc3 Quebec chilled 675 39.6
## 42 Qc3 Quebec chilled 1000 41.4
## 43 Mn1 Mississippi nonchilled 95 10.6
## 44 Mn1 Mississippi nonchilled 175 19.2
## 45 Mn1 Mississippi nonchilled 250 26.2
## 46 Mn1 Mississippi nonchilled 350 30.0
## 47 Mn1 Mississippi nonchilled 500 30.9
## 48 Mn1 Mississippi nonchilled 675 32.4
## 49 Mn1 Mississippi nonchilled 1000 35.5
## 50 Mn2 Mississippi nonchilled 95 12.0
## 51 Mn2 Mississippi nonchilled 175 22.0
## 52 Mn2 Mississippi nonchilled 250 30.6
## 53 Mn2 Mississippi nonchilled 350 31.8
## 54 Mn2 Mississippi nonchilled 500 32.4
## 55 Mn2 Mississippi nonchilled 675 31.1
## 56 Mn2 Mississippi nonchilled 1000 31.5
## 57 Mn3 Mississippi nonchilled 95 11.3
## 58 Mn3 Mississippi nonchilled 175 19.4
## 59 Mn3 Mississippi nonchilled 250 25.8
## 60 Mn3 Mississippi nonchilled 350 27.9
## 61 Mn3 Mississippi nonchilled 500 28.5
## 62 Mn3 Mississippi nonchilled 675 28.1
## 63 Mn3 Mississippi nonchilled 1000 27.8
## 64 Mc1 Mississippi chilled 95 10.5
## 65 Mc1 Mississippi chilled 175 14.9
## 66 Mc1 Mississippi chilled 250 18.1
## 67 Mc1 Mississippi chilled 350 18.9
## 68 Mc1 Mississippi chilled 500 19.5
## 69 Mc1 Mississippi chilled 675 22.2
## 70 Mc1 Mississippi chilled 1000 21.9
## 71 Mc2 Mississippi chilled 95 7.7
## 72 Mc2 Mississippi chilled 175 11.4
## 73 Mc2 Mississippi chilled 250 12.3
## 74 Mc2 Mississippi chilled 350 13.0
## 75 Mc2 Mississippi chilled 500 12.5
## 76 Mc2 Mississippi chilled 675 13.7
## 77 Mc2 Mississippi chilled 1000 14.4
## 78 Mc3 Mississippi chilled 95 10.6
## 79 Mc3 Mississippi chilled 175 18.0
## 80 Mc3 Mississippi chilled 250 17.9
## 81 Mc3 Mississippi chilled 350 17.9
## 82 Mc3 Mississippi chilled 500 17.9
## 83 Mc3 Mississippi chilled 675 18.9
## 84 Mc3 Mississippi chilled 1000 19.9#Esta muestra la absorción de CO2 en plantas en diferentes lugares y con tratamiento de enfriamiento. Vamos a trabajar con la variable uptake
A esta variable se le adiciona un poco de ruido con la ayuda de (runif) para generar variabilidad
#upatake2 = uptake + ruido
ruido = runif(n=84,0.02,0.04)
uptake2 = CO2$uptake + ruido
print(uptake2)## [1] 16.030128 30.433578 34.822584 37.229631 35.323230 39.227688 39.736748
## [8] 13.621314 27.333905 37.134534 41.825210 40.631532 41.436327 44.322050
## [15] 16.229495 32.420357 40.333937 42.127929 42.937961 43.921575 45.528202
## [22] 14.225966 24.126640 30.325413 34.631917 32.534359 35.428772 38.727273
## [29] 9.337069 27.321252 35.032647 38.828974 38.622903 37.527127 42.435219
## [36] 15.139540 21.021917 38.121975 34.027715 38.935183 39.626032 41.431095
## [43] 10.627710 19.221778 26.222157 30.031944 30.938281 32.431195 35.530366
## [50] 12.024993 22.036789 30.633208 31.836024 32.436642 31.135216 31.533650
## [57] 11.332929 19.438922 25.834227 27.934525 28.536579 28.120536 27.822741
## [64] 10.536651 14.923521 18.121386 18.920102 19.534854 22.231464 21.927602
## [71] 7.738414 11.436318 12.328901 13.023921 12.526749 13.724507 14.424849
## [78] 10.630388 18.037466 17.928960 17.920037 17.935075 18.923248 19.936569#son 84 datos en totalVisualización: se hace un grafico con el fin de identificar como se distribuyen los resultados por Type (lugar) y Trearment (enfriada o no)
library(ggplot2)
ggplot(CO2, aes(x = conc, y = uptake2, color = Type)) +
geom_point(aes(shape = Treatment)) +theme_bw()
En este caso se observa que las plantas de Quebec tienen mayor absorción de CO2 que las plantas de Mississippi
Tambien se tiene que en mayores concentraciones hay mayor variabilidad
library(ggplot2)
ggplot(CO2, aes(x = conc, y = uptake2 , color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment)) +
theme_bw()
#agrupadas por planta, tipo de linea según el tratamiento)Observamos planta como variable aleatoria, el objetivo del estudio no es detectar las diferencias de los niveles del factor PLANTA.
Hay algunas tendencias en ciertos resultados que no eran esperados, esto quiere decir que hay factores que la afectan pero que no se encuentran dentro de las variables de efecto fijo (lugar, tratamiento)
A partir de esta información se tiene que cada planta se comporta como un individuo diferente y puede tener un efecto en la captación (variable respuesta)
las plantas que presentan mayor captacion de dioxido de carbono son aquellas a las que no se les trato con enfriamiento tanto para Quebec como para Mississippi. Sin embargo, globalmente las que tuvieron mejor desempeño fueron las plantas que se encontraban en Quebec.
En el ejemplo se observa que la concentración toma forma logaritmica por lo que se plantean los modelos teniendo en cuenta esto:
Modelos
library(lme4)
#MODELO LINEAL SIMPLE
Fit1 <- lm(uptake2 ~ I(log(conc)) + Type:Treatment, data = CO2)
summary(Fit1)##
## Call:
## lm(formula = uptake2 ~ I(log(conc)) + Type:Treatment, data = CO2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.7090 -2.9015 0.5859 2.7530 8.8674
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -33.5247 4.0773 -8.222 3.22e-12 ***
## I(log(conc)) 8.4835 0.6784 12.504 < 2e-16 ***
## TypeQuebec:Treatmentnonchilled 19.5189 1.4405 13.550 < 2e-16 ***
## TypeMississippi:Treatmentnonchilled 10.1404 1.4405 7.040 6.27e-10 ***
## TypeQuebec:Treatmentchilled 15.9380 1.4405 11.064 < 2e-16 ***
## TypeMississippi:Treatmentchilled NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.668 on 79 degrees of freedom
## Multiple R-squared: 0.8227, Adjusted R-squared: 0.8137
## F-statistic: 91.63 on 4 and 79 DF, p-value: < 2.2e-16#se utiliza unicamente para variables de efecto fijo.
#siendo uptake 2 la variable respuesta que depende de la concentración de CO2, y la interacción entre el tratamiento y el lugar.
#Planteamos interacción en el ejemplo al observar la gráfica.
#MODELO CON EL FACTOR ALEATORIO (PLANTA)
Fit2 <- lmer(uptake2 ~ I(log(conc)) + Type : Treatment + (1 | Plant), data = CO2)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientsummary(Fit2)## Linear mixed model fit by REML ['lmerMod']
## Formula: uptake2 ~ I(log(conc)) + Type:Treatment + (1 | Plant)
## Data: CO2
##
## REML criterion at convergence: 482.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6914 -0.5385 0.1042 0.6894 1.8911
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant (Intercept) 2.157 1.469
## Residual 20.258 4.501
## Number of obs: 84, groups: Plant, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -33.5247 4.0221 -8.335
## I(log(conc)) 8.4835 0.6542 12.968
## TypeQuebec:Treatmentnonchilled 19.5189 1.8350 10.637
## TypeMississippi:Treatmentnonchilled 10.1404 1.8350 5.526
## TypeQuebec:Treatmentchilled 15.9380 1.8350 8.686
##
## Correlation of Fixed Effects:
## (Intr) I(l()) TypQbc:Trtmntn TypM:T
## I(log(cnc)) -0.947
## TypQbc:Trtmntn -0.228 0.000
## TypMssssp:T -0.228 0.000 0.500
## TypQbc:Trtmntc -0.228 0.000 0.500 0.500
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient#se define como factor aleatorio (1 | Plant)
#Se tiene en cuenta errores aleatorios por individuo.
#Cada planta con un intercepto distinto. Se observa que los interceptos y pendientes de cada modelos son muy similares
En este caso para las variables aleatorias la varianza es de 2.15
# MODELO SIMPLE
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.6 7.08e-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.22e-12
## 2 I(log(conc)) 8.48 0.678 12.5 2.05e-20
## 3 TypeQuebec:Treatmentnonchilled 19.5 1.44 13.6 2.63e-22
## 4 TypeMississippi:Treatmentnonchilled 10.1 1.44 7.04 6.27e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 1.00e-17
## 6 TypeMississippi:Treatmentchilled NA NA NA NA# MODELOS MIXTOS
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.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 NALos resultados entre el modelo simple y el modelo aleatorio no son muy distintos y se analizan de la misma forma
En los valores mixtos NO existe el parametro de r cuadrado.
Se observan efectos fijos y no aleatorios
los resultados indican lo que se habia planteado en la interpretacion grafica, donde el tratamiento con enfriamiento en Mississippi es el que tiene menor captacion de dioxido de carbono.
Además, tambien demostró que los mejores reultados son en Quebec, teniendo en cuenta que el tratamiento con enfriamietno disminuye la captación de dioxido de carbono. Sin embargo, sigue siendo mayor a los resultados de Mississippi con el tratamiento de no enfriar
Efectos glmer
library(ggplot2)
#Vamos a recurrir a otros datos.
#El peso es explicado por la dieta en interacción con la dieta
data("ChickWeight")
Fit1_Poisson <- glm(weight ~ Diet:Time, data =ChickWeight, family = poisson)
#Planteamos interacción entre dieta y tiempo.
#Podemos utilizar poisson porque los valores de peso siempre son positivos
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?#Errores aleatorios de la variable pollo.
ggplot(ChickWeight, aes(x = Time, y = weight)) +
geom_point(aes(color = Diet)) +
geom_path(aes(color = Diet, group = Chick))
#El peso aumenta con el tiempo.
#Tenemos distintas pendientes para cada dieta.
#La linea lo de chick (pollo), es la variable aleatoria.
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#plot de tipo efecto de Se observan diferencias en el intercepto entre cada modelo.
La identidad del individuo genera alguna variabilidad.
A pesar de similaridades en la dieta y tiempo los pollos tienen crecimientos diferentes.
Qué importancia tiene este modelo en la agricultura
Es util debido a que este contribuye a un mejor estudio de las plantas dentro de un cultivo, al interpretarlas como individuos que tienen comportamientos diferentes y pueden influenciar la variable respuesta de la investigación.
es util en ensayos que incluyan factores fijos o controlados por el investigador y factores aleatorios como el individuo para cuantificar la variación de la relación variables fijas -variable respuesta entre la variable aleatoria.