## -----------------------------------------------------------------------------
library(agro)
ff <- get_data_from_uri("https://dataverse.harvard.edu/dataset.xhtml?persistentId=hdl:1902.1/18485", ".")
ff


## ----import-------------------------------------------------------------------
library(readxl)
d <- read_xlsx(ff, skip = 22)


## ----organize-----------------------------------------------------------------
colnames(d)[c(12, 30)]
colnames(d)[12] <- "height"
colnames(d)[30] <- "yield"
colnames(d)[c(12, 30)]


## -----------------------------------------------------------------------------
d$yield <- d$yield / 1000


## -----------------------------------------------------------------------------
table(d$N)
d$N <- c(0,30,60,120)[d$N + 1]
table(d$N)

table(d$V)
d$V <- c("Nerica 1", "Nerica 2", "Nerica 3", "WAB 56-104")[d$V]
table(d$V)

table(d$D)
d$D <- c("dib30", "dib20", "drill")[d$D]
table(d$D)


## ----table1a------------------------------------------------------------------
d$D <- as.factor(d$D)
d$V <- as.factor(d$V)


## ----table1b------------------------------------------------------------------
m <- lm(height ~ N * D * V + I(N^2), data = d)
anova(m)


## ----figure2a-----------------------------------------------------------------
boxplot(height ~ N, data = d)


## ----out1---------------------------------------------------------------------
nols.2q <- function(x, var="height") {
    m <- median(x[[var]], na.rm=TRUE)
    r <- quantile(x[[var]], c(0.25, 0.75), na.rm=TRUE)
    r <- 2 * diff(r)
    i <- which(abs(x[[var]] - m) < r)
    x[i,]
}

Nlev <- unique(d$N)
out <- vector(length=length(Nlev), mode="list")
for (i in 1:length(Nlev)) {
  dd <- d[d$N == Nlev[i], ]
  out[[i]] <- nols.2q(dd, var="height")
}
dd <- do.call(rbind, out)

boxplot(height ~ N, data = dd)


## ----nooutliers---------------------------------------------------------------
height_N <- aggregate(dd[, 'height', drop=FALSE], dd[, 'N', drop=FALSE], mean)

# Caculate standard errors (can also use aggregate function here)
height_N_sem <- tapply(d$height, d$N, sd) / sqrt(tapply(d$height, d$N, length))
  
# Make the bar plot
xp <- barplot(height~N, data=height_N, xlab="Nitrogen (kg/ha)", ylab="Plant height (cm)", ylim = c(0,100), space = 1) 

box(bty="l")

# Add error bars
arrows (p, (height_N$height - height_N_sem), p, (height_N$height + height_N_sem), code = 3, angle = 90, length = 0.05)



## ----table2a------------------------------------------------------------------
# Caculate mean grain yield for each N level
Yield_N <- aggregate(d[, "yield", drop=FALSE], d[, "N", drop=FALSE], mean)

# Perform the simple linear regression between yield and nitrogen level using the lm function
Nyield_analysis <- lm (yield ~ as.factor(N), data = d)


## ----table2b------------------------------------------------------------------
TukeyHSD(aov(Nyield_analysis))


## ----figure3------------------------------------------------------------------
yNV <- aggregate(d[, "yield", drop=FALSE], d[, c("N","V"), drop=FALSE], mean)
cols <- c("red", "blue", "orange", "black")

plot(yNV$N, yNV$yield, xlab = "N (kg/ha)", ylab = "Grain yield (ton/ha)", pch = c(15, 16, 17, 18)[yNV$V], col=cols[yNV$V], ylim=c(200, 1000), cex=1.5)

legend ("bottomleft", legend = levels(yNV$V), pch = c(15:18), col=cols, horiz=TRUE, cex=1.25)