Including:
Load Libraries
library(ezknitr)
library(knitr)
library(metafor)
library(ggplot2)
Clear environment
remove(list=ls())
Document settings
opts_chunk$set(fig.width = 6, fig.height = 4)
load(file="data/output_data/data_cleaned.R")
Rust
rust.v.year <- rma.uni(yi = yi, vi = vi, data=rust.data.SMD[rust.data.SMD$studyYear!=2009&
rust.data.SMD$studyYear!=2010&
rust.data.SMD$studyYear!=2011,],
method = "REML",
mods = ~studyYear)
rust.v.year.cat <- rma.uni(yi = yi, vi = vi, data=rust.data.SMD[rust.data.SMD$studyYear!=2009&
rust.data.SMD$studyYear!=2010&
rust.data.SMD$studyYear!=2011,],
method = "REML",
mods = ~as.character(studyYear)-1)
summary(rust.v.year)
##
## Mixed-Effects Model (k = 324; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -518.7071 1037.4142 1043.4142 1054.7378 1043.4896
##
## tau^2 (estimated amount of residual heterogeneity): 0.6334 (SE = 0.1024)
## tau (square root of estimated tau^2 value): 0.7959
## I^2 (residual heterogeneity / unaccounted variability): 49.28%
## H^2 (unaccounted variability / sampling variability): 1.97
## R^2 (amount of heterogeneity accounted for): 12.24%
##
## Test for Residual Heterogeneity:
## QE(df = 322) = 643.7100, p-val < .0001
##
## Test of Moderators (coefficient(s) 2):
## QM(df = 1) = 16.1462, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -208.1511 51.4684 -4.0442 <.0001 -309.0274 -107.2749 ***
## studyYear 0.1030 0.0256 4.0182 <.0001 0.0528 0.1533 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yield
yield.v.year <- rma.uni(yi = yi, vi = vi, data=yield.data.SMD,
method = "REML",
mods = ~studyYear)
yield.v.year.cat <- rma.uni(yi = yi, vi = vi, data=yield.data.SMD,
method = "REML",
mods = ~as.character(studyYear)-1)
summary(yield.v.year)
##
## Mixed-Effects Model (k = 513; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -536.9279 1073.8559 1079.8559 1092.5650 1079.9032
##
## tau^2 (estimated amount of residual heterogeneity): 0.0000 (SE = 0.0338)
## tau (square root of estimated tau^2 value): 0.0006
## I^2 (residual heterogeneity / unaccounted variability): 0.00%
## H^2 (unaccounted variability / sampling variability): 1.00
## R^2 (amount of heterogeneity accounted for): 78.09%
##
## Test for Residual Heterogeneity:
## QE(df = 511) = 428.9720, p-val = 0.9965
##
## Test of Moderators (coefficient(s) 2):
## QM(df = 1) = 27.5830, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 130.4474 24.7084 5.2795 <.0001 82.0198 178.8750 ***
## studyYear -0.0646 0.0123 -5.2520 <.0001 -0.0887 -0.0405 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
100-seed-weight
seedwt.v.year <- rma.uni(yi = yi, vi = vi, data=seedwt.data.SMD,
method = "REML",
mods = ~studyYear)
seedwt.v.year.cat <- rma.uni(yi = yi, vi = vi, data=seedwt.data.SMD,
method = "REML",
mods = ~as.character(studyYear)-1)
summary(seedwt.v.year)
##
## Mixed-Effects Model (k = 207; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -177.8704 355.7407 361.7407 371.7098 361.8601
##
## tau^2 (estimated amount of residual heterogeneity): 0 (SE = 0.0527)
## tau (square root of estimated tau^2 value): 0
## I^2 (residual heterogeneity / unaccounted variability): 0.00%
## H^2 (unaccounted variability / sampling variability): 1.00
## R^2 (amount of heterogeneity accounted for): NA%
##
## Test for Residual Heterogeneity:
## QE(df = 205) = 105.1426, p-val = 1.0000
##
## Test of Moderators (coefficient(s) 2):
## QM(df = 1) = 15.5938, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 155.9732 39.3587 3.9629 <.0001 78.8316 233.1149 ***
## studyYear -0.0774 0.0196 -3.9489 <.0001 -0.1158 -0.0390 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Combine data for ggplot2
# Yield data
pred.yr.yield <- predict(yield.v.year, newmods = seq(from=2004, to=2014, by=0.1))
pred.yr.yield <- as.data.frame(cbind(seq(from=2004, to=2014, by=0.1),
pred.yr.yield$pred,
pred.yr.yield$ci.lb,
pred.yr.yield$ci.ub))
colnames(pred.yr.yield) <- c("x", "mean", "pLL", "pUL")
pred.yr.yield$dependent <- "Yield"
coef.yield.year.cat <- as.data.frame(cbind(yield.v.year.cat$b,
yield.v.year.cat$ci.lb,
yield.v.year.cat$ci.ub))
coef.yield.year.cat$dependent <- "Yield"
# Rust data
pred.yr.rust <- predict(rust.v.year, newmods = seq(from=2004, to=2014, by=0.1))
pred.yr.rust <- as.data.frame(cbind(seq(from=2004, to=2014, by=0.1),
pred.yr.rust$pred,
pred.yr.rust$ci.lb,
pred.yr.rust$ci.ub))
colnames(pred.yr.rust) <- c("x", "mean", "pLL", "pUL")
pred.yr.rust$dependent <- "Rust"
coef.rust.year.cat <- as.data.frame(cbind(rust.v.year.cat$b,
rust.v.year.cat$ci.lb,
rust.v.year.cat$ci.ub))
coef.rust.year.cat$dependent <- "Rust"
# 100-seed-weight
pred.yr.seedwt <- predict(seedwt.v.year, newmods = seq(from=2004, to=2014, by=0.1))
pred.yr.seedwt <- as.data.frame(cbind(seq(from=2004, to=2014, by=0.1),
pred.yr.seedwt$pred,
pred.yr.seedwt$ci.lb,
pred.yr.seedwt$ci.ub))
colnames(pred.yr.seedwt) <- c("x", "mean", "pLL", "pUL")
pred.yr.seedwt$dependent <- "100sw"
pred.app <- rbind(pred.yr.yield, pred.yr.rust, pred.yr.seedwt)
coef.seedwt.year.cat <- as.data.frame(cbind(seedwt.v.year.cat$b,
seedwt.v.year.cat$ci.lb,
seedwt.v.year.cat$ci.ub))
coef.seedwt.year.cat$dependent <- "100sw"
# Bind categorical data together
year.cat <- rbind(coef.seedwt.year.cat,
coef.rust.year.cat,
coef.yield.year.cat)
colnames(year.cat) <- c("mean", "LL", "UL", "dependent")
year.cat$Year[year.cat$dependent=="Yield"] <-
sort(unique(yield.data.SMD$studyYear))
year.cat$Year[year.cat$dependent=="Rust"] <-
sort(unique(rust.data.SMD$studyYear[rust.data.SMD$studyYear!=2009&
rust.data.SMD$studyYear!=2010&
rust.data.SMD$studyYear!=2011]))
year.cat$Year[year.cat$dependent=="100sw"] <-
sort(unique(seedwt.data.SMD$studyYear))
ggplot(data=pred.app, aes(x = x, y=mean))+
geom_line(aes(color=dependent))+
geom_ribbon(aes(ymin=pLL, ymax=pUL, fill=dependent), alpha=0.4)+
geom_point(data=year.cat,
aes(x=as.numeric(Year), y=mean, colour=dependent, shape=dependent))+
geom_hline(yintercept = 0, lty=2, color="grey")+
geom_errorbar(data=year.cat,
aes(x=Year, ymin=LL, ymax=UL, color=dependent) ,width=0.4)+
theme_bw()+
xlab("Study Year")+
ylab("Standardized Mean Difference")+
theme(legend.position="none")
Rust
# Model with all of the data points included
overall.std.rust <- rma.uni(yi = yi, vi = vi, data = rust.data.SMD, method="ML")
summary(overall.std.rust)
##
## Random-Effects Model (k = 324; tau^2 estimator: ML)
##
## logLik deviance AIC BIC AICc
## -528.6158 581.6922 1061.2317 1068.7932 1061.2691
##
## tau^2 (estimated amount of total heterogeneity): 0.7169 (SE = 0.1088)
## tau (square root of estimated tau^2 value): 0.8467
## I^2 (total heterogeneity / total variability): 52.39%
## H^2 (total variability / sampling variability): 2.10
##
## Test for Heterogeneity:
## Q(df = 323) = 680.9467, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -1.3496 0.0659 -20.4871 <.0001 -1.4787 -1.2205 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Testing the influence of each data point
inf.rust <- influence.rma.uni(overall.std.rust) # Test the influence
plot(inf.rust) # Plot the influence
rust.data.SMD[inf.rust$is.infl==TRUE,c(1:3,47:50)] # Records that are influential
## FID Reference ReferenceNumb growthStateClean scale m1i m2i
## 1 582 Sikora et al ref 74 74 2+ Scale 0-8 5.0 100
## 2 583 Sikora et al ref 74 74 2+ Scale 0-8 5.0 100
## 3 584 Sikora et al ref 74 74 2+ Scale 0-8 2.5 100
## 4 585 Sikora et al ref 74 74 2+ Scale 0-8 2.5 100
## 5 588 Sikora et al ref 74 74 2+ Scale 0-8 5.0 100
## 6 589 Sikora et al ref 74 74 2+ Scale 0-8 5.0 100
rust.is.infl <- rust.data.SMD$FID[inf.rust$is.infl==TRUE]
# Re-run analysis without the influential datapoints
overall.std.rust.noinfl <- rma.uni(yi = yi, vi = vi,
data = rust.data.SMD[!rust.data.SMD$FID %in% rust.is.infl,], method="ML")
summary(overall.std.rust.noinfl)
##
## Random-Effects Model (k = 318; tau^2 estimator: ML)
##
## logLik deviance AIC BIC AICc
## -496.5243 526.1320 997.0486 1004.5727 997.0867
##
## tau^2 (estimated amount of total heterogeneity): 0.4778 (SE = 0.0901)
## tau (square root of estimated tau^2 value): 0.6913
## I^2 (total heterogeneity / total variability): 42.31%
## H^2 (total variability / sampling variability): 1.73
##
## Test for Heterogeneity:
## Q(df = 317) = 572.4466, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -1.2544 0.0604 -20.7795 <.0001 -1.3727 -1.1361 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yield
# Model with all of the data points included
overall.std.yield <- rma.uni(yi = yi, vi = vi, data = yield.data.SMD, method="ML")
summary(overall.std.yield)
##
## Random-Effects Model (k = 513; tau^2 estimator: ML)
##
## logLik deviance AIC BIC AICc
## -551.9523 456.5556 1107.9046 1116.3852 1107.9282
##
## tau^2 (estimated amount of total heterogeneity): 0.0000 (SE = 0.0337)
## tau (square root of estimated tau^2 value): 0.0025
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 1.00
##
## Test for Heterogeneity:
## Q(df = 512) = 456.5550, p-val = 0.9623
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6800 0.0328 20.7501 <.0001 0.6158 0.7443 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Testing the influence of each data point
inf.yield <- influence.rma.uni(overall.std.yield) # Test the influence
plot(inf.yield) # Plot the influence
yield.data.SMD[inf.yield$is.infl==TRUE,c(1:3,48:51)] # Records that are influential
## FID Reference ReferenceNumb yield.kg.ha yieldCont.kg.ha
## 1 583 Sikora et al ref 74 74 6032.325 3517.175
## 2 584 Sikora et al ref 74 74 5548.125 3517.175
## 3 585 Sikora et al ref 74 74 5097.550 3517.175
## 4 586 Sikora et al ref 74 74 5400.175 3517.175
## 5 587 Sikora et al ref 74 74 5198.425 3517.175
## 6 588 Sikora et al ref 74 74 5655.725 3517.175
## 7 589 Sikora et al ref 74 74 6193.725 3517.175
## m1i sd1i
## 1 6032.325 1206.465
## 2 5548.125 1109.625
## 3 5097.550 1019.510
## 4 5400.175 1080.035
## 5 5198.425 1039.685
## 6 5655.725 1131.145
## 7 6193.725 1238.745
yield.is.infl <- yield.data.SMD$FID[inf.yield$is.infl==TRUE]
# Re-run analysis without the influential datapoints
overall.std.yield.noinfl <- rma.uni(yi = yi, vi = vi,
data = yield.data.SMD[!yield.data.SMD$FID %in% yield.is.infl,], method="ML")
summary(overall.std.yield.noinfl)
##
## Random-Effects Model (k = 506; tau^2 estimator: ML)
##
## logLik deviance AIC BIC AICc
## -520.9695 398.1567 1045.9390 1054.3921 1045.9629
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0348)
## tau (square root of estimated tau^2 value): 0
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 1.00
##
## Test for Heterogeneity:
## Q(df = 505) = 398.1567, p-val = 0.9998
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6377 0.0333 19.1759 <.0001 0.5725 0.7029 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
100-seed-weight
# Model with all of the data points included
overall.std.seedwt <- rma.uni(yi = yi, vi = vi, data = seedwt.data.SMD, method="ML")
summary(overall.std.seedwt)
##
## Random-Effects Model (k = 207; tau^2 estimator: ML)
##
## logLik deviance AIC BIC AICc
## -186.8635 120.7364 377.7269 384.3924 377.7858
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0524)
## tau (square root of estimated tau^2 value): 0
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 1.00
##
## Test for Heterogeneity:
## Q(df = 206) = 120.7364, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.5499 0.0509 10.8022 <.0001 0.4502 0.6497 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Testing the influence of each data point
inf.seedwt <- influence.rma.uni(overall.std.seedwt) # Test the influence
plot(inf.seedwt) # Plot the influence
seedwt.data.SMD[inf.seedwt$is.infl==TRUE,c(1:3,47:50)] # Records that are influential
## [1] FID Reference ReferenceNumb growthStateClean
## [5] m1i sd1i m2i
## <0 rows> (or 0-length row.names)
# No record were found to be influential!
Cercospora
# Model with all of the data points included
overall.std.cerco <- rma.uni(yi = yi, vi = vi, data = cerco.data.MD, method="ML")
summary(overall.std.cerco)
##
## Random-Effects Model (k = 40; tau^2 estimator: ML)
##
## logLik deviance AIC BIC AICc
## -154.8307 14.4324 313.6613 317.0391 313.9856
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 56.4276)
## tau (square root of estimated tau^2 value): 0
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 1.00
##
## Test for Heterogeneity:
## Q(df = 39) = 14.4324, p-val = 0.9999
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -9.6508 2.5198 -3.8300 0.0001 -14.5894 -4.7121 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Testing the influence of each data point
inf.cerco <- influence.rma.uni(overall.std.cerco) # Test the influence
plot(inf.cerco) # Plot the influence
cerco.data.MD[inf.cerco$is.infl==TRUE,c(1:3,47:50)] # Records that are influential
## [1] FID Reference ReferenceNumb growthStateClean
## [5] scale m1i m2i
## <0 rows> (or 0-length row.names)
# No record were found to be influential!
Target Spot
# Model with all of the data points included
overall.std.target.spot <- rma.uni(yi = yi, vi = vi, data = target.spot.data.MD, method="ML")
summary(overall.std.target.spot)
##
## Random-Effects Model (k = 31; tau^2 estimator: ML)
##
## logLik deviance AIC BIC AICc
## -139.3531 57.4842 282.7062 285.5742 283.1348
##
## tau^2 (estimated amount of total heterogeneity): 269.9597 (SE = 119.3699)
## tau (square root of estimated tau^2 value): 16.4304
## I^2 (total heterogeneity / total variability): 57.44%
## H^2 (total variability / sampling variability): 2.35
##
## Test for Heterogeneity:
## Q(df = 30) = 72.8437, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -18.9129 3.8936 -4.8575 <.0001 -26.5442 -11.2816 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Testing the influence of each data point
inf.target.spot <- influence.rma.uni(overall.std.target.spot) # Test the influence
plot(inf.target.spot) # Plot the influence
target.spot.data.MD[inf.target.spot$is.infl==TRUE,c(1:3,47:50)] # Records that are influential
## [1] FID Reference ReferenceNumb growthStateClean
## [5] scale m1i m2i
## <0 rows> (or 0-length row.names)
target.spot.is.infl <- target.spot.data.MD$FID[inf.target.spot$is.infl==TRUE]
# No record were found to be influential!
Do analysis results change if random-effects model with study-level clustering is incorporated?
Rust
rust.pseudorep <- rma.mv(yi = yi, V = vi, data=rust.data.SMD,
random = ~1|Reference)
summary(rust.pseudorep)
##
## Multivariate Meta-Analysis Model (k = 324; method: REML)
##
## logLik Deviance AIC BIC AICc
## -370.8220 741.6441 745.6441 753.1994 745.6816
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 1.2148 1.1022 34 no Reference
##
## Test for Heterogeneity:
## Q(df = 323) = 680.9467, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -1.4105 0.1964 -7.1824 <.0001 -1.7954 -1.0256 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yield
yield.pseudorep <- rma.mv(yi = yi, V = vi, data=yield.data.SMD,
random = ~1|Reference)
summary(yield.pseudorep)
##
## Multivariate Meta-Analysis Model (k = 513; method: REML)
##
## logLik Deviance AIC BIC AICc
## -446.1894 892.3788 896.3788 904.8555 896.4024
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.3453 0.5876 57 no Reference
##
## Test for Heterogeneity:
## Q(df = 512) = 456.5550, p-val = 0.9623
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6693 0.0864 7.7434 <.0001 0.4999 0.8387 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
100-seed-weight
seedwt.pseudorep <- rma.mv(yi = yi, V = vi, data=seedwt.data.SMD,
random = ~1|Reference)
summary(seedwt.pseudorep)
##
## Multivariate Meta-Analysis Model (k = 207; method: REML)
##
## logLik Deviance AIC BIC AICc
## -166.5192 333.0384 337.0384 343.6941 337.0975
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.1965 0.4433 22 no Reference
##
## Test for Heterogeneity:
## Q(df = 206) = 120.7364, p-val = 1.0000
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.5165 0.1096 4.7136 <.0001 0.3017 0.7313 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Cercospora
cerco.pseudorep <- rma.mv(yi = yi, V = vi, data=cerco.data.MD,
random = ~1|Reference)
summary(cerco.pseudorep)
##
## Multivariate Meta-Analysis Model (k = 40; method: REML)
##
## logLik Deviance AIC BIC AICc
## -150.0619 300.1238 304.1238 307.4509 304.4571
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 43.1053 6.5655 4 no Reference
##
## Test for Heterogeneity:
## Q(df = 39) = 14.4324, p-val = 0.9999
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -9.1675 4.1492 -2.2095 0.0271 -17.2997 -1.0353 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Target spot
target.spot.pseudorep <- rma.mv(yi = yi, V = vi, data=target.spot.data.MD,
random = ~1|Reference)
summary(target.spot.pseudorep)
##
## Multivariate Meta-Analysis Model (k = 31; method: REML)
##
## logLik Deviance AIC BIC AICc
## -113.7684 227.5369 231.5369 234.3392 231.9813
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 509.8941 22.5808 4 no Reference
##
## Test for Heterogeneity:
## Q(df = 30) = 72.8437, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -15.7447 11.5871 -1.3588 0.1742 -38.4551 6.9656
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Spun with ezspin(“programs/sensitivity_analysis.R”, out_dir=“output”, fig_dir=“figures”, keep_md=FALSE)
Session Info:
sessionInfo()
## R version 3.3.0 (2016-05-03)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.10.5 (Yosemite)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_2.1.0 doBy_4.5-15 metafor_1.9-8 Matrix_1.2-6 knitr_1.13
## [6] ezknitr_0.4
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.5 lattice_0.20-33 R.methodsS3_1.7.1
## [4] MASS_7.3-45 plyr_1.8.4 grid_3.3.0
## [7] gtable_0.2.0 formatR_1.4 magrittr_1.5
## [10] scales_0.4.0 evaluate_0.9 stringi_1.1.1
## [13] R.oo_1.20.0 R.utils_2.3.0 labeling_0.3
## [16] tools_3.3.0 stringr_1.0.0 munsell_0.4.3
## [19] markdown_0.7.7 colorspace_1.2-6