data(edwards.oats)
dat <- edwards.oats
str(dat)
## 'data.frame': 3694 obs. of 7 variables:
## $ eid : int 1 1 1 1 1 1 1 1 1 1 ...
## $ year : int 1997 1997 1997 1997 1997 1997 1997 1997 1997 1997 ...
## $ loc : Factor w/ 5 levels "Ame","Cra","Lew",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ block : Factor w/ 3 levels "B1","B2","B3": 1 1 1 1 1 1 1 1 1 1 ...
## $ gen : Factor w/ 80 levels "ACStewart","Belle",..: 2 4 6 7 8 9 10 11 13 14 ...
## $ yield : num 151 172 159 152 156 ...
## $ testwt: num 32 33.3 33.7 32.7 34.4 34.6 33.9 36.9 33.1 34.3 ...
dat$year <- as.factor(dat$year)
head(dat)
## eid year loc block gen yield testwt
## 1 1 1997 Ame B1 Belle 151.1 32.0
## 2 1 1997 Ame B1 Brawn 171.9 33.3
## 3 1 1997 Ame B1 Burton 158.9 33.7
## 4 1 1997 Ame B1 Chairman 151.7 32.7
## 5 1 1997 Ame B1 Chaps 155.6 34.4
## 6 1 1997 Ame B1 Classic 153.8 34.6
model <- lmer(yield ~ gen +
(1|block:loc) +
(1|loc:year) +
(1|gen:loc:year), data = dat)
summary(model)
## Linear mixed model fit by REML ['lmerMod']
## Formula: yield ~ gen + (1 | block:loc) + (1 | loc:year) + (1 | gen:loc:year)
## Data: dat
##
## REML criterion at convergence: 29079.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4692 -0.5328 0.0155 0.5502 3.8282
##
## Random effects:
## Groups Name Variance Std.Dev.
## gen:loc:year (Intercept) 54.365 7.373
## loc:year (Intercept) 793.383 28.167
## block:loc (Intercept) 6.216 2.493
## Residual 121.744 11.034
## Number of obs: 3694, groups: gen:loc:year, 1235; loc:year, 34; block:loc, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 119.6767 7.7884 15.366
## genBelle -6.9627 6.3009 -1.105
## genBlaze 10.5086 6.3009 1.668
## genBrawn 10.7049 6.3009 1.699
## genBuff -30.0823 7.5414 -3.989
## genBurton -5.2953 6.4636 -0.819
## genChairman -5.0141 7.4325 -0.675
## genChaps 8.9295 6.3009 1.417
## genClassic 6.3055 6.3009 1.001
## genDane -0.6115 6.3009 -0.097
## genDon -5.3498 6.3009 -0.849
## genEbeltoft 5.7662 6.4779 0.890
## genGem 1.9641 6.3009 0.312
## genHazel -16.4270 8.2378 -1.994
## genHyproBlend -18.2222 8.4685 -2.152
## genIA91462-4-1-6 -3.2949 6.6018 -0.499
## genIA91483-7-18-2 -8.5516 7.5237 -1.137
## genIA93227-2-1 2.8359 7.5414 0.376
## genIA93359-3 6.9674 6.8485 1.017
## genIAACE2-16-1 6.3702 7.5414 0.845
## genIAH61-3-3-4 3.1297 7.5237 0.416
## genIAK993-7-5 -24.0677 11.2543 -2.139
## genIAL24-5-1 -8.3777 8.2198 -1.019
## genIAN979-5-1-22 -1.7646 6.5940 -0.268
## genIAR56-5 5.9555 6.4985 0.916
## genIAR56-8 -0.6991 7.8412 -0.089
## genIAR66-6 12.0140 6.9349 1.732
## genIL95-1241 20.8310 6.8524 3.040
## genIN09201 7.2628 6.3020 1.152
## genIda -0.3243 6.6338 -0.049
## genJay 8.4748 6.2931 1.347
## genJerry 0.8321 6.3009 0.132
## genJim 6.6048 6.3020 1.048
## genJud 4.0873 6.3009 0.649
## genKilldeer 8.9962 6.4028 1.405
## genLeonard 7.9280 6.8524 1.157
## genLoyal -4.0575 6.6018 -0.615
## genMN 93272 -6.0007 7.4325 -0.807
## genMN98236 10.2224 7.5414 1.356
## genMO1022-1 4.0502 6.8485 0.591
## genMO8701 0.2635 6.8485 0.038
## genMilton -1.6920 6.4636 -0.262
## genMoraine -2.3979 6.4028 -0.375
## genND9816784 -6.0029 7.5414 -0.796
## genNewdak -0.3296 7.4325 -0.044
## genOA1021-1 14.5421 7.5414 1.928
## genOA21-2 9.3174 7.5296 1.237
## genOH1086 -4.0418 7.8412 -0.515
## genOH1087 -5.4430 7.4325 -0.732
## genOH1128 -6.5585 7.8412 -0.836
## genOH1149 -12.8409 6.8628 -1.871
## genOgle 5.7586 6.3009 0.914
## genP9213A2-2-4-4-3 7.5660 6.8579 1.103
## genP971A9-7-4-1 15.9820 7.5414 2.119
## genP973A38-3-6 6.2070 7.5296 0.824
## genP976A6-9-6 16.7725 7.5296 2.228
## genPal -25.0325 7.4325 -3.368
## genPrairie 8.1847 6.8628 1.193
## genPremier -11.7215 7.4325 -1.577
## genReeves -3.8906 6.6075 -0.589
## genRichard -5.7073 6.3584 -0.898
## genRiser -15.3922 6.3009 -2.443
## genRodeo 8.9898 6.3365 1.419
## genSD0000731 -3.6779 7.5414 -0.488
## genSD96024 12.0065 6.6075 1.817
## genSesqui 8.3360 6.4804 1.286
## genSheldon -5.5488 6.3365 -0.876
## genStarter -14.2558 6.3009 -2.262
## genStarter/Amag 2.1059 7.8412 0.269
## genTroy -5.6983 6.3009 -0.904
## genValley -5.4289 6.4636 -0.840
## genVista -3.9432 6.3365 -0.622
## genWIX6356-1 -3.7290 6.8628 -0.543
## genWIX7541-3 6.8116 7.5272 0.905
## genWIX7822-3 10.7024 7.5414 1.419
## genWIX8177-1 5.2981 7.5414 0.703
## genWIX8179-2 21.5484 7.5414 2.857
## genWabasha 2.0297 6.4837 0.313
## genWebster -0.4745 8.2378 -0.058
## genYoungs 0.1839 6.4779 0.028
##
## Correlation matrix not shown by default, as p = 80 > 12.
## Use print(value, correlation=TRUE) or
## vcov(value) if you need it
print(model, correlation=TRUE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: yield ~ gen + (1 | block:loc) + (1 | loc:year) + (1 | gen:loc:year)
## Data: dat
## REML criterion at convergence: 29079.87
## Random effects:
## Groups Name Std.Dev.
## gen:loc:year (Intercept) 7.373
## loc:year (Intercept) 28.167
## block:loc (Intercept) 2.493
## Residual 11.034
## Number of obs: 3694, groups: gen:loc:year, 1235; loc:year, 34; block:loc, 15
## Fixed Effects:
## (Intercept) genBelle genBlaze genBrawn
## 119.6767 -6.9627 10.5086 10.7049
## genBuff genBurton genChairman genChaps
## -30.0823 -5.2953 -5.0141 8.9295
## genClassic genDane genDon genEbeltoft
## 6.3055 -0.6115 -5.3498 5.7662
## genGem genHazel genHyproBlend genIA91462-4-1-6
## 1.9641 -16.4270 -18.2222 -3.2949
## genIA91483-7-18-2 genIA93227-2-1 genIA93359-3 genIAACE2-16-1
## -8.5516 2.8359 6.9674 6.3702
## genIAH61-3-3-4 genIAK993-7-5 genIAL24-5-1 genIAN979-5-1-22
## 3.1297 -24.0677 -8.3777 -1.7646
## genIAR56-5 genIAR56-8 genIAR66-6 genIL95-1241
## 5.9555 -0.6991 12.0140 20.8310
## genIN09201 genIda genJay genJerry
## 7.2628 -0.3243 8.4748 0.8321
## genJim genJud genKilldeer genLeonard
## 6.6048 4.0873 8.9962 7.9280
## genLoyal genMN 93272 genMN98236 genMO1022-1
## -4.0575 -6.0007 10.2224 4.0502
## genMO8701 genMilton genMoraine genND9816784
## 0.2635 -1.6920 -2.3979 -6.0029
## genNewdak genOA1021-1 genOA21-2 genOH1086
## -0.3296 14.5421 9.3174 -4.0418
## genOH1087 genOH1128 genOH1149 genOgle
## -5.4430 -6.5585 -12.8409 5.7586
## genP9213A2-2-4-4-3 genP971A9-7-4-1 genP973A38-3-6 genP976A6-9-6
## 7.5660 15.9820 6.2070 16.7725
## genPal genPrairie genPremier genReeves
## -25.0325 8.1847 -11.7215 -3.8906
## genRichard genRiser genRodeo genSD0000731
## -5.7073 -15.3922 8.9898 -3.6779
## genSD96024 genSesqui genSheldon genStarter
## 12.0065 8.3360 -5.5488 -14.2558
## genStarter/Amag genTroy genValley genVista
## 2.1059 -5.6983 -5.4289 -3.9432
## genWIX6356-1 genWIX7541-3 genWIX7822-3 genWIX8177-1
## -3.7290 6.8116 10.7024 5.2981
## genWIX8179-2 genWabasha genWebster genYoungs
## 21.5484 2.0297 -0.4745 0.1839
# Residual vs Fitted plot
residuals <- resid(model)
fitted_values <- fitted(model)
res_fit_plot <- ggplot(data = data.frame(fitted = fitted_values,
residuals = residuals),
aes(x = fitted, y = residuals)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(title = "Residual vs Fitted Plot", x = "Fitted Values", y = "Residuals") +
theme_minimal()
res_fit_plot
# Extract random effects from model
ranef_model <- ranef(model)$loc
# Plotting random effects
ranef_plot <- ggplot(data = data.frame(location = rownames(ranef_model),
effect = ranef_model[,1]),
aes(x = location, y = effect)) +
geom_bar(stat = "identity") +
labs(title = "Random Effects of Location By Year", x = "Location x Year", y = "Random Effect") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
ranef_plot
# Predicted values (fitted values)
predicted_values <- predict(model)
pred_val <- ggplot(dat, aes(x = loc, y = yield)) +
geom_point(aes(color = year), size = 3) +
geom_line(aes(y = predicted_values, group = gen),
linetype = "dashed", color = "blue") +
labs(title = "Observed vs Predicted yield by Genotype",
x = "Location", y = "Number Kernel Ears") +
theme_minimal()
pred_val
plt_mdl11 <- plot_model(model, type = "diag") # For fixed effects "diag"
plt_mdl11
## [[1]]
## `geom_smooth()` using formula = 'y ~ x'
##
## [[2]]
## [[2]]$`gen:loc:year`
## `geom_smooth()` using formula = 'y ~ x'
##
## [[2]]$`loc:year`
## `geom_smooth()` using formula = 'y ~ x'
##
## [[2]]$`block:loc`
## `geom_smooth()` using formula = 'y ~ x'
##
##
## [[3]]
##
## [[4]]
## `geom_smooth()` using formula = 'y ~ x'
plt_mdl <- plot_model(model, type = "re") # For random effects
plt_mdl
## [[1]]
##
## [[2]]
##
## [[3]]
