## ----download-----------------------------------------------------------------
library(agro)
ff <- get_data_from_uri("https://doi.org/10.21223/P3/MIVGMU", ".")
ff


## ----importdata1--------------------------------------------------------------
library(openxlsx)
getSheetNames(ff[1])


## ----importdata2--------------------------------------------------------------
# tuber yield
ty <- read.xlsx(ff[1], sheet="TY")

# water use efficiency
wue <- read.xlsx(ff[1], sheet = "WUE", startRow = 6)

# crop water demand
cwd <- read.xlsx(ff[1], sheet = "CWD", startRow = 2)



## -----------------------------------------------------------------------------
rwc6 <- read.xlsx(ff[1], sheet = "RWC", rows=2:20)
rwc8 <- read.xlsx(ff[1], sheet = "RWC", startRow=22)


## ----table1-------------------------------------------------------------------
ty_WRIT6W <- ty[ty$T == "WRIT6w", ]
colnames(ty_WRIT6W)[4] <- "TY"

ty_WRIT8W <- ty[ty$T == "WRIT8w", ]
colnames(ty_WRIT8W)[4] <- "TY"


## ----tableout1----------------------------------------------------------------
par(mfrow=c(2,2))
boxplot(TY ~ TRT, ty_WRIT6W, outline = F, ylim=c(0,30))
boxplot(TY ~ B, ty_WRIT6W, outline = F, ylim=c(0,30))
boxplot(TY ~ TRT, ty_WRIT8W, outline = F, ylim=c(0,30))
boxplot(TY ~ B, ty_WRIT8W, outline = F, ylim=c(0,30))


## ----tableout5----------------------------------------------------------------
it6yield_analysis <- lm (ty_WRIT6W$TY ~ ty_WRIT6W$TRT + ty_WRIT6W$B)
anova(it6yield_analysis) 

it8yield_analysis <- lm (ty_WRIT8W$TY ~ ty_WRIT8W$TRT + ty_WRIT8W$B)
anova(it8yield_analysis) 


## ----tableout6----------------------------------------------------------------
colnames(wue)[4]
colnames(wue)[4] <- "WUE"

wue_WRIT6W <- wue[wue$T == "WRIT6w", ]
wue_WRIT8W <- wue[wue$T == "WRIT8w", ]


## ----tableout7----------------------------------------------------------------
it6wue_analysis <- lm (wue_WRIT6W$WUE ~ wue_WRIT6W$TRT + wue_WRIT6W$B)
anova(it6wue_analysis) 

it8wue_analysis <- lm (wue_WRIT8W$WUE ~ wue_WRIT8W$TRT + wue_WRIT8W$B)
anova(it8wue_analysis) 
### Results are similar with table1 in paper 


## ----table2-------------------------------------------------------------------
# Change the data colnames 
colnames(cwd)[4] <- "CWD"

# Read data for two groups WRIT6w and WRIT8w 
cwd_WRIT6W <- cwd[cwd$T == "WRIT6w", ]

cwd_WRIT8W <- cwd[cwd$T == "WRIT8w", ]


## -----------------------------------------------------------------------------
CWD6_TRT <- tapply(cwd_WRIT6W$CWD, cwd_WRIT6W$TRT, mean)
CWD6_TRT
CWD8_TRT <- tapply(cwd_WRIT8W$CWD, cwd_WRIT8W$TRT, mean)
CWD8_TRT


## ----figure1------------------------------------------------------------------
# average yield for each irrigation treatment
yield6_TRT <- tapply(ty_WRIT6W$TY, ty_WRIT6W$TRT, mean)

# standard deviation
yield6_TRT_sd <- tapply(ty_WRIT6W$TY, ty_WRIT6W$TRT, sd)

# number of observation per irrigation treatment 
yield6_TRT_n <- tapply(ty_WRIT6W$TY, ty_WRIT6W$TRT, length)

# standard error
yield6_TRT_sem <- yield6_TRT_sd / sqrt(yield6_TRT_n)


## ----figure1b-----------------------------------------------------------------
# Mean
yield8_TRT <- tapply(ty_WRIT8W$TY,ty_WRIT8W$TRT, mean)

# Standard error
yield8_TRT_sem <- tapply(ty_WRIT8W$TY,ty_WRIT8W$TRT,sd)/sqrt(tapply(ty_WRIT8W$TY,ty_WRIT8W$TRT, length))


## ----figure1c-----------------------------------------------------------------
yield_TRT <- cbind(WRIT6w=yield6_TRT, WRIT8w=yield8_TRT)
yield_TRT_sem <- cbind(yield6_TRT_sem, yield8_TRT_sem)


## ----figure1d-----------------------------------------------------------------
mids <- barplot(t(yield_TRT), beside=T, ylab="Tuber Yield (t/ha)", 
          cex.names=0.8, las=1, ylim=c(0,30))

# Add the bottom line
box(bty="l")

# Add legend
legend("topright", fill = c("black", "grey"), c("WRIT6w", "WRIT8w"), horiz = F)

# Add standard erros to the plot (only the postive part)
arrows (mids, t(yield_TRT), mids, t(yield_TRT + yield_TRT_sem), code = 3, angle = 90, length = .1)


## -----------------------------------------------------------------------------
library(agricolae)
yield6TRT_aov <- aov(ty_WRIT6W$TY ~ ty_WRIT6W$TRT)
a <- HSD.test(yield6TRT_aov, trt = "ty_WRIT6W$TRT" )

yield8TRT_aov <- aov(ty_WRIT8W$TY ~ ty_WRIT8W$TRT)

# Method 1
b <- TukeyHSD(yield8TRT_aov)

# Method 2 
b1 <- HSD.test(yield8TRT_aov, trt = "ty_WRIT8W$TRT" )

### plot values are similar but not the siginificance. 


## ----table3-------------------------------------------------------------------

# In the table, there are observations for different time, need to reorgnize the data, to make time as a factor
rwc = reshape(rwc8[,-1], varying = c("10DAT", "28DAT", "39DAT", "52DAT"), idvar=c("B", "TRT"), direction="long", v.names="RWC", times=c(10, 28, 39, 52))

rwc$time <- as.factor(rwc$time)


## -----------------------------------------------------------------------------
model <- lm(RWC ~ block + TRT + time + TRT*time, data = rwc)
anova(model)


## ----figure2------------------------------------------------------------------

times <- c(10, 28, 39, 52)
lapply(times, function(i) anova(lm(RWC ~ TRT, data=rwc[rwc$time==i, ])))

f <- function(i) { 
  d <- rwc[rwc$time==i, ]
  anov <- aov(RWC ~ TRT, data=d)
  HSD.test(anov, trt = "TRT")$groups
}
lapply(times, f)


## ----agg----------------------------------------------------------------------
rwc_agg <- aggregate(rwc[, 'RWC', drop=FALSE], rwc[, c('time',"TRT")], function(i) cbind(mean(i), sd(i), length(i)))

rwc_agg <- cbind(rwc_agg[,-3], rwc_agg[,3])
colnames(rwc_agg)[3:5] <- c("mean", "sd", "n")
dim(rwc_agg)

# Standard error
rwc_agg$sem <- rwc_agg$sd / sqrt(rwc_agg$n)

# Make a qqplot this time
library(ggplot2)
qplot(data = rwc_agg, x = as.factor(time), y = mean, pch= TRT, xlab = "DAT", ylab = "RWC(%)", ylim = c(70,90))
#+ geom_errorbar(aes(ymin = rwc_agg$mean - rwc_agg$sem, ymax = rwc_agg$sem + rwc_agg$mean), colour="black", width=.05)