#' Stability analysis. SHUKLA'S STABILITY VARIANCE AND KANG'S
#'
#' This procedure calculates the stability variations as well as the statistics
#' of selection for the yield and the stability. The averages of the genotype
#' through the different environment repetitions are required for the
#' calculations. The mean square error must be calculated from the joint
#' variance analysis.
#'
#' Stable (i) determines the contribution of each genotype to GE interaction by
#' calculating var(i); (ii) assigns ranks to genotypes from highest to lowest
#' yield receiving the rank of 1; (iii) calculates protected LSD for mean yield
#' comparisons; (iv) adjusts yield rank according to LSD (the adjusted rank
#' labeled Y); (v) determines significance of var(i) usign an aproximate
#' F-test; (vi) assigns stability rating (S) as follows: -8, -4 and -2 for
#' var(i) significant at the 0.01, 0.05 and 0.10 probability levels, and 0 for
#' nonsignificant var(i) ( the higher the var(i), the less stable the
#' genotype); (vii) sums adjusted yield rank, Y, and stability rating, S, for
#' each genotype to determine YS(i) statistic; and (viii) calculates mean YS(i)
#' and identifies genotypes (selection) with YS(i) > mean YS(i).
#'
#' @param data matrix of averages, by rows the genotypes and columns the
#' environment
#' @param rep Number of repetitions
#' @param MSerror Mean Square Error
#' @param alpha Label significant
#' @param main Title
#' @param cova Covariable
#' @param name.cov Name covariable
#' @param file.cov Data covariable
#' @param console logical, print output
#' @return \item{analysis}{ Analysis of variance } \item{statistics}{
#' Statistics of the model } \item{stability}{ summary stability analysis }
#' @author Felipe de Mendiburu
#' @seealso \code{\link{stability.nonpar} }
#' @references Kang, M. S. 1993. Simultaneous selection for yield and
#' stability: Consequences for growers. Agron. J. 85:754-757. Manjit S. Kang
#' and Robert Mangari. 1995. Stable: A basic program for calculating stability
#' and yield-stability statistics. Agron. J. 87:276-277
#' @keywords models
#' @importFrom stats qf
#' @export
#' @examples
#'
#' library(agricolae)
#' # example 1
#' # Experimental data,
#' # replication rep= 4
#' # Mean square error, MSerror = 1.8
#' # 12 environment
#' # 17 genotype = 1,2,3,.., 17
#' # yield averages of 13 genotypes in localities
#' f <- system.file("external/dataStb.csv", package="agricolae")
#' dataStb<-read.csv(f)
#' stability.par(dataStb, rep=4, MSerror=1.8, alpha=0.1, main="Genotype",console=TRUE)
#'
#' #example 2 covariable. precipitation
#' precipitation<- c(1000,1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,2100)
#' stability.par(dataStb, rep=4, MSerror=1.8, alpha=0.1, main="Genotype",
#' cova=TRUE, name.cov="Precipitation", file.cov=precipitation,console=TRUE)
#'
stability.par <-
function(data,rep,MSerror,alpha=0.1,main=NULL,cova=FALSE,
name.cov=NULL,file.cov=0,console=FALSE){
#-----------
KK <- "Environmental index"
y<-data
if (cova) {
x<- as.matrix(file.cov)
KK <- name.cov # cuando hay covariable
}
A <- nrow(y) # Number of genotypes
M <- ncol(y) # Number of environments
N <- rep # Number of replications
MKE <- MSerror # Pooled error mean square
RR <- main # "Stability" Title of the experiments
#
FM0 <- qf(1-alpha,M-1,M * (A - 1) * (N - 1))
# Data
#dimension for genotype (A)
dimA <- rep(0,A); dim(dimA)<-A
SHV=MV=SU=MV1=FF=GY=U1=G=SI=B=SA=S=FS=FSS=R=GYS=F1=GY=NN=X1=X1M=W=GYY=MMM=dimA
#dimension for environment (M)
dimA <- rep(0,M); dim(dimA)<-M
SHM= MM= II= X2= X2M=dimA
#dimension (A,M)
dimA <- rep(0,A*M); dim(dimA)<-c(A,M)
U= G1=dimA
SV = SM = GG1 = L= SMES = 0
SS<- sum(y^2)
SHT<- sum(y)
#..............
SHV<-apply(y,1,sum)
MV<-SHV / M
SV<-sum(SHV^2)/M
FK <- SHT ^ 2 / (A * M)
SKV <- (SV - FK) * N
MKV = SKV / (A - 1)
#..............
SHM<-apply(y,2,sum)
MM<-SHM / A
SM <-sum(SHM^2)/A
II<-MM - SHT / (A * M)
#..............
SKM <- (SM - FK) * N
MKM <- SKM / (M - 1)
SKT <- (SS - FK) * N
SKVM <-SKT - SKV - SKM; MKVM <- SKVM / ((A - 1) * (M - 1))
y<-as.matrix(y)
for ( i in 1: A) {
U[i, ] <- y[i, ] - MM
}
SU<-apply(U,1,sum)
U1 <- SU/ M
#...............
b <- 1 / ((M - 1) * (A - 1) * (A - 2))
for ( i in 1:A) {
G[i] <- sum((U[i,] - U1[i])^2)
}
GG1 <- sum(G)
SI <- b * (A * (A - 1) * G - GG1) * N
if ( cova ) ZZ <- sum(x^2)
if (!cova ) IN <- sum(II^2)
for ( i in 1: A) {
if ( cova ) B[i] <- sum( (U[i,] - U1[i])* x)/ZZ
if (!cova ) B[i] <- sum( (U[i,] - U1[i])*II)/IN
}
#...............
for ( i in 1: A) {
if ( cova ) SA[i] <- sum((U[i,]-U1[i]- B[i]*x)^2)
if (!cova ) SA[i] <- sum((U[i,]-U1[i]- B[i]*II)^2)
}
#...............
L <- sum(SA)
S <- (A / ((A - 2) * (M - 2))) * (SA - L / (A * (A - 1))) * N
KI <- ((A - 1) * (M - 2)) / A
SS<-sum(S * KI)
#...............
SKMB <- SS
MS <- SKMB / ((A - 1) * (M - 2))
SKH <- SKVM - SKMB
MKH <- SKH / (A - 1)
SHKLT <- M * (A - 1) * (N - 1)
T05 <- qt(0.95,SHKLT)
DMV05 <- T05 * sqrt(2 * MKE / (N * M))
SMES <- sum(MV)
MES <- SMES / A
for ( i in 1: A) {
if( MV[i] > MES) MV1[i] <- 1
if( MV[i] >= (MES + DMV05) ) MV1[i] <- 2
if( MV[i] >= (MES + 2 * DMV05)) MV1[i] <- 3
if( MV[i] < MES ) MV1[i] <- -1
if( MV[i] <= (MES - DMV05)) MV1[i] <- -2
if( MV[i] <= (MES - 2 * DMV05)) MV1[i] <- -3
}
#-------
FV <- MKV / MKVM; FM <- MKM / MKE
FVM <-MKVM / MKE; FH <- MKH / MS
FMS <- MS / MKE
FS <- SI / MKE
FSS <- S / MKE
SHF2 <-(A - 1) * (M - 1); SHF1 = (A - 1)
# value F and t-student
F05 <- qf(0.95,SHF1, SHF2)
F01 <- qf(0.99,SHF1, SHF2)
pvalue <- round(1-pf( FV ,SHF1, SHF2),3)
DD<-paste("",pvalue)
if (pvalue <0.001) DD<-"<0.001"
SHF2 <- M * (A - 1) * (N - 1)
SHF1 <- M - 1
F05 <- qf(0.95,SHF1, SHF2)
F01 <- qf(0.99,SHF1, SHF2)
pvalue<- round(1-pf( FM ,SHF1, SHF2),3)
NNN<-paste("",pvalue)
if (pvalue <0.001) NNN<-"<0.001"
SHF2 <- M * (A - 1) * (N - 1); SHF1 <- (A - 1) * (M - 1)
F05 <- qf(0.95,SHF1, SHF2)
F01 <- qf(0.99,SHF1, SHF2)
pvalue <- round(1-pf( FVM ,SHF1, SHF2),3)
LL<-paste("",pvalue)
if (pvalue <0.001) LL<-"<0.001"
SHF2 <- (A - 1) * (M - 2); SHF1 <- A - 1
F05 <- qf(0.95,SHF1, SHF2)
F01 <- qf(0.99,SHF1, SHF2)
pvalue <- round(1-pf( FH ,SHF1, SHF2),3)
HH<-paste("",pvalue)
if (pvalue <0.001) HH<-"<0.001"
SHF2 <- M * (A - 1) * (N - 1); SHF1 <- (A - 1) * (M - 2)
F05 <- qf(0.95,SHF1, SHF2)
F01 <- qf(0.99,SHF1, SHF2)
pvalue <- round(1-pf( FMS ,SHF1, SHF2),4)
BB<-paste("",pvalue)
if (pvalue <0.001) BB<-"<0.001"
SHF2 <- M * (A - 1) * (N - 1)
SHF1 <- M - 1
F05 <- qf(0.95,SHF1, SHF2)
F01 <- qf(0.99,SHF1, SHF2)
for ( i in 1: A) {
if(FS[i] >= F01 ) NN[i] <- "**"
if( (FS[i] < F01) & (FS[i] >= F05) ) NN[i] <- "*"
if( FS[i] < F05 ) NN[i] <- "ns"
}
SHF2 <- M * (A - 1) * (N - 1)
SHF1 <- M - 2
F05 <- qf(0.95,SHF1, SHF2)
F01 <- qf(0.99,SHF1, SHF2)
for ( i in 1: A) {
if( FSS[i] >= F01 ) MMM[i] <- "**"
if( (FSS[i] < F01) & (FSS[i] >= F05) ) MMM[i] <- "*"
if( FSS[i] < F05 ) MMM[i] <- "ns"
}
#-------
if(console){
cat("\n","INTERACTIVE PROGRAM FOR CALCULATING SHUKLA'S STABILITY VARIANCE AND KANG'S")
cat("\n"," YIELD - STABILITY (YSi) STATISTICS")
cat("\n",RR,"\n",KK," - covariate \n")
cat("\n","Analysis of Variance\n")
#cat("\n")
#cat("\n","Source d.f. Sum of Squares Mean Squares F p.value")
#cat("\n",rep("-",35))
}
fuentes<- c( "Total ","Genotypes", "Environments","Interaction",
"Heterogeneity" , "Residual","Pooled Error")
gl <- c(A*M-1, A-1, M-1, (M - 1) * (A - 1), A - 1,(A - 1) * (M - 2),M * (A - 1) * (N - 1))
SC <- round(c(SKT, SKV, SKM , SKVM ,SKH, SKMB),4)
SC <- c(SC,0)
CM <- round(c(MKV, MKM, MKVM ,MKH ,MS,MKE),4)
CM <- c(0,CM)
Fcal<-round(c(FV,FM,FVM,FH,FMS),2)
Fcal<-c(0,Fcal,0)
resul<-c(" ",DD,NNN,LL,HH,BB," ")
Z<-data.frame(gl,SC,CM,Fcal, resul)
names(Z)<-c("Df","Sum Sq","Mean Sq","F value","Pr(>F)")
rownames(Z)<- fuentes
Z[7,c(2,4)]<-" ";Z[1,c(3,4)]<-" ";
if(console){
cat("\n")
print(Z)
}
Z1<-Z # Analysis
#Z<-as.matrix(Z)
#Z[1,4:6] <-" "
#Z[7,c(3,5,6)]<-" "
#for ( i in 1: 7) {
#cat("\n", Z[i,1],"\t",Z[i,2],"\t", Z[i,3],"\t", Z[i,4], "\t",Z[i,5], Z[i,6])
#}
#cat("\n",rep("-",35))
X1<-apply(y,1,sum)
X1M<- X1/ M
X2<-apply(y,2,sum)
X2M <- X2 / A
SH <- sum(X2)
MM1 <- SH / (A * M)
for ( i in 1: A) {
W[i] <- sum((y[i,] - X1M[i] - X2M + MM1) ^ 2) * N
}
#................
MV<-round(MV,6); SI<-round(SI,6); S<-round(S,6); W<-round(W,6);
Z<-data.frame(MV, SI, NN, S,MMM,W)
names(Z)<-c("Mean","Sigma-square",".","s-square",".","Ecovalence")
rownames(Z)<-rownames(y)
if(console){
cat("\n","Genotype. Stability statistics\n\n")
print(Z)
}
Z2<-Z # statistics
Z<-as.matrix(Z)
#for ( i in 1: A) {
#cat("\n",i, "\t", Z[i,1],"\t", Z[i,2], Z[i,3], "\t",Z[i,4], Z[i,5],"\t", Z[i,6])
#}
FF <- SI / MKE # each genotype
SHF2 <- M * (A - 1) * (N - 1); SHF1 <- M - 1
F05 <- qf(0.95,SHF1, SHF2)
F01 <- qf(0.99,SHF1, SHF2)
for ( i in 1: A) {
if( FF[i] < FM0) F1[i] <- 0
if( FF[i] >= FM0) F1[i] <- -2
if( FF[i] >= F05) F1[i] <- -4
if( FF[i] >= F01) F1[i] <- -8
}
for ( i in 1: A) {
R0 <- 1
for ( j in 1: A) {
if( MV[j] < MV[i]) R0 <- R0 + 1
}
R[i] <- R0
}
for ( i in 1: A) GY[i] <- R[i] + MV1[i]
for ( i in 1: A) GYS[i] <- GY[i] + F1[i]
SGYS<-0
for ( i in 1: A) SGYS <- SGYS + GYS[i]
MGYS <- SGYS / A
for ( i in 1: A) {
GYY[i]<-""
if( GYS[i] > MGYS ) GYY[i] <- "+"
}
names<-c("Yield","Rank","Adj.rank","Adjusted","Stab.var","Stab.rating","YSi","..." )
Z<-data.frame(MV, R, MV1,GY,SI,F1,GYS,GYY)
rownames(Z)<-rownames(y)
names(Z)<-names
Z3<-Z # stability
if(console){
cat("\n\nSignif. codes: 0 '**' 0.01 '*' 0.05 'ns' 1\n\n")
cat("Simultaneous selection for yield and stability (++)\n\n")
print(Z)
cat("\n","Yield Mean:", MES)
cat("\n","YS Mean:", MGYS)
cat("\n","LSD (0.05):", DMV05)
cat("\n",rep("-",11))
cat("\n","+ selected genotype" )
cat("\n","++ Reference: Kang, M. S. 1993. Simultaneous selection for yield")
cat("\n","and stability: Consequences for growers. Agron. J. 85:754-757." )
cat("\n")
}
out<-list(analysis=Z1,statistics =Z2,stability =Z3)
invisible(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.