## ----data---------------------------------------------------------------------
ff <- agro::get_data_from_uri("doi:10.34725/DVN/1BLSQY", ".")
ff


## ----sheets-------------------------------------------------------------------
readxl::excel_sheets(ff)


## ----readxl-------------------------------------------------------------------
grain <- readxl::read_excel(ff, sheet = "Grain")
straw <- readxl::read_excel(ff, sheet = "Straw")
dim(straw)
dim(grain)
colnames(grain)


## ----compare0-----------------------------------------------------------------
raw <- readxl::read_excel(ff, sheet = "Rawdata")
dim(raw)
dim(straw) + dim(grain)
table(raw$Component)


## ----compare1-----------------------------------------------------------------
grain2 <- raw[tolower(raw$Component) == "grain", ]


## ----compare2-----------------------------------------------------------------
grain  <- grain[order(grain$PubN, grain$Treatment), ]
grain2 <- grain2[order(grain2$PubN, grain2$Treatment), ]


## ----compare------------------------------------------------------------------
all(grain == grain2, na.rm=TRUE)


## ----fix----------------------------------------------------------------------
fixme <- function(x) {
  x$Technology[x$Technology == "Coppiced G. manure"] <- "Coppicing"
  x$Technology[x$Technology == "Rot-Ass"] <- "Rotations"
  x$Technology[x$Technology == "SWC"] <- "Soil water"
  x
}

grain <- fixme(grain)
straw <- fixme(straw)


## ----newvars------------------------------------------------------------------

newvars <- function(x) {
  x$D <- x$`Treatment yield` - x$`Control yield` 
  x$productivity = cut(x$`Control yield`, c(0, 0.5, 1, 1.5, 2, Inf))
}
grain <- newvars(grain)
straw <- newvars(straw)


## ----tab1B--------------------------------------------------------------------
make_table <- function(x) {
  # unique number of publications by technology
  N <- tapply(x$PubN, x$Technology, function(i) length(unique(i)))
  # observations by technology
  k <- table(x$Technology)
  # mean yield gain by technology
  mean <- tapply(x$D, x$Technology, mean)
  # quartiles of yield gain by technology
  q <- tapply(x$D, x$Technology, function(i) quantile(i, c(.75, .5, .25)))
  q <- do.call(cbind, q)
  iqr <- q[1,] - q[3,]
  ## count cases within classes
  # classify
  ycls <- cut(x$D, c(-Inf, 0, 1, 2, Inf))
  # tabulate
  z <- table(ycls, x$Technology)
  # compute percentage. Note that we need to transpose to divide
  # correctly, and then transpose back
  z <- 100 * t(t(z) / colSums(z))
  
  ## combine results
  tab <- rbind(N, k, mean, q, iqr, z)
  tab <- round(tab, 2)
  ## ordering
  tab <- tab[, c("Parkland", "Coppicing", "Green manure", "Mulching", "Rotations", "Soil water")]
  tab
}


## ----tab1C--------------------------------------------------------------------
tab1 <- make_table(grain)
tab1


## ----tab2A--------------------------------------------------------------------
tab2 <- make_table(straw)
tab2


## -----------------------------------------------------------------------------
# Check the crop types
sort(unique(grain$Crop))


## -----------------------------------------------------------------------------
### Select data for sorghum
grain$cropname <- ""
grain$cropname[substr(grain$Crop, 1, 7) == "Sorghum"] <- "Sorghum"
grain$cropname[grain$Crop %in% c("White sorghum", "Red sorghum")] <- "Sorghum"
grain$cropname[substr(grain$Crop, 1, 5) == "Maize"] <- "Maize"
grain$cropname[grain$Crop == "Millet"] <- "Millet"

# check by making a table
tab <- table(grain$Crop, grain$cropname)
head(tab, 15)


## -----------------------------------------------------------------------------
par(mfrow=c(1,3))
for (crop in c("Sorghum", "Millet", "Maize")) {
    plotmeans(D~Technology, data=grain, subset=cropname==crop, connect=FALSE, 
          xlab= "", ylab="Grain yield", las=2, ylim=c(-1, 1.6))
    text(1, 1.5, crop)
}



## ----xyz----------------------------------------------------------------------
g <- grain[grain$cropname != "", ]
g$productivity <- as.integer(g$productivity)
sciplot::lineplot.CI(productivity, D, group=cropname, data=g, 
          ci.fun=function(x) c(mean(x)-1.96 * sciplot::se(x), mean(x)+1.96 * sciplot::se(x)),
          xlab= "site productivity", ylab="Mean difference in yield", lwd=2)
abline(h=0, lwd=2, col="gray")



## -----------------------------------------------------------------------------

## Seperate data for grain yield of parkland based on rainfall class

The methods section of the paper states that "Rainfall (long-term average of total annual) of sites was also classified as low (<600 mm), medium (600-800 mm) and high (>800 mm)."


## ----rain---------------------------------------------------------------------
rain <- cut(grain_pt$Rainfall, c(0, 600, 800, Inf))
levels(rain) <- c("low", "medium", "high")

grain_pt$raincls <- rain