R/stability.nonpar.R

In agricolae: Statistical Procedures for Agricultural Research

stability.nonpar <-
function(data, variable= NULL, ranking = FALSE,console=FALSE)
{
row.names(data)<-data[,1]
data<-data[,-1]
audpc<-as.matrix(data)
rnames<-rownames(data)
cnames<-colnames(data)
n<-ncol(audpc); k<-nrow(audpc)
lineas<-matrix(0,k)
ldev <-matrix(0,k)
s2 <- matrix(0,k)
s1 <- matrix(0,k)
mr <- matrix(0,k)
ambientes<-matrix(0,1,n)
Maudpc <- mean(audpc)
for(i in 1:k) {
lineas[i]<-mean(audpc[i,])
ldev[i] <-lineas[i] - Maudpc }
cyld<-matrix(0,k,n)
rcyld<-matrix(0,k,n)
raudpc<-matrix(0,k,n)
dimnames(raudpc)<-list(rnames,cnames)
for(j in 1:n) ambientes[j]<-mean(audpc[,j])

for(i in 1:k){
for(j in 1:n) cyld[i,j] <- audpc[i,j] - ldev[i] }
for(j in 1:n) {
rcyld[,j]<-rank(cyld[,j])
raudpc[,j]<-rank(audpc[,j])}
#=========================================
nn=n*(n-1)/2
for (i in 1:k) {
sumq<-0
for(j in 1:(n-1)) for(jj in (j+1):n) sumq <- sumq+abs(rcyld[i,j]-rcyld[i,jj])
s1[i]=sumq/nn
}
#=========================================
for(i in 1:k) {
s2[i]<-var(rcyld[i,])
# mr[i]<-mean(rcyld[i,]) 
}

#Calculation of expectation and variance

k2<-k**2;                     
es1<-(k2-1)/(3*k) 

es2<-(k2-1)/12

vs1<-k2*((k**2-4)*(n+3)+30)/(45*k**2*n*(n-1))
vs2<-(k2-1)*(2*(k2-4)*(n-1)+5*(k2-1))/(360*n*(n-1))
alphas<-1-0.05/k
chi.ind<-qchisq(alphas,1)
chi.sum<-qchisq(0.95,k)
#=========================================
z1<-(s1-es1)**2/vs1
z2<-(s2-es2)**2/vs2
suma.z1 <- sum(z1)
suma.z2 <- sum(z2)
mr<-rank(lineas)
stat1 <- data.frame(lineas,mr,s1,z1,s2,z2)
row.names(stat1)<-rnames
names(stat1)<- c("Mean","Rank", "s1", "Z1", "s2", "Z2")
stat2 <- data.frame(MEAN=Maudpc,es1,es2,vs1,vs2,chi.ind,chi.sum)
if(console){
cat("\n")
cat("Nonparametric Method for Stability Analysis\n")
cat("-------------------------------------------\n\n")
cat("Estimation and test of nonparametric measures\n")
cat("Variable:",variable,"\n\n")
if (ranking) {
cat("Ranking...\n")
print(raudpc)
cat("\n")
}
cat("Statistics...\n")
print(round(stat1,2))
cat("------------------------")
cat("\nSum of Z1: ",suma.z1)
cat("\nSum of Z2: ",suma.z2)
cat("\n------------------------\n\n")
cat("Test...\n")
cat("The Z-statistics are measures of stability. The test for the significance\n")
cat("of the sum of Z1 or Z2 are compared to a Chi-Square value of chi.sum. \n")
cat("individual Z1 or Z2 are compared to a Chi-square value of chi.ind.\n\n")
print(stat2)
cat("---\n")
cat("expectation and variance: es1, es2, vs1, vs2\n")
}
out<-list(ranking=stat1,statistics =stat2)
invisible(out)
}