#---------------------------------
# Prueba de Durbin para bloques incompletos
# k = nro de unidades experimentales por bloque
# r = numero de veces que aparece el tratamiento
# Lambda = numero de bloques en el que dos tratamientos estan juntos (dise‰o)
durbin.comp<-function(juez,trat,eval,k,r,alpha=0.05,main=" ") {
x<-data.frame(juez,trat,eval)
t<-length(unique(x$trat))
b<-length(unique(x$juez))
# Determina el rango dentro de cada juez
z <- by(x,x$juez,function(x) rank(x$eval))
y<-data.frame(c(z))
m<-dim(y)
n<-m[1]*m[2]
rango <- 1:n
for (i in 1:m[1]) {
for (j in 1:m[2]) {
kk=i+m[1]*(j-1)
rango[kk]<-y[i,j]
}
}
x<-data.frame(x,rango)
z <-by(x,x$trat,function(x) sum(x$rango))
y<-as.vector(c(z))
s <- (y-r*(k+1)/2)^2
s1 <- sum(s)
# determina el valor de Durbin
gl1<-t-1 ;gl2<-b*k-t-b+1
s <- 12*(t-1)*s1/(r*t*(k-1)*(k+1))
prob<-1-pchisq(s,gl1); Tprob<-qt(1-alpha/2,gl2)
Mc <-Tprob*sqrt( r*k*(k+1)*(b*(k-1)-s)/(6*gl2))

# s,prob,Tprob,Mc,gl1,gl2)
# Impresion de resultados
cat("\nEstudio:",main,"\n")
cat("Prueba de Durbin\n\n")
cat("Valor calculado    = ",s,"\n")
cat("Grados de libertad = ",gl1,"\n")
cat("P-valor            = ",prob,"\n\n")
cat("t-Student          = ",Tprob,"\n")
cat("Grados de libertad = ",gl2,"\n\n")
cat("Diferencia minima\npara la diferencia entre suma de rangos = ",Mc,"\n\n")

# comparacion de tratamientos.
comb <-combn(t,2)
n<-ncol(comb)
diff<-array(0,n)
stat<-array(" ns  ",n)
for (kk in 1:n) {
diff[kk]<-abs(y[comb[1,kk]]-y[comb[2,kk]])
if (diff[kk] >= Mc) stat[kk]<-" *   "
}
tcomb<-t(comb)
tstat<-cbind(stat)
tr.comp<-cbind(diff)
cat("Suma de rangos\n\Trat Suma\n----------\n")
for(i in 1:t){
a <-c(i)
b <-c(" = ",y[i])
cat(a,b,"\n")
}
compara <- data.frame(tcomb,tr.comp,tstat)
return(compara)
}

#---------------------------------
# Prueba de Durbin para bloques incompletos
# k = nro de unidades experimentales por bloque
# r = numero de veces que aparece el tratamiento
# Lambda = numero de bloques en el que dos tratamientos estan juntos (dise‰o)
durbin.grupos<-function(juez,trat,eval,k,r,alpha=0.05, main=" ") {
name.y <- paste(deparse(substitute(eval)))
x<-data.frame(juez,trat,eval)
t<-length(unique(x$trat))
b<-length(unique(x$juez))
# Determina el rango dentro de cada juez
z <- by(x,x$juez,function(x) rank(x$eval))
y<-data.frame(c(z))
m<-dim(y)
n<-m[1]*m[2]
rango <- 1:n
for (i in 1:m[1]) {
for (j in 1:m[2]) {
kk=i+m[1]*(j-1)
rango[kk]<-y[i,j]
}
}
x<-data.frame(x,rango)
z <-by(x,x$trat,function(x) sum(x$rango))
y<-as.vector(c(z))
nombre<-as.character(dimnames(z)$"x$trat")
s <- (y-r*(k+1)/2)^2
s1 <- sum(s)
# determina el valor de Durbin
gl1<-t-1 ;gl2<-b*k-t-b+1
s <- 12*(t-1)*s1/(r*t*(k-1)*(k+1))
prob<-1-pchisq(s,gl1); Tprob<-qt(1-alpha/2,gl2)
Mc <-Tprob*sqrt( r*k*(k+1)*(b*(k-1)-s)/(6*gl2))

# s,prob,Tprob,Mc,gl1,gl2)
# Impresion de resultados
cat("\nEstudio:",main)
cat("\nPrueba de Durbin")
cat("\n..............")
cat("\nChi Cuadrado : ",s)
cat("\nGl.          : ",gl1)
cat("\nP-valor      : ",prob)
cat("\n..............")
cat("\nComparaci♭n de tratamientos\n")
cat("\nAlpha        : ",alpha)
cat("\nGl.          : ",gl2)
cat("\nt-Student    : ",Tprob)
cat("\nDiferencia minima\npara la diferencia entre suma de rangos = ",Mc)
# comparacion de tratamientos.
cat("\n\nGrupos, Tratamientos y la Suma de sus rangos\n")
y<-as.numeric(y)
resultado<-orden.est(nombre,y,Mc)
#
return(resultado)
}

#-----------------
# establece el orden estadistico en la comparacion de tratamientos
# Nombre de tratamientos : trat
# promedios : prom
# valor minimo de comparacion

orden.est<-function(tratamiento,promedio,minimo) {
n<-length(promedio)
z<-data.frame(tratamiento,promedio)
# ordena tratamientos
w<-z[order(z[,2],decreasing=T),]
M<-rep("",n)
k<-1
k1<-0
j<-1
i<-1
cambio<-n
cambio1<-0
chequeo=0
M[1]<-letters[k]
while(j<n) {
chequeo<-chequeo+1
if (chequeo > n) break
for(i in j:n) {
s<-abs(w[i,2]-w[j,2])<=minimo
if(s) {
if(ultimo.c(M[i]) != letters[k])M[i]<-paste(M[i],letters[k],sep="")
}
else {
k<-k+1
cambio<-i
cambio1<-0
ja<-j
for(jj in cambio:n) M[jj]<-paste(M[jj]," ",sep="")
M[cambio]<-paste(M[cambio],letters[k],sep="")
for( v in ja:cambio) {
if(abs(w[v,2]-w[cambio,2])>minimo) {j<-j+1
cambio1<-1
}
else break
}
break
}
}
if (cambio1 ==0 )j<-j+1
}
#-----------
w<-data.frame(w,stat=M)
trat<-as.character(w$tratamiento)
prom<-as.numeric(w$promedio)
for(i in 1:n){
cat(M[i],"\t",trat[i],"\t",prom[i],"\n")
}
resultado<-data.frame(trat,prom,M)
return(resultado)
}
ultimo.c<-function(x) {
y<-sub(" +$", "",x)
p1<-nchar(y)
cc<-substr(y,p1,p1)
return(cc)
}
##
