Rust severity as dependent variable
Load Libraries
library(ezknitr)
library(knitr)
library(metafor)
library(devtools)
Clear environment and set seed
remove(list=ls())
set.seed(18837)
load(file="data/output_data/data_cleaned.R")
remove(list=c(#"rust.data.ROM",
"seedwt.data.ROM",
"yield.data.ROM"
))
Number of bootstrap simulations:
nsims <- 5000
Empty data frame to hold results
results.rust <- as.data.frame(matrix(NA,ncol=24,nrow=nsims))
colnames(results.rust) <- c("OVERALL","tau2",
"FLUT","MIXED","PYR","TEBU",
"Strobilurin","Triaz_Strob","Triazole",
"R1+", "R2+",
"R3","R5",
"Application Intrcpt","Application Slope",
"AZO_PROP",
"high","low","medium",
"Year Intrcpt|2004","Year Slope",
"2006","2007","2013")
Running the bootstraps
for(ii in 1:nsims){
studyIDS <- data.frame(ReferenceNumb = sample(rust.data.ROM$ReferenceNumb,
length(unique(rust.data.ROM$ReferenceNumb)),
replace=T))
newdata <- merge(studyIDS, rust.data.ROM)
# Overall analysis
meta <- rma.uni(yi = yi,
vi = 1/m2i,
data = newdata,
method = "REML")
results.rust$OVERALL[ii] <- meta$b
results.rust$tau2[ii] <- meta$tau2
# Active Ingredients
ai <- rma.uni(yi = yi,
vi = 1/m2i,
data = newdata,
method = "REML",
mods = ~category_ai-1)
results.rust$FLUT[ii] <- ai$b[rownames(ai$b)=="category_aiFLUT"]
results.rust$MIXED[ii] <- ai$b[rownames(ai$b)=="category_aiMIXED"]
results.rust$PYR[ii] <- ai$b[rownames(ai$b)=="category_aiPYR"]
results.rust$TEBU[ii] <- ai$b[rownames(ai$b)=="category_aiTEBU"]
# Classes
class <- rma.uni(yi = yi,
vi = 1/m2i,
data = newdata,
method = "REML",
mods = ~category_class-1)
results.rust$Strobilurin[ii] <-
class$b[rownames(class$b)=="category_classstrobilurin"]
results.rust$Triaz_Strob[ii] <-
class$b[rownames(class$b)=="category_classtriaz + strob"]
results.rust$Triazole[ii] <-
class$b[rownames(class$b)=="category_classtriazole"]
# Growth stage
rstage <- rma.uni(yi = yi,
vi = 1/m2i,
data = newdata,
method = "REML",
mods = ~category_rstage-1)
results.rust$`R1+`[ii] <-
rstage$b[rownames(rstage$b)=="category_rstage1+"]
results.rust$`R2+`[ii] <-
rstage$b[rownames(rstage$b)=="category_rstage2+"]
results.rust$R3[ii] <-
rstage$b[rownames(rstage$b)=="category_rstage3"]
results.rust$R5[ii] <-
rstage$b[rownames(rstage$b)=="category_rstage5"]
# Disease pressure
pressure <- rma.uni(yi = yi,
vi = (n1i+n2i)/(n1i*n2i),
data = newdata,
method = "REML",
mods = ~category_pressure-1)
results.rust$low[ii] <- pressure$b[rownames(pressure$b)=="category_pressurelow"]
if(nrow(newdata[newdata$category_pressure=="medium",])>0){
results.rust$medium[ii] <- pressure$b[rownames(pressure$b)=="category_pressuremedium"]
}
results.rust$high[ii] <- pressure$b[rownames(pressure$b)=="category_pressurehigh"]
# Number of applications
applications <- rma.uni(yi = yi,
vi = 1/m2i,
data = newdata,
method = "REML",
mods = ~number_applications)
results.rust$`Application Intrcpt`[ii] <-
applications$b[rownames(applications$b)=="intrcpt"]
results.rust$`Application Slope`[ii] <-
applications$b[rownames(applications$b)=="number_applications"]
# Number of applications as category
apps.categ <- rma.uni(yi = yi,
vi = 1/m2i,
data = newdata,
method = "REML",
mods = ~as.character(number_applications)-1)
results.rust$`1 Application`[ii] <-
apps.categ$b[rownames(apps.categ$b)=="as.character(number_applications)1"]
results.rust$`2 Applications`[ii] <-
apps.categ$b[rownames(apps.categ$b)=="as.character(number_applications)2"]
# Study year
# Condition on year1 = 2005
newdata$category_yearc <- newdata$category_year - 2004
year <- rma.uni(yi = yi,
vi = (n1i+n2i)/(n1i*n2i),
data = newdata,
method = "REML",
mods = ~category_yearc)
results.rust$`Year Intrcpt|2004`[ii] <- year$b[rownames(year$b)=="intrcpt"]
results.rust$`Year Slope`[ii] <- year$b[rownames(year$b)=="category_yearc"]
# Study year as category
year.categ <- rma.uni(yi = yi,
vi = (n1i+n2i)/(n1i*n2i),
data = newdata,
method = "REML",
mods = ~as.character(category_year)-1)
if("2006" %in% unique(newdata$category_year)){
results.rust$`2006`[ii] <-
year.categ$b[rownames(year.categ$b)=="as.character(category_year)2006"]
}
if("2007" %in% unique(newdata$category_year)){
results.rust$`2007`[ii] <-
year.categ$b[rownames(year.categ$b)=="as.character(category_year)2007"]
}
if("2013" %in% unique(newdata$category_year)){
results.rust$`2013`[ii] <-
year.categ$b[rownames(year.categ$b)=="as.character(category_year)2013"]
}
if(length(unique(newdata$category_year))==2 &
unique(newdata$category_year[!is.na(newdata$category_year)]=="2006")){
year.categ <- rma.uni(yi = yi,
vi = (n1i+n2i)/(n1i*n2i),
data = newdata[!is.na(newdata$category_year),],
method = "REML")
results.rust$`2006`[ii] <- year.categ$b
}
# Mixed active ingredients
if("AZO + PROP" %in% newdata$alphaIngred){
mixed <- rma.uni(yi = yi,
vi = 1/m2i,
data = newdata[newdata$alphaIngred=="AZO + PROP",],
method = "REML")
results.rust$`AZO_PROP`[ii] <- mixed$b
}
}
save(results.rust, file="data/output_data/results_rust.R")
Spun with ezspin(“programs/data_analysis_bootstrap_rust.R”, out_dir=“output”, fig_dir=“figures”, keep_md=FALSE)
Session Info:
devtools::session_info()
## Session info --------------------------------------------------------------
## setting value
## version R version 3.3.2 (2016-10-31)
## system x86_64, darwin13.4.0
## ui RStudio (0.99.902)
## language (EN)
## collate en_US.UTF-8
## tz America/Chicago
## date 2017-01-25
## Packages ------------------------------------------------------------------
## package * version date source
## devtools * 1.12.0 2016-06-24 CRAN (R 3.3.0)
## digest 0.6.10 2016-08-02 CRAN (R 3.3.0)
## doBy * 4.5-15 2016-03-31 CRAN (R 3.3.0)
## evaluate 0.10 2016-10-11 CRAN (R 3.3.0)
## ezknitr * 0.6 2016-09-16 CRAN (R 3.3.0)
## knitr * 1.15.1 2016-11-22 CRAN (R 3.3.2)
## lattice 0.20-34 2016-09-06 CRAN (R 3.3.2)
## magrittr 1.5 2014-11-22 CRAN (R 3.3.0)
## markdown 0.7.7 2015-04-22 CRAN (R 3.3.0)
## MASS 7.3-45 2016-04-21 CRAN (R 3.3.2)
## Matrix * 1.2-7.1 2016-09-01 CRAN (R 3.3.2)
## memoise 1.0.0 2016-01-29 CRAN (R 3.3.0)
## metafor * 1.9-9 2016-09-25 CRAN (R 3.3.0)
## R.methodsS3 1.7.1 2016-02-16 CRAN (R 3.3.0)
## R.oo 1.21.0 2016-11-01 CRAN (R 3.3.0)
## R.utils 2.5.0 2016-11-07 CRAN (R 3.3.0)
## rstudioapi 0.6 2016-06-27 CRAN (R 3.3.0)
## stringi 1.1.2 2016-10-01 CRAN (R 3.3.0)
## stringr 1.1.0 2016-08-19 CRAN (R 3.3.0)
## withr 1.0.2 2016-06-20 CRAN (R 3.3.0)