MODELO LINEAL MIXTO
Para llevar a cabo el análisis se descargan las siguientes
librerias:
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 errors
library(broom.mixed)
## Warning: package 'broom.mixed' was built under R version 4.3.1
library(lme4)
## 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, unpack
1. CREACIÓN DE BASE DE DATOS
Se usa una base de datos preestablecida donde se tienen datos de la
captación de CO2 de las plantas sometidas a 7 diferentes tipos de
concentraciones de CO2 (95,150,250,350,500,675,1000) y a distintas
temperaturas, que sería el tratamiento experimental (chilled y
nonchilled). No obstante, se le introduce ruido a los datos con la
función “runif”.
#Base de datos original.
data("CO2")
view(CO2)
#Ruido de los datos.
ruido <- runif(84,0.02,0.04)
#Variable con el ruido con la que se va a trabajar.
uptake2= CO2$uptake + ruido
datos = data.frame(CO2$Plant, CO2$Treatment, CO2$Type, CO2$conc, uptake2)
datos
## CO2.Plant CO2.Treatment CO2.Type CO2.conc uptake2
## 1 Qn1 nonchilled Quebec 95 16.038466
## 2 Qn1 nonchilled Quebec 175 30.433421
## 3 Qn1 nonchilled Quebec 250 34.834564
## 4 Qn1 nonchilled Quebec 350 37.226392
## 5 Qn1 nonchilled Quebec 500 35.328485
## 6 Qn1 nonchilled Quebec 675 39.235105
## 7 Qn1 nonchilled Quebec 1000 39.730434
## 8 Qn2 nonchilled Quebec 95 13.633556
## 9 Qn2 nonchilled Quebec 175 27.339666
## 10 Qn2 nonchilled Quebec 250 37.121054
## 11 Qn2 nonchilled Quebec 350 41.832767
## 12 Qn2 nonchilled Quebec 500 40.633768
## 13 Qn2 nonchilled Quebec 675 41.422844
## 14 Qn2 nonchilled Quebec 1000 44.321721
## 15 Qn3 nonchilled Quebec 95 16.224039
## 16 Qn3 nonchilled Quebec 175 32.420252
## 17 Qn3 nonchilled Quebec 250 40.330528
## 18 Qn3 nonchilled Quebec 350 42.126811
## 19 Qn3 nonchilled Quebec 500 42.920128
## 20 Qn3 nonchilled Quebec 675 43.936873
## 21 Qn3 nonchilled Quebec 1000 45.525696
## 22 Qc1 chilled Quebec 95 14.235426
## 23 Qc1 chilled Quebec 175 24.134800
## 24 Qc1 chilled Quebec 250 30.337718
## 25 Qc1 chilled Quebec 350 34.620581
## 26 Qc1 chilled Quebec 500 32.531210
## 27 Qc1 chilled Quebec 675 35.422487
## 28 Qc1 chilled Quebec 1000 38.720335
## 29 Qc2 chilled Quebec 95 9.338500
## 30 Qc2 chilled Quebec 175 27.333540
## 31 Qc2 chilled Quebec 250 35.020420
## 32 Qc2 chilled Quebec 350 38.837153
## 33 Qc2 chilled Quebec 500 38.626433
## 34 Qc2 chilled Quebec 675 37.536330
## 35 Qc2 chilled Quebec 1000 42.437641
## 36 Qc3 chilled Quebec 95 15.136937
## 37 Qc3 chilled Quebec 175 21.028653
## 38 Qc3 chilled Quebec 250 38.139335
## 39 Qc3 chilled Quebec 350 34.034995
## 40 Qc3 chilled Quebec 500 38.931279
## 41 Qc3 chilled Quebec 675 39.629196
## 42 Qc3 chilled Quebec 1000 41.437804
## 43 Mn1 nonchilled Mississippi 95 10.621642
## 44 Mn1 nonchilled Mississippi 175 19.228127
## 45 Mn1 nonchilled Mississippi 250 26.234757
## 46 Mn1 nonchilled Mississippi 350 30.039172
## 47 Mn1 nonchilled Mississippi 500 30.930955
## 48 Mn1 nonchilled Mississippi 675 32.432354
## 49 Mn1 nonchilled Mississippi 1000 35.533045
## 50 Mn2 nonchilled Mississippi 95 12.026602
## 51 Mn2 nonchilled Mississippi 175 22.034150
## 52 Mn2 nonchilled Mississippi 250 30.626697
## 53 Mn2 nonchilled Mississippi 350 31.829506
## 54 Mn2 nonchilled Mississippi 500 32.428483
## 55 Mn2 nonchilled Mississippi 675 31.136885
## 56 Mn2 nonchilled Mississippi 1000 31.527150
## 57 Mn3 nonchilled Mississippi 95 11.327043
## 58 Mn3 nonchilled Mississippi 175 19.424919
## 59 Mn3 nonchilled Mississippi 250 25.838789
## 60 Mn3 nonchilled Mississippi 350 27.929741
## 61 Mn3 nonchilled Mississippi 500 28.530645
## 62 Mn3 nonchilled Mississippi 675 28.121267
## 63 Mn3 nonchilled Mississippi 1000 27.826076
## 64 Mc1 chilled Mississippi 95 10.536184
## 65 Mc1 chilled Mississippi 175 14.929098
## 66 Mc1 chilled Mississippi 250 18.128453
## 67 Mc1 chilled Mississippi 350 18.926850
## 68 Mc1 chilled Mississippi 500 19.535367
## 69 Mc1 chilled Mississippi 675 22.234847
## 70 Mc1 chilled Mississippi 1000 21.921293
## 71 Mc2 chilled Mississippi 95 7.722971
## 72 Mc2 chilled Mississippi 175 11.427182
## 73 Mc2 chilled Mississippi 250 12.325896
## 74 Mc2 chilled Mississippi 350 13.026802
## 75 Mc2 chilled Mississippi 500 12.520670
## 76 Mc2 chilled Mississippi 675 13.721461
## 77 Mc2 chilled Mississippi 1000 14.429965
## 78 Mc3 chilled Mississippi 95 10.628221
## 79 Mc3 chilled Mississippi 175 18.031431
## 80 Mc3 chilled Mississippi 250 17.922757
## 81 Mc3 chilled Mississippi 350 17.921374
## 82 Mc3 chilled Mississippi 500 17.927478
## 83 Mc3 chilled Mississippi 675 18.935650
## 84 Mc3 chilled Mississippi 1000 19.935284
#Variable aleatoria: Plant.
view(uptake2)
2. ANÁLISIS DESCRIPTIVO
#Grafica
ggplot(CO2, aes(x = conc, y = uptake2, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment)) +
theme_bw()
La grafica muestra la captación de CO2 de cada planta al ser
sometidas a diferentes concentraciones y a diferentes tratamientos
(nonchilled y chilled), asimismo, se podría observar que en general las
plantas chilled captaron menos CO2 en comparación con las plantas
nonchilled. No obstante, se observa una diferencia entre lugares, donde,
Quebec presenta plantas nonchilled con mayor captación de CO2.
Por otro lado, se observa una diferencia entre las plantas de Quebec
y las plantas de Mississippi, en donde las de Quebec presentan una mayor
captación de CO2 que las de Mississippi.
3. MODELO LINEAL VS MODELO LINEAL MIXTO
#NOTA: Como se tiene un modelo exponencial se usa la función (Log).
#Modelo lineal
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.711 -2.890 0.587 2.762 8.882
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -33.5238 4.0768 -8.223 3.20e-12 ***
## I(log(conc)) 8.4832 0.6784 12.505 < 2e-16 ***
## TypeQuebec:Treatmentnonchilled 19.5203 1.4403 13.553 < 2e-16 ***
## TypeMississippi:Treatmentnonchilled 10.1399 1.4403 7.040 6.26e-10 ***
## TypeQuebec:Treatmentchilled 15.9420 1.4403 11.069 < 2e-16 ***
## TypeMississippi:Treatmentchilled NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.667 on 79 degrees of freedom
## Multiple R-squared: 0.8227, Adjusted R-squared: 0.8138
## F-statistic: 91.67 on 4 and 79 DF, p-value: < 2.2e-16
Para analizar los modelos aleatorios se usa:
#Modelo lineal mixto
#Aleatorizar los datos, factor aleatorio (1| Plant).Esto hace que cada una de las plantas tenga un intercepto distinto.
Fit2 <- lmer(uptake2 ~ I(log(conc)) + Type : Treatment + (1 | Plant), data = CO2)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
summary(Fit2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: uptake2 ~ I(log(conc)) + Type:Treatment + (1 | Plant)
## Data: CO2
##
## REML criterion at convergence: 482.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6926 -0.5375 0.1045 0.6900 1.8943
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant (Intercept) 2.162 1.47
## Residual 20.249 4.50
## Number of obs: 84, groups: Plant, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -33.5238 4.0214 -8.336
## I(log(conc)) 8.4832 0.6541 12.970
## TypeQuebec:Treatmentnonchilled 19.5203 1.8357 10.634
## TypeMississippi:Treatmentnonchilled 10.1399 1.8357 5.524
## TypeQuebec:Treatmentchilled 15.9420 1.8357 8.685
##
## Correlation of Fixed Effects:
## (Intr) I(l()) TypQbc:Trtmntn TypM:T
## I(log(cnc)) -0.946
## 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
La diferencia entre los interceptos indica que las plantas responden
de manera distinta ante los tratamientos.
MODELO LINEAL.
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.00e-29 4 -246. 504. 519.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>
Debido a que se tiene un modelo lineal se analiza el r-squared.Para
decidir si es un buen modelo ajustado, el valor nos debe cercano a 1. En
este caso, nuestro valor da 0.82, lo que se considera un buen ajuste
lineal.
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.20e-12
## 2 I(log(conc)) 8.48 0.678 12.5 2.04e-20
## 3 TypeQuebec:Treatmentnonchilled 19.5 1.44 13.6 2.60e-22
## 4 TypeMississippi:Treatmentnonchilled 10.1 1.44 7.04 6.26e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 9.83e-18
## 6 TypeMississippi:Treatmentchilled NA NA NA NA
#El tidy nos indica cual sera el estimador de cada uno de los factores.
Modelo lineal mixto
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. 77
#En este modelo no se tienen R cuadrados.
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.84 10.6
## 4 fixed <NA> TypeMississippi:Treatmentnonch… 10.1 1.84 5.52
## 5 fixed <NA> TypeQuebec:Treatmentchilled 15.9 1.84 8.68
## 6 ran_pars Plant sd__(Intercept) 1.47 NA NA
## 7 ran_pars Residual sd__Observation 4.50 NA NA
Se observa los valores de las interacciones para la columna estimate
y se puede decir que se cumple con el análisis descriptivo, debido a que
Quebec-nonchilled presenta mayor valor (19.52) y en la grafica concuerda
con que son las plantas que más captan CO2. Por otro lado, Mississipi -
nonchilled presentan menor valor (10.14) y en la grafica concuerda que
captan menos CO2 en comparación con las de Quebec.
EFECTOS FIJOS GLMER.
data("ChickWeight")
#Total pollos
ChickWeight$Chick %>% length()
## [1] 578
#Tipos de pollos
ChickWeight$Chick %>% unique ()%>% length
## [1] 50
head(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
NOTA: Como se tienen valores positivos en el peso y no se tiene
decimales, por tal motivo, se podría analizar como un gml tipo
poisson.
ggplot(ChickWeight, aes(x = Time, y = weight)) +
geom_point(aes(color = Diet)) +
geom_path(aes(color = Diet, group = Chick))
Hay pollos donde a pesar de que se tienen dietas simialres tienen
comportamientos distintos. Esto indica que si se tiene estas medidades
repetitas se puede hacer un glmer.
#gml tipo poisson
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?
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. 0
broom.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 observan diferencias entre los modelos, porque, en el modelo
lineal da un estimate del intercepto de 3.85 y cuando se aleatorizan la
variable individuo en el modelo mixtose observaun ligero cambio de ese
intercepto a 3.83.
Por otro lado, se observa que la identidad del pollo capta alguna
variabilidad y esto genera que no sean iguales los modelos.
Para concluir si se toma en cuenta la variablidad entre individuos
esto mejora el modelo, asimismo, cuando se toma la indentidad del pollo
es necesario que se tenga más de una medida por cada individuo para
poder hacer el modelo aleatorio.
CONCLUSIÓN:
Se usa el modelo lineal mixto porque se observa que la identidad de
los individuos genera un cambio en los resultados, por esta razón, este
modelo permite dividir entre variables fijas (controlados) y aleatorias
(no controlados). En cambio, el modelo lineal solo toma en cuenta un
efecto.
Ventajas en el uso del modelo lineal mixto en la
agricultura
1) Con este modelo se pueden establecer datos con estrucutra
jerarquica para poder observar niveles, que pueden ser parcelas,
regiones, lotes, etc y así poder determinar los efectos fijos y los
efectos aleatorios de cada nivel, por otro lado, este modelo permite
modelar la variabilidad de otros factores no explicados dentro de
nuestro modelo.
2) Es utilizado en programas de mejoramiento genético enfocado en
tres direcciones: predicción del comportamiento familiar, estimación de
componentes de varianza y ensayos multiambientales.