Skip to content

Commit

Permalink
version 1.0-5
Browse files Browse the repository at this point in the history
  • Loading branch information
Felipe de Mendiburu authored and gaborcsardi committed Jul 25, 2008
1 parent b13ada1 commit 47c1125
Show file tree
Hide file tree
Showing 176 changed files with 1,614 additions and 1,242 deletions.
1 change: 0 additions & 1 deletion .cvsignore

This file was deleted.

18 changes: 0 additions & 18 deletions .project

This file was deleted.

8 changes: 4 additions & 4 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
Package: agricolae
Type: Package
Title: Statistical Procedures for Agricultural Research
Version: 1.0-4
Date: 2007-09-11
Version: 1.0-5
Date: 2008-07-25
Author: Felipe de Mendiburu
Maintainer: Felipe de Mendiburu <f.mendiburu@cgiar.org>
Suggests: akima, klaR, SuppDists, corpcor
Description: These functions are currently utilized by the International Potato Center Research (CIP), the Statistics and Informatics Instructors and the Students of the Universidad Nacional Agraria La Molina Peru, and the Specialized Master in "Bosques y Gestion de Recursos Forestales" (Forest Resource Management). This package contains functionality for the statistical analysis of experimental designs applied specially for field experiments in agriculture and plant breeding. Planning of field experiments: Lattice, factorial, RCBD, CRD, Latin Square, Greaco, BIB, PBIB, Alpha design. Comparison of multi-location trials: AMMI (biplot and triplot), Stability. Comparison between treatments: LSD, Bonferroni, HSD, Waller, Kruskal, Friedman, Durbin, Van Der Waerden. Resampling and simulation: resampling.model, simulation.model, analysis Mother and baby trials, Ecology: Indices Biodiversity, path analysis, consensus cluster, Uniformity Soil: Index Smith's.
Description: Agricolae is a project in order to obtain the degree of master in Systems Engineering in the National University of Engineering in Lima-Peru (UNI in Spanish). These functions are currently utilized by the International Potato Center Research (CIP), the Universidad Nacional Agraria La Molina (UNALM-PERU), and the Instituto Nacional de Investigacion Agricola (INIA-PERU). It comprises the functionality of statistical analysis into experimental designs applied specially for field experiments in agriculture and plant breeding: Lattice, factorial, complete and incomplete block, Latin Square, Greaco, Alpha designs, Cyclic designs, comparison of multi-location trials, comparison between treatments, resampling, simulation, biodiversity indexes and consensus cluster.
License: GPL
URL: http://tarwi.lamolina.edu.pe/~fmendiburu
Packaged: Tue Sep 11 20:55:49 2007; hornik
Packaged: Fri Jul 25 19:04:19 2008; fdemendiburu
45 changes: 0 additions & 45 deletions Para instalar agricola en DOS.txt

This file was deleted.

17 changes: 8 additions & 9 deletions R/AMMI.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
AMMI <-
function (ENV, GEN, REP, Y, MSE = 0, number = TRUE, graph = "biplot",
`AMMI` <-
function (ENV, GEN, REP, Y, MSE = 0, number = TRUE, graph = "biplot",
...)
{
name.y <- paste(deparse(substitute(Y)))
Expand Down Expand Up @@ -134,15 +134,15 @@
nssammi<-nrow(SSAMMI)
SSAMMI<-SSAMMI[SSAMMI$Df>0,]
nss<-nrow(SSAMMI)
row.names(SSAMMI) <- paste("CP", 1:nss, sep = "")
row.names(SSAMMI) <- paste("PC", 1:nss, sep = "")
cat("\nAnalysis\n")

print(SSAMMI)
LL <- sqrt(diag(L))
SCOREG <- U %*% LL
SCOREE <- V %*% LL
SCORES <- rbind(SCOREG, SCOREE)
colnames(SCORES) <- paste("CP", 1:nssammi, sep = "")
colnames(SCORES) <- paste("PC", 1:nssammi, sep = "")
MSCORES <- SCORES[1:ngen, ]
NSCORES <- SCORES[(ngen + 1):(ngen + nenv), ]
MGEN <- data.frame(type = "GEN", Y = apply(OUTMED, 1, mean),
Expand All @@ -164,13 +164,13 @@
if (number == TRUE) {
plot(MGEN[, 3], MGEN[, 4], cex = 0, text(MGEN[, 3],
MGEN[, 4], labels = as.character(1:nrow(MGEN)),
col = "blue"), xlab = "CP 1", ylab = "CP 2",
col = "blue"), xlab = "PC 1", ylab = "PC 2",
frame = TRUE, ...)
}
if (number == FALSE) {
plot(MGEN[, 3], MGEN[, 4], cex = 0, text(MGEN[, 3],
MGEN[, 4], labels = row.names(MGEN), col = "blue"),
xlab = "CP 1", ylab = "CP 2", frame = TRUE, ...)
xlab = "PC 1", ylab = "PC 2", frame = TRUE, ...)
}
points(MENV[, 3], MENV[, 4], cex = 0, text(MENV[, 3],
MENV[, 4], labels = row.names(MENV), col = "brown"))
Expand All @@ -179,7 +179,7 @@
arrows(0, 0, 0.9 * MENV[, 3][s], 0.9 * MENV[, 4][s],
col = "brown", lwd = 1.8, length = 0.1, code = 2)
legend("topleft", NULL, pch = c("1", "2"), cp.per[1:2],
, title = "CP %", lty = 0)
, title = "PC %", lty = 0)
}
if (graph == "triplot") {
y <- bplot
Expand Down Expand Up @@ -209,7 +209,7 @@
suppressWarnings(warning(text(tritrafo(point3), cp.name,
adj = c(0.5, 0), cex = 1)))
legend("topleft", NULL, pch = c("1", "2", "3"), cp.per,
, title = "CP %", lty = 0)
, title = "PC %", lty = 0)
trilines(centerlines(3), lty = 2.5, col = "green", lwd = 2)
for (i in 1:nlugar) {
suppressWarnings(warning(trilines(c(point2[i, 1],
Expand All @@ -220,4 +220,3 @@
return(list(genXenv=OUTRES2, analysis=SSAMMI, means=MEDIAS, biplot=bplot))
}


5 changes: 3 additions & 2 deletions R/AMMI.contour.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@
function(model,distance,shape,...)
{
G<- subset(model$biplot,model$biplot$type=="GEN")
x<-G$CP1
y<-G$CP2
x<-G$PC1
y<-G$PC2
d<-sqrt(x^2+y^2)
r <-distance*max(d)
x<-seq(-r,r,length=shape)
Expand All @@ -21,3 +21,4 @@ cat("\nGenotype out:",length(GEN.out),"\n\n")
distance<-data.frame(row.names=row.names(G),distance=d)
return(list("GENOTYPE IN"=GEN.in, "GENOTYPE OUT"=GEN.out,Distance=distance))
}

206 changes: 110 additions & 96 deletions R/BIB.test.R
Original file line number Diff line number Diff line change
@@ -1,100 +1,114 @@
`BIB.test` <-
function(block,trt,y, method="lsd", alpha=0.05,group=TRUE)
function (block, trt, y, method = "lsd", alpha = 0.05, group = TRUE)
{
block.unadj<-as.factor(block)
trt.adj <-as.factor(trt)
name.y <- paste(deparse(substitute(y)))
model<-lm(y ~ block.unadj+trt.adj)
DFerror<-df.residual(model)
MSerror<-deviance(model)/DFerror
k <-unique(table(block.unadj))
r <-unique(table(trt.adj))
b <-nlevels(block.unadj)
ntr <-nlevels(trt.adj)
lambda<-r*(k-1)/(ntr-1)
tabla<-suppressWarnings(mxyz(block,trt,y ))
AA<-!is.na(tabla)
BB<-tapply(y,block.unadj,sum)
B<-BB%*%AA
Y<-tapply(y,trt.adj,sum)
Q<-Y-as.numeric(B)/k

SStrt.adj<-sum(Q^2)*k/(lambda*ntr)
MStrt.adj<- SStrt.adj/(ntr-1)
sdtdif<-sqrt(2*k*MSerror/(lambda*ntr))
Fvalue<- MStrt.adj/MSerror
# mean adjusted.

mean.adj<-mean(y)+Q*k/(lambda*ntr)
sdmean.adj <- sqrt(MSerror*(1+k*r*(ntr-1)/(lambda*ntr))/(r*ntr))
cat("\nANALYSIS BIB: ", name.y, "\nClass level information\n")
cat("\nBlock: ", unique(as.character(block)))
cat("\nTrt : ", unique(as.character(trt)))
cat("\n\nNumber of observations: ", length(y), "\n\n")
print(anova(model))
cat("coefficient of variation:",round(cv.model(model),1),"%\n")
cat("\nTreatments\n")
print(data.frame( row.names=NULL,trt=row.names(Y),means=Y/r,mean.adj,sdmean.adj))
parameter<- k/(lambda*ntr)
if (method=="lsd") {
Tprob<-qt(1-alpha/2,DFerror)
cat("\nLSD test\n")
cat("\nAlpha :",alpha)
cat("\nLSD :",Tprob*sdtdif)
}
if (method=="tukey") {
Tprob <- qtukey(1-alpha, ntr, DFerror)
cat("\nTukey\n")
cat("\nAlpha :",alpha)
cat("\nHSD :",Tprob*sdtdif)
parameter<-parameter/2
}
if (method=="waller") {
K<-650-16000*alpha+100000*alpha^2
Tprob<-waller(K,ntr-1,DFerror,Fvalue)
cat("\nWaller-Duncan K-ratio\n")
cat("\nThis test minimizes the Bayes risk under additive")
cat("\nloss and certain other assumptions.\n")
cat("\nk Ratio: ",K)
cat("\nMSD :",Tprob*sdtdif)
}
E<-lambda*ntr/(r*k)
cat("\nParameters BIB")
cat("\nLambda :",lambda)
cat("\ntreatmeans :",ntr)
cat("\nBlock size :",k)
cat("\nBlocks :",b)
cat("\nReplication:",r,"\n")
cat("\nEfficiency factor",E,"\n\n<<< Book >>>\n")
if (group) {
cat("\nMeans with the same letter are not significantly different.")
cat("\n\nComparison of treatments\n\nGroups, Treatments and means\n")
output <- order.group(names(mean.adj), as.numeric(mean.adj), rep(1,ntr),
MSerror, Tprob,std.err=sdmean.adj,parameter)
output[,4]<-r
block.unadj <- as.factor(block)
trt.adj <- as.factor(trt)
name.y <- paste(deparse(substitute(y)))
modelo <- formula(paste(name.y,"~ block.unadj + trt.adj"))
model <- lm(modelo)
DFerror <- df.residual(model)
MSerror <- deviance(model)/DFerror
k <- unique(table(block.unadj))
r <- unique(table(trt.adj))
b <- nlevels(block.unadj)
ntr <- nlevels(trt.adj)
lambda <- r * (k - 1)/(ntr - 1)
datos <- data.frame(block, trt, y)
tabla <- by(datos[,3], datos[,1:2], function(x) mean(x,na.rm=TRUE))
tabla <-as.data.frame(tabla[,])
AA <- !is.na(tabla)
BB <- tapply(y, block.unadj, sum)
B <- BB %*% AA
Y <- tapply(y, trt.adj, sum)
Q <- Y - as.numeric(B)/k
SStrt.adj <- sum(Q^2) * k/(lambda * ntr)
MStrt.adj <- SStrt.adj/(ntr - 1)
sdtdif <- sqrt(2 * k * MSerror/(lambda * ntr))
Fvalue <- MStrt.adj/MSerror
mean.adj <- mean(y) + Q * k/(lambda * ntr)
StdError.adj <- sqrt(MSerror * (1 + k * r * (ntr - 1)/(lambda *
ntr))/(r * ntr))
cat("\nANALYSIS BIB: ", name.y, "\nClass level information\n")
cat("\nBlock: ", unique(as.character(block)))
cat("\nTrt : ", unique(as.character(trt)))
cat("\n\nNumber of observations: ", length(y), "\n\n")
print(anova(model))
cat("\ncoefficient of variation:", round(cv.model(model), 1),
"%\n")
cat(name.y, "Means:", mean(y,na.rm=TRUE), "\n")
cat("\nTreatments\n")
print(data.frame(row.names = 1:ntr, trt = row.names(Y), means = Y/r,
mean.adj, StdError.adj))
parameter <- k/(lambda * ntr)
if (method == "lsd") {
Tprob <- qt(1 - alpha/2, DFerror)
cat("\nLSD test")
cat("\nStd.diff :", sdtdif)
cat("\nAlpha :", alpha)
cat("\nLSD :", Tprob * sdtdif)
}
if (method == "tukey") {
Tprob <- qtukey(1 - alpha, ntr, DFerror)
cat("\nTukey")
cat("\nAlpha :", alpha)
cat("\nStd.err :", sdtdif)
cat("\nHSD :", Tprob * sdtdif)
parameter <- parameter/2
}
if (method == "waller") {
K <- 650 - 16000 * alpha + 1e+05 * alpha^2
Tprob <- waller(K, ntr - 1, DFerror, Fvalue)
cat("\nWaller-Duncan K-ratio")
cat("\nThis test minimizes the Bayes risk under additive")
cat("\nloss and certain other assumptions.\n")
cat("\nk Ratio : ", K)
cat("\nMSD :", Tprob * sdtdif)
}
E <- lambda * ntr/(r * k)
cat("\n\nParameters BIB")
cat("\nLambda :", lambda)
cat("\ntreatmeans :", ntr)
cat("\nBlock size :", k)
cat("\nBlocks :", b)
cat("\nReplication:", r, "\n")
cat("\nEfficiency factor", E, "\n\n<<< Book >>>\n")
if (group) {
cat("\nMeans with the same letter are not significantly different.")
cat("\n\nComparison of treatments\n\nGroups, Treatments and means\n")
output <- order.group(names(mean.adj), as.numeric(mean.adj),
rep(1, ntr), MSerror, Tprob, std.err = StdError.adj,
parameter)
output[, 4] <- r
}
if (!group) {
comb <- combn(ntr, 2)
nn <- ncol(comb)
dif <- rep(0, nn)
pvalue <- rep(0, nn)
for (k in 1:nn) {
i <- comb[1, k]
j <- comb[2, k]
dif[k] <- abs(mean.adj[i] - mean.adj[j])
if (method == "lsd")
pvalue[k] <- 2 * round(1 - pt(dif[k]/sdtdif,
DFerror), 4)
if (method == "tukey")
pvalue[k] <- round(1 - ptukey(dif[k] * sqrt(2)/sdtdif,
ntr, DFerror), 4)
}
if (method == "waller")
significant = dif > Tprob * sdtdif
tr.i <- comb[1, ]
tr.j <- comb[2, ]
cat("\nComparison between treatments means\n")
if (method == "waller")
print(data.frame(row.names = NULL, tr.i, tr.j, diff = dif,
significant))
else print(data.frame(row.names = NULL, tr.i, tr.j, diff = dif,
pvalue))
output <- data.frame(trt = names(mean.adj), means = as.numeric(mean.adj),
M = "", N = r, std.err = StdError.adj)
}
return(output)
}
if (!group) {
comb <-combn(ntr,2)
nn<-ncol(comb)
dif<-rep(0,nn)
pvalue<-rep(0,nn)
for (k in 1:nn) {
i<-comb[1,k]
j<-comb[2,k]
dif[k]<-abs(mean.adj[i]-mean.adj[j])
if (method=="lsd") pvalue[k]<- 2*round(1-pt(dif[k]/sdtdif,DFerror),4)
if (method=="tukey") pvalue[k]<- round(1-ptukey(dif[k]*sqrt(2)/sdtdif,ntr,DFerror),4)
}
if (method=="waller") significant = dif > Tprob*sdtdif

tr.i<-comb[1,]
tr.j<-comb[2,]

cat("\nComparison between treatments means\n")
if (method=="waller") print(data.frame(row.names=NULL,tr.i,tr.j,diff=dif,significant))
else print(data.frame(row.names=NULL,tr.i,tr.j,diff=dif,pvalue))
output<-data.frame(trt= names(mean.adj),means= as.numeric(mean.adj),M="",
N=r,std.err=sdmean.adj)
}
return(output)
}
9 changes: 5 additions & 4 deletions R/HSD.test.R
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
"HSD.test" <-
`HSD.test` <-
function (y, trt, DFerror, MSerror, alpha=0.05, group=TRUE,main = NULL)
{
name.y <- paste(deparse(substitute(y)))
name.t <- paste(deparse(substitute(trt)))
junto <- subset(data.frame(y, trt), is.na(y) == FALSE)
means <- tapply.stat(junto[,2],junto[,1],stat="mean")
sds <- tapply.stat(junto[,2],junto[,1],stat="sd")
nn <- tapply.stat(junto[,2],junto[,1],stat="length")
means <- tapply.stat(junto[,1],junto[,2],stat="mean") # change
sds <- tapply.stat(junto[,1],junto[,2],stat="sd") #change
nn <- tapply.stat(junto[,1],junto[,2],stat="length") # change
means<-data.frame(means,std.err=sds[,2]/sqrt(nn[,2]),replication=nn[,2])
names(means)[1:2]<-c(name.t,name.y)
# row.names(means)<-means[,1]
Expand Down Expand Up @@ -62,3 +62,4 @@ output<-data.frame(trt= means[,1],means= means[,2],M="",N=means[,4],std.err=mean
}
return(output)
}

Loading

0 comments on commit 47c1125

Please sign in to comment.