lu.stability.RdMulti-environment trial to illustrate stability statistics
data("lu.stability")
A data frame with 120 observations on the following 4 variables.
yieldyield
gengenotype factor, 5 levels
envenvironment factor, 6 levels
blockblock factor, 4 levels
Data for 5 maize genotypes in 2 years x 3 sites = 6 environments.
H.Y. Lu and C. T. Tien. (1993) Studies on nonparametric method of phenotypic stability: II. Selection for stability of agroeconomic concept. J. Agric. Assoc. China 164:1-17.
Hsiu Ying Lu. 1995. PC-SAS Program for Estimating Huehn's Nonparametric Stability Statistics. Agron J. 87:888-891.
Kae-Kang Hwu and Li-yu D Liu. (2013) Stability Analysis Using Multiple Environment Trials Data by Linear Regression. (In Chinese) Crop, Environment & Bioinformatics 10:131-142.
# \dontrun{ library(agridat) data(lu.stability) dat <- lu.stability # GxE means. Match Lu 1995 table 1 libs(reshape2) datm <- acast(dat, gen~env, fun=mean, value.var='yield') round(datm, 2)#> E1 E2 E3 E4 E5 E6 #> G1 3288.75 2162.50 3692.75 3920.75 1588.75 1805.25 #> G2 3487.25 1659.50 3148.25 3897.25 1936.50 1404.00 #> G3 2651.25 1996.25 2374.50 2919.75 982.50 1158.00 #> G4 2467.00 2095.25 2574.50 3743.00 1845.75 1488.75 #> G5 2752.50 1761.25 3912.75 4178.75 1889.50 1678.00#> G1 G2 G3 G4 G5 #> 2743.125 2588.792 2013.708 2369.042 2695.458#> E1 E2 E3 E4 E5 E6 #> 2929.35 1934.95 3140.55 3731.90 1648.60 1506.80# Traditional ANOVA. Match Hwu table 2 # F value for gen,env m1 = aov(yield~env+gen+Error(block:env+env:gen), data=dat)#> Warning: Error() model is singular#> #> Error: block:env #> Df Sum Sq Mean Sq F value Pr(>F) #> env 5 82817889 16563578 102.9 1.36e-12 *** #> Residuals 18 2897821 160990 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Error: env:gen #> Df Sum Sq Mean Sq F value Pr(>F) #> gen 4 8573084 2143271 4.436 0.00995 ** #> Residuals 20 9662873 483144 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Error: Within #> Df Sum Sq Mean Sq F value Pr(>F) #> Residuals 72 12912246 179337# F value for gen:env, block:env m2 <- aov(yield ~ gen + env + gen:env + block:env, data=dat) summary(m2)#> Df Sum Sq Mean Sq F value Pr(>F) #> gen 4 8573084 2143271 11.951 1.68e-07 *** #> env 5 82817889 16563578 92.360 < 2e-16 *** #> gen:env 20 9662873 483144 2.694 0.00114 ** #> env:block 18 2897821 160990 0.898 0.58302 #> Residuals 72 12912246 179337 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Finlay Wilkinson regression coefficients # First, calculate env mean, merge in libs(dplyr) dat2 <- group_by(dat, env) dat2 <- mutate(dat2, locmn=mean(yield)) m4 <- lm(yield ~ gen -1 + gen:locmn, data=dat2) coef(m4) # Match Hwu table 4#> genG1 genG2 genG3 genG4 genG5 genG1:locmn #> 24.4698520 -176.3280420 0.5186791 366.5962974 -215.2567865 1.0953375 #> genG2:locmn genG3:locmn genG4:locmn genG5:locmn #> 1.1140580 0.8111077 0.8067789 1.1727179# Table 6: Shukla's heterogeneity test dat2$ge = paste0(dat2$gen, dat2$env) # Create a separate ge interaction term m6 <- lm(yield ~ gen + env + ge + ge:locmn, data=dat2) m6b <- lm( yield ~ gen + env + ge + locmn, data=dat2) anova(m6, m6b) # Non-significant difference#> Analysis of Variance Table #> #> Model 1: yield ~ gen + env + ge + ge:locmn #> Model 2: yield ~ gen + env + ge + locmn #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 90 15810067 #> 2 90 15810067 0 0# Table 7 - Shukla stability # First, environment means emn <- group_by(dat2, env) emn <- summarize(emn, ymn=mean(yield))#>#> env #> gen E1 E2 E3 E4 E5 E6 #> G1 98.3000 -33.5500 291.1000 -72.2500 -320.9500 37.3500 #> G2 451.1333 -382.2167 -99.0667 58.5833 181.1333 -209.5667 #> G3 190.2167 529.6167 -297.7333 -343.8333 -197.7833 119.5167 #> G4 -349.3667 273.2833 -453.0667 124.0833 310.1333 94.9333 #> G5 -390.2833 -387.1333 558.7667 233.4167 27.4667 -42.2333for (ll in 1:nrow(getab)){ m7l <- lm(getab[ll, ] ~ emn$ymn) cat("\n\n*************** Gen ",ll," ***************\n") cat("Regression coefficient: ",round(coefficients(m7l)[2],5),"\n") print(anova(m7l)) } # Match Hwu table 7.#> #> #> *************** Gen 1 *************** #> Regression coefficient: 0.09534 #> Analysis of Variance Table #> #> Response: getab[ll, ] #> Df Sum Sq Mean Sq F value Pr(>F) #> emn$ymn 1 37638 37638 0.8987 0.3968 #> Residuals 4 167514 41879 #> #> #> *************** Gen 2 *************** #> Regression coefficient: 0.11406 #> Analysis of Variance Table #> #> Response: getab[ll, ] #> Df Sum Sq Mean Sq F value Pr(>F) #> emn$ymn 1 53870 53870 0.5586 0.4963 #> Residuals 4 385715 96429 #> #> #> *************** Gen 3 *************** #> Regression coefficient: -0.18889 #> Analysis of Variance Table #> #> Response: getab[ll, ] #> Df Sum Sq Mean Sq F value Pr(>F) #> emn$ymn 1 147748 147748 1.377 0.3057 #> Residuals 4 429197 107299 #> #> #> *************** Gen 4 *************** #> Regression coefficient: -0.19322 #> Analysis of Variance Table #> #> Response: getab[ll, ] #> Df Sum Sq Mean Sq F value Pr(>F) #> emn$ymn 1 154598 154598 1.6804 0.2646 #> Residuals 4 368004 92001 #> #> #> *************** Gen 5 *************** #> Regression coefficient: 0.17272 #> Analysis of Variance Table #> #> Response: getab[ll, ] #> Df Sum Sq Mean Sq F value Pr(>F) #> emn$ymn 1 123529 123529 0.9018 0.3961 #> Residuals 4 547906 136976# } # dontrun