Random Effects Models
Erik David Aldana
26/8/2021
Modelo efectos Fijos
library(readxl)
Random_Effects_Models <- read_excel("C:/Users/ERICK/Desktop/Semestre 2021-I/Disenio de Experimentos/Ejercicios/Random Effects Models.xlsx",
sheet = "Modelo de Efectos Fijos")
Rend= Random_Effects_Models$Rendimiento
Gen= Random_Effects_Models$Genotipos; as.factor(Gen)
Amb= Random_Effects_Models$Ambientes; as.factor(Amb)
Rep=Random_Effects_Models$Repeticiones; as.factor(Rep)
df<-Random_Effects_Models;dflibrary(collapsibleTree)
collapsibleTreeSummary(Random_Effects_Models, hierarchy=c ("Ambientes", "Genotipos", "Repeticiones", "Rendimiento"))boxplot( Rend~Gen,las=2, df)
boxplot( Rend~Amb,las=2, df)
medias<-tapply(Rend, list(Gen, Amb), mean)
Desv<-tapply(Rend, list(Gen, Amb), sd)
cv=100*Desv/mediasAnálisis de Varianza
Mod1<-aov(Rend~Gen*Amb,df)
summary(Mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## Gen 29 39289 1355 0.838 0.709
## Amb 9 6084 676 0.418 0.925
## Gen:Amb 261 409660 1570 0.970 0.598
## Residuals 300 485283 1618Revisión de Supuestos
#Comparación de medias
library(agricolae)
HSD_test = with(df, HSD.test (Rend, Gen, DFerror = 300,MSerror = 1618))
HSD_test## $statistics
## MSerror Df Mean CV MSD
## 1618 300 286.9354 14.01861 48.15775
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey Gen 30 5.354168 0.05
##
## $means
## Rend std r Min Max Q25 Q50 Q75
## Col-11 276.0715 29.82083 20 201.6180 338.9759 257.5360 280.5669 292.7454
## Col-13 276.2776 27.47255 20 230.5845 328.2326 261.1193 274.5290 287.1297
## Col-14 292.6078 31.47335 20 238.2013 372.4844 268.6544 286.9483 308.5266
## Col-15 300.0565 43.19628 20 249.3834 383.7265 262.6716 285.3840 338.8865
## Col-16 289.1938 52.15930 20 189.1805 418.2087 253.3516 288.8335 318.4868
## Col-17 276.9876 43.64813 20 170.7534 347.4957 258.8673 276.3349 300.7460
## Col-19 285.6930 44.55251 20 194.5887 355.8221 267.7497 281.2020 324.5590
## Col-20 286.2532 44.28823 20 222.1054 406.8233 257.1418 277.0342 314.3994
## Col-21 282.3140 44.02041 20 211.4012 398.1853 253.3804 278.2097 308.0645
## Col-23 290.9390 42.96054 20 180.0171 375.0906 274.2403 291.0797 306.3484
## Col-27 284.5826 41.03584 20 212.9674 356.3278 259.9157 287.2120 320.1600
## Col-3 289.1546 46.59129 20 193.4529 400.2604 264.6375 279.6279 308.6763
## Col-30 273.8719 36.28836 20 212.3708 365.2514 252.1911 262.0110 284.5737
## Col-31 283.2879 36.97459 20 220.4327 365.3183 260.1640 276.1711 310.2211
## Col-32 283.4009 37.29838 20 218.1608 366.4654 262.8878 284.9818 299.8092
## Col-33 292.0710 43.95785 20 222.6497 356.4889 252.4947 289.8415 331.6462
## Col-34 295.7183 39.35450 20 200.6801 361.8237 278.3977 305.0901 322.4632
## Col-35 288.4559 46.48229 20 205.9930 378.4565 249.8825 295.8397 319.8337
## Col-36 268.0548 26.75652 20 213.4396 317.5826 250.8101 264.6241 284.0954
## Col-37 299.0274 42.05575 20 240.5460 376.7750 259.5521 298.2148 330.9174
## Col-38 290.1597 42.76571 20 162.0717 343.0395 268.1100 294.5173 321.9192
## Col-4 286.4697 36.79078 20 197.3375 369.8054 267.5245 286.7927 308.0851
## Col-40 292.5114 38.90496 20 196.5161 355.9223 269.4490 289.8617 315.9759
## Col-41 285.7607 35.76239 20 191.5233 340.8131 265.3794 290.4587 310.0349
## Col-42 301.8878 36.79634 20 215.3350 384.5530 283.8936 301.4930 322.5690
## Col-5 284.6546 33.35554 20 221.0930 348.9119 260.7482 283.0406 306.7781
## Col-6 294.1657 30.10734 20 239.9994 348.9012 273.2686 288.3108 317.0407
## Col-7 284.3262 36.03634 20 203.4413 360.5467 260.0466 287.3190 308.8660
## Col-8 276.1299 43.95701 20 193.7119 349.5078 243.6677 271.2431 306.2384
## Col-9 297.9775 43.94826 20 215.0179 357.3131 272.9975 304.5212 334.2741
##
## $comparison
## NULL
##
## $groups
## Rend groups
## Col-42 301.8878 a
## Col-15 300.0565 a
## Col-37 299.0274 a
## Col-9 297.9775 a
## Col-34 295.7183 a
## Col-6 294.1657 a
## Col-14 292.6078 a
## Col-40 292.5114 a
## Col-33 292.0710 a
## Col-23 290.9390 a
## Col-38 290.1597 a
## Col-16 289.1938 a
## Col-3 289.1546 a
## Col-35 288.4559 a
## Col-4 286.4697 a
## Col-20 286.2532 a
## Col-41 285.7607 a
## Col-19 285.6930 a
## Col-5 284.6546 a
## Col-27 284.5826 a
## Col-7 284.3262 a
## Col-32 283.4009 a
## Col-31 283.2879 a
## Col-21 282.3140 a
## Col-17 276.9876 a
## Col-13 276.2776 a
## Col-8 276.1299 a
## Col-11 276.0715 a
## Col-30 273.8719 a
## Col-36 268.0548 a
##
## attr(,"class")
## [1] "group"#Normalidad de residuales
res<-Mod1$residuals
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.99784, p-value = 0.6474plot(Mod1,2)
Modelo efectos Aleatorios
Random_Effects_Models2 <- read_excel("C:/Users/ERICK/Desktop/Semestre 2021-I/Disenio de Experimentos/Ejercicios/Random Effects Models.xlsx",
sheet = "Modelo de Efectos Aleatorios")
Rend2= Random_Effects_Models2$Rendimiento
Gen2= Random_Effects_Models2$Genotipos; as.factor(Gen2)
Amb2= Random_Effects_Models2$Ambientes; as.factor(Amb2)
Rep2=Random_Effects_Models2$Repeticiones; as.factor(Rep2)
df2<-Random_Effects_Models2; df2collapsibleTreeSummary(Random_Effects_Models2, hierarchy=c ("Ambientes", "Genotipos", "Repeticiones", "Rendimiento"))library(lattice)
boxplot( Rend2~Gen2,las=3,font = 4,font.lab = 4,font.axis = 3, df2)
boxplot( Rend2~Amb2,las=3,font = 4,font.lab = 4,font.axis = 3, df2)
bwplot(Rend2~Gen2|Amb2,las=2, df2)
library(lme4)## Loading required package: MatrixMod2 <- lmer(Rend2 ~ (1 | Gen2) + (1 | Amb2) + (1 | Amb2:Gen2), data = df2)## boundary (singular) fit: see ?isSingularsummary(Mod2)## Linear mixed model fit by REML ['lmerMod']
## Formula: Rend2 ~ (1 | Gen2) + (1 | Amb2) + (1 | Amb2:Gen2)
## Data: df2
##
## REML criterion at convergence: 391.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.46201 -0.75337 -0.00318 0.77977 1.81198
##
## Random effects:
## Groups Name Variance Std.Dev.
## Amb2:Gen2 (Intercept) 0 0.00
## Gen2 (Intercept) 0 0.00
## Amb2 (Intercept) 0 0.00
## Residual 1213 34.82
## Number of obs: 40, groups: Amb2:Gen2, 20; Gen2, 5; Amb2, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 286.412 5.506 52.02
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularvar_amb_gen2 = 0.00**2
var_amb2 = 0.00**2
var_gen2 = 0.00**2
var_res2 = 34.82**2
var_tot2 = var_amb_gen2 + var_amb2 + var_gen2 + var_res2; var_tot2## [1] 1212.432100 * var_amb2/var_tot2## [1] 0100 * var_gen2/var_tot2## [1] 0library(agricolae)
HSD_test = with(df2, HSD.test (Rend2, Gen2, DFerror = 300,MSerror = 202.4))
HSD_test## $statistics
## MSerror Df Mean CV MSD
## 202.4 300 286.4116 4.967234 19.52222
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey Gen2 5 3.881226 0.05
##
## $means
## Rend2 std r Min Max Q25 Q50 Q75
## Col-13 283.3235 30.99016 8 236.3346 328.2326 268.1952 281.2495 301.3427
## Col-21 285.4014 30.38124 8 243.5142 324.8634 262.1829 282.4404 310.7889
## Col-34 285.7805 47.38611 8 200.6801 330.8739 267.4040 306.4182 315.3705
## Col-40 294.1548 28.12727 8 257.6374 341.5524 278.6421 291.9509 311.8766
## Col-8 283.3980 41.76410 8 228.0917 349.5078 255.3630 271.2431 315.9077
##
## $comparison
## NULL
##
## $groups
## Rend2 groups
## Col-40 294.1548 a
## Col-34 285.7805 a
## Col-21 285.4014 a
## Col-8 283.3980 a
## Col-13 283.3235 a
##
## attr(,"class")
## [1] "group"## Plot data ####
with(df2, interaction.plot(x.factor = Gen2, trace.factor = Amb2, response = Rend2))
Modelo efectos Mixtos
Random_Effects_Models3 <- read_excel("C:/Users/ERICK/Desktop/Semestre 2021-I/Disenio de Experimentos/Ejercicios/Random Effects Models.xlsx",
sheet = "Modelo de Efectos Mixtos")
Rend3= Random_Effects_Models3$Rendimiento
Gen3= Random_Effects_Models3$Genotipos; as.factor(Gen3)
Amb3= Random_Effects_Models3$Ambientes; as.factor(Amb3)
Rep3=Random_Effects_Models3$Repeticiones; as.factor(Rep3)
df3<-Random_Effects_Models2; df3collapsibleTreeSummary(Random_Effects_Models3, hierarchy=c ("Ambientes", "Genotipos", "Repeticiones", "Rendimiento"))boxplot( Rend3~Gen3,col.lab = 3,las= 3, df3)
boxplot( Rend3~Amb3,col.lab = 4,las= 3, df3)
bwplot(Rend3~Gen3|Amb3,las=3, df3)
Mod3 <- lmer(Rend3 ~ Amb3 + (1 | Gen3) + (1 | Gen3:Amb3), data = df3)## boundary (singular) fit: see ?isSingularsummary(Mod3)## Linear mixed model fit by REML ['lmerMod']
## Formula: Rend3 ~ Amb3 + (1 | Gen3) + (1 | Gen3:Amb3)
## Data: df3
##
## REML criterion at convergence: 944.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4274 -0.6861 0.0620 0.6384 2.8448
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gen3:Amb3 (Intercept) 0.000 0.000
## Gen3 (Intercept) 1.543 1.242
## Residual 1631.325 40.390
## Number of obs: 100, groups: Gen3:Amb3, 50; Gen3, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 277.780 12.784 21.728
## Amb3A10 22.342 18.063 1.237
## Amb3A2 14.133 18.063 0.782
## Amb3A3 6.200 18.063 0.343
## Amb3A4 3.851 18.063 0.213
## Amb3A5 -3.947 18.063 -0.219
## Amb3A6 6.716 18.063 0.372
## Amb3A7 5.549 18.063 0.307
## Amb3A8 -0.730 18.063 -0.040
## Amb3A9 13.993 18.063 0.775
##
## Correlation of Fixed Effects:
## (Intr) Am3A10 Amb3A2 Amb3A3 Amb3A4 Amb3A5 Amb3A6 Amb3A7 Amb3A8
## Amb3A10 -0.706
## Amb3A2 -0.706 0.500
## Amb3A3 -0.706 0.500 0.500
## Amb3A4 -0.706 0.500 0.500 0.500
## Amb3A5 -0.706 0.500 0.500 0.500 0.500
## Amb3A6 -0.706 0.500 0.500 0.500 0.500 0.500
## Amb3A7 -0.706 0.500 0.500 0.500 0.500 0.500 0.500
## Amb3A8 -0.706 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## Amb3A9 -0.706 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularvar_amb_gen3 = 0.000**2
var_amb3 = 0.000**2
var_gen3 = 1.242**2
var_res3 = 40.390**2
var_tot3 = var_amb_gen3 + var_amb3 + var_gen3 + var_res3; var_tot3## [1] 1632.895100 * var_amb3/var_tot3## [1] 0100 * var_gen3/var_tot3## [1] 0.09446807