Entrega Final DE Modelos lineales mixtos - variables aleatorias
Andrés Felipe Hincapié Castañeda
23/06/2023
Modelos lineales mixtos - variables aleatorias
# Librerias
library(tidyverse)
library(broom.mixed)
library(lme4)# Cargar Base de Datos
data("CO2")
# Vista Gráfica de los Datos
ggplot(CO2,
aes(x = conc, y = uptake, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment))
En términos generales las plantas de Quebec fueron capaces de fijar mayores cantidades de nitrógeno para todas las concentraciones de CO2. Hay una aparente relación directamente proporcional entre la capacidad de fijación y la cantidad de CO2 inicial.
La gráfica sugiere que la variabilidad está asociada al tipo de planta.
Para poder utilizar este modelo con plantas que no estén dentro de este experimento y adaptar esa variable a aleatorio vamos a hacer uso de lme4
# Para modelo lineal simple
fit1 <- lm(uptake ~ conc + Type + Treatment, data = CO2)
# Para agregar el Factor aleatorio
fit2 <- lmer(uptake ~ conc + Type + Treatment + (1 | Plant), data = CO2)Cada una de las plantas va a tener un intercepto distinto
# Para modelo lineal Simple
broom::glance(fit1)# Para modelo lineal Simple
broom::tidy(fit1)# Para modelo lineal MIXTO
broom.mixed::glance(fit2)# Para modelo lineal MIXTO
broom.mixed::tidy(fit2)GLM PARA EFECTOS FIJOS
Agregar una familia
# Cargar Datos
# 2 Factores Tiempo y Dieta
# Identidad de Pollo (50 pollos diferentes)
# Tiempo: 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 21
# Dieta: 1, 2, 3, 4
# Variable Respuesta: Peso
data("ChickWeight")# Vista Gráfica de los datos
ggplot(ChickWeight,
aes(x = Time, y = weight)) +
geom_point(aes(color = Diet)) +
geom_path(aes(color = Diet, group = Chick))
De la gráfica se puede apreciar hay pollos iguales que tienen la misma dieta pero los valores del peso difieren mucho entre sí.
La variable respuesta está dada por un número entero (sin decimales) positivo que va desde 35 hasta 373; por lo que puede ser pensado como un gml de tipo poisson.
# Definir que el peso es el resultado de la interacción entre el tiempo y la dieta
fit1_poisson <- glm(weight ~ Diet:Time, data = ChickWeight, family = poisson())
# Para agregar el Factor aleatorio Chick
fit2_poisson <- glmer(weight ~ Diet:Time + (1|Chick), data = ChickWeight, family = poisson())broom::glance(fit1_poisson)broom::tidy(fit1_poisson)Existe un aumento del peso con el paso del tiempo y la interacción del tiempo con la dieta es diferente para cada una de ellas.
broom.mixed::glance(fit2_poisson)broom.mixed::tidy(fit2_poisson)NALa identidad de los Pollos estaba generando alguna variabilidad sobre la variable respuesta, es por esto que se obtienen diferencias entre el estimado del intercepto en el fit1(glm) y en el fit2 (glmer pollo aleatoreo)
La dieta que mejor está dando resultados es la dieta 3, con una pendiente mayor que las otras (0.086). Es importante mencionar que el valor de las pendientes SI cambia con respecto al realizado sin aleatorizar la identidad del Pollo.
Conclusión:
En este tipo de modelos donde se tiene la identidad de los individuos experimentados es importante reconocer la capacidad que tienen de alterar los valores de la variable respuesta, por eso es muy recomendado realizar un tratamiento estadistico y aleatorizar dichos valores para asegurarse que el analisis estadistico éste tomando únicamente las interferencias causadas por los tratamientos a analizar o de interes.
Los modelos mixtos, con un “factor” aleatorizado solo funcionan cuando para cada identidad del individuo hay varias medidas.
Ruido para diferenciar el grupo en el Primer Ejemplo Uptake de CO2
ruido <- runif(84, 0.02, 0.04)
ruido [1] 0.03864546 0.03238910 0.03973308 0.03303593 0.03466091 0.02222489 0.03228484 0.03627802
[9] 0.03095276 0.03382782 0.02659018 0.02265830 0.03899295 0.03421237 0.02642007 0.02767638
[17] 0.02752098 0.02847850 0.02033391 0.03378538 0.03900894 0.02362247 0.03108338 0.02162671
[25] 0.03612366 0.03520916 0.03170573 0.02978904 0.02269815 0.03935512 0.02413953 0.02368893
[33] 0.02580169 0.03313340 0.02214627 0.03802851 0.02782775 0.03084036 0.03911506 0.03877985
[41] 0.02895324 0.03962587 0.03369271 0.02322698 0.02103089 0.03157449 0.02970703 0.03124269
[49] 0.03175258 0.03359871 0.02751353 0.02788774 0.02252121 0.03384963 0.03408366 0.02685025
[57] 0.02233944 0.02168466 0.03905663 0.03800303 0.03468761 0.03064149 0.03157105 0.02483939
[65] 0.02535009 0.02986647 0.03294344 0.02786299 0.02929221 0.03280463 0.03150133 0.02287919
[73] 0.02065500 0.03512132 0.03785358 0.03233157 0.03157752 0.02147002 0.02990935 0.02021458
[81] 0.02153992 0.02527841 0.02210986 0.02302687CO2
Grouped Data: uptake ~ conc | Plant
Plant Type Treatment conc uptake uptaker
1 Qn1 Quebec nonchilled 95 16.0 16.038645
2 Qn1 Quebec nonchilled 175 30.4 30.432389
3 Qn1 Quebec nonchilled 250 34.8 34.839733
4 Qn1 Quebec nonchilled 350 37.2 37.233036
5 Qn1 Quebec nonchilled 500 35.3 35.334661
6 Qn1 Quebec nonchilled 675 39.2 39.222225
7 Qn1 Quebec nonchilled 1000 39.7 39.732285
8 Qn2 Quebec nonchilled 95 13.6 13.636278
9 Qn2 Quebec nonchilled 175 27.3 27.330953
10 Qn2 Quebec nonchilled 250 37.1 37.133828
11 Qn2 Quebec nonchilled 350 41.8 41.826590
12 Qn2 Quebec nonchilled 500 40.6 40.622658
13 Qn2 Quebec nonchilled 675 41.4 41.438993
14 Qn2 Quebec nonchilled 1000 44.3 44.334212
15 Qn3 Quebec nonchilled 95 16.2 16.226420
16 Qn3 Quebec nonchilled 175 32.4 32.427676
17 Qn3 Quebec nonchilled 250 40.3 40.327521
18 Qn3 Quebec nonchilled 350 42.1 42.128479
19 Qn3 Quebec nonchilled 500 42.9 42.920334
20 Qn3 Quebec nonchilled 675 43.9 43.933785
21 Qn3 Quebec nonchilled 1000 45.5 45.539009
22 Qc1 Quebec chilled 95 14.2 14.223622
23 Qc1 Quebec chilled 175 24.1 24.131083
24 Qc1 Quebec chilled 250 30.3 30.321627
25 Qc1 Quebec chilled 350 34.6 34.636124
26 Qc1 Quebec chilled 500 32.5 32.535209
27 Qc1 Quebec chilled 675 35.4 35.431706
28 Qc1 Quebec chilled 1000 38.7 38.729789
29 Qc2 Quebec chilled 95 9.3 9.322698
30 Qc2 Quebec chilled 175 27.3 27.339355
31 Qc2 Quebec chilled 250 35.0 35.024140
32 Qc2 Quebec chilled 350 38.8 38.823689
33 Qc2 Quebec chilled 500 38.6 38.625802
34 Qc2 Quebec chilled 675 37.5 37.533133
35 Qc2 Quebec chilled 1000 42.4 42.422146
36 Qc3 Quebec chilled 95 15.1 15.138029
37 Qc3 Quebec chilled 175 21.0 21.027828
38 Qc3 Quebec chilled 250 38.1 38.130840
39 Qc3 Quebec chilled 350 34.0 34.039115
40 Qc3 Quebec chilled 500 38.9 38.938780
41 Qc3 Quebec chilled 675 39.6 39.628953
42 Qc3 Quebec chilled 1000 41.4 41.439626
43 Mn1 Mississippi nonchilled 95 10.6 10.633693
44 Mn1 Mississippi nonchilled 175 19.2 19.223227
45 Mn1 Mississippi nonchilled 250 26.2 26.221031
46 Mn1 Mississippi nonchilled 350 30.0 30.031574
47 Mn1 Mississippi nonchilled 500 30.9 30.929707
48 Mn1 Mississippi nonchilled 675 32.4 32.431243
49 Mn1 Mississippi nonchilled 1000 35.5 35.531753
50 Mn2 Mississippi nonchilled 95 12.0 12.033599
51 Mn2 Mississippi nonchilled 175 22.0 22.027514
52 Mn2 Mississippi nonchilled 250 30.6 30.627888
53 Mn2 Mississippi nonchilled 350 31.8 31.822521
54 Mn2 Mississippi nonchilled 500 32.4 32.433850
55 Mn2 Mississippi nonchilled 675 31.1 31.134084
56 Mn2 Mississippi nonchilled 1000 31.5 31.526850
57 Mn3 Mississippi nonchilled 95 11.3 11.322339
58 Mn3 Mississippi nonchilled 175 19.4 19.421685
59 Mn3 Mississippi nonchilled 250 25.8 25.839057
60 Mn3 Mississippi nonchilled 350 27.9 27.938003
61 Mn3 Mississippi nonchilled 500 28.5 28.534688
62 Mn3 Mississippi nonchilled 675 28.1 28.130641
63 Mn3 Mississippi nonchilled 1000 27.8 27.831571
64 Mc1 Mississippi chilled 95 10.5 10.524839
65 Mc1 Mississippi chilled 175 14.9 14.925350
66 Mc1 Mississippi chilled 250 18.1 18.129866
67 Mc1 Mississippi chilled 350 18.9 18.932943
68 Mc1 Mississippi chilled 500 19.5 19.527863
69 Mc1 Mississippi chilled 675 22.2 22.229292
70 Mc1 Mississippi chilled 1000 21.9 21.932805
71 Mc2 Mississippi chilled 95 7.7 7.731501
72 Mc2 Mississippi chilled 175 11.4 11.422879
73 Mc2 Mississippi chilled 250 12.3 12.320655
74 Mc2 Mississippi chilled 350 13.0 13.035121
75 Mc2 Mississippi chilled 500 12.5 12.537854
76 Mc2 Mississippi chilled 675 13.7 13.732332
77 Mc2 Mississippi chilled 1000 14.4 14.431578
78 Mc3 Mississippi chilled 95 10.6 10.621470
79 Mc3 Mississippi chilled 175 18.0 18.029909
80 Mc3 Mississippi chilled 250 17.9 17.920215
81 Mc3 Mississippi chilled 350 17.9 17.921540
82 Mc3 Mississippi chilled 500 17.9 17.925278
83 Mc3 Mississippi chilled 675 18.9 18.922110
84 Mc3 Mississippi chilled 1000 19.9 19.923027# Vista Gráfica de los Datos
ggplot(CO2,
aes(x = conc, y = uptaker, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment))
El ruido generó cambios evidenciables en la distribución de CO2 fijado por las plantas, sin embargo aún se logra evidenciar que las plantas oriundas de Quebec tienen una mayor capacidad de fijación en comparación a las plantas de Mississipi. Adicionalmente, es posible ver como el tratamiento juega un papel importante en la capacidad de fijación, en la mayoría de las situaciones las plantas que fueron chilled o enfriadas presentaron una menor capacidad de fijación.
Como era de esperar la gráfica sugiere que la variabilidad está asociada al tipo de planta, por esta razón vamos a hacer uso de lme4, así adaptamos la variable planta a aleatorio y asegurar que dicha variable no genere efectos sobre la variable respuesta.
# Agregar el Factor aleatorio
fitsh <- lmer(uptaker ~ conc + Type + Treatment, data = CO2)
Error: No random effects terms specified in formula# Para modelo lineal MIXTO
broom.mixed::glance(fith)# Para modelo lineal MIXTO
broom.mixed::tidy(fith)Los valores para los coeficientes estimados provocados por la concentración de co2 (conc) y el origen de la planta (Mississippi o Quebec) son apreciables en la columna de estimate; éstos son los valores que nos permitirán usar el modelo en contextos diferentes a éste experimento, pues no están inflenciados por algo específico como la planta expermentada.
Podemos apreciar como el factor planta fue intercepatado y no se tiene en cuenta el coeficiente de variación que estaría causando de no ser así, en caso de no haber aleatorizado este factor, estaría alterando nuestra variable respuesta a razón de 1.708.
Ésta tabla también nos permite confirmar lo propuesto en el analisis gráfico, “TypeMississippi” muestra una diferencia esperada de -12.66185859 unidades en la variable de respuesta, y “Treatmentchilled” muestra una diferencia esperada de -6.86105965. Lo que explica que las plantas de Mississippi y las tratadas con chilled tengan una capacidad de fijación menor.
Conclusión Final
Un modelo de efectos mixtos es de gran utilidad para los analisis estadisticos ya que permite establecer algunas variables como fijas y evitar que éstas tengan algún efecto sobre la variable respuesta a la hora de hacer el analisis estadistico y una potencial predicción de datos; como lo explicaban en el video, en algunas ocasiones los experimentos pueden ser muy especificos con las variables (como la identidad de los pollos), pero éstas variables no ser relevantes o de interes para conocer los efectos de los tratamientos que se están estudiando.
En agricultura sabemos que hay muchas variables que pueden afectar el rendimiento de un cultivo, los modelos mixtos son una herramienta muy útil ya que permiten segregar o interceptar las variaciones causadas por variables no deseadas, así capturar la variación no explicada y establecer una correcta relación entre la variable respuesta y los tratamientos de interes.
Esto conduce a un análisis más preciso y a una mejor comprensión de los factores que influyen en los resultados de los experimentos agronomicos.