library(readr)
refined_rainfall <- read_csv("~/Dr Mupangwa/refined_rainfall.csv")
## Parsed with column specification:
## cols(
## country = col_character(),
## site = col_character(),
## agroecoregion = col_character(),
## farmer = col_character(),
## `Farmer Code` = col_logical(),
## Rep = col_double(),
## original_treatment = col_character(),
## treatment = col_character(),
## association = col_character(),
## cropping_system = col_character(),
## tillage_practice = col_character(),
## season = col_character(),
## season_since = col_double(),
## rainfall = col_double(),
## maize_grain = col_double()
## )
View(refined_rainfall)
attach(refined_rainfall)
names(refined_rainfall)
## [1] "country" "site" "agroecoregion"
## [4] "farmer" "Farmer Code" "Rep"
## [7] "original_treatment" "treatment" "association"
## [10] "cropping_system" "tillage_practice" "season"
## [13] "season_since" "rainfall" "maize_grain"
library(ggplot2)
library(maps)
library(ggalt)
library(extrafontdb)
library(MASS)
library(pscl)
## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis
library(psych)
##
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
library(gridExtra)
library(repr) ### adjusting the length and width of your plot
library(beanplot)
library("devtools")
library("yarrr")
## Loading required package: jpeg
## Loading required package: BayesFactor
## Loading required package: coda
## Loading required package: Matrix
## ************
## Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
##
## Type BFManual() to open the manual.
## ************
## Loading required package: circlize
## ========================================
## circlize version 0.4.5
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: http://jokergoo.github.io/circlize_book/book/
##
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
## in R. Bioinformatics 2014.
## ========================================
## yarrr v0.1.5. Citation info at citation('yarrr'). Package guide at yarrr.guide()
## Email me at Nathaniel.D.Phillips.is@gmail.com
##
## Attaching package: 'yarrr'
## The following object is masked from 'package:ggplot2':
##
## diamonds
library(agricolae)
library(easynls)
library(MVN)
## sROC 0.1-2 loaded
library(lme4)
country=as.factor(country)
site=as.factor(site)
agroecoregion=as.factor(agroecoregion)
treatment=as.factor(treatment)
cropping_system=as.factor(cropping_system)
season=as.factor(season)
mod<-aov(maize_grain~cropping_system+rainfall+season_since+agroecoregion+cropping_system+rainfall:season_since+cropping_system:rainfall,data=refined_rainfall)
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## cropping_system 4 1.163e+08 2.908e+07 10.023 4.83e-08 ***
## rainfall 1 3.437e+07 3.437e+07 11.845 0.000588 ***
## season_since 1 1.143e+07 1.143e+07 3.940 0.047255 *
## agroecoregion 1 1.034e+09 1.034e+09 356.446 < 2e-16 ***
## rainfall:season_since 1 1.263e+07 1.263e+07 4.354 0.037038 *
## cropping_system:rainfall 4 1.847e+07 4.617e+06 1.591 0.173911
## Residuals 2333 6.770e+09 2.902e+06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod)
## Analysis of Variance Table
##
## Response: maize_grain
## Df Sum Sq Mean Sq F value Pr(>F)
## cropping_system 4 116334822 29083705 10.0230 4.833e-08 ***
## rainfall 1 34371126 34371126 11.8451 0.0005883 ***
## season_since 1 11434054 11434054 3.9405 0.0472550 *
## agroecoregion 1 1034304193 1034304193 356.4464 < 2.2e-16 ***
## rainfall:season_since 1 12633062 12633062 4.3537 0.0370378 *
## cropping_system:rainfall 4 18468489 4617122 1.5912 0.1739108
## Residuals 2333 6769689214 2901710
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
x<- LSD.test(mod,"cropping_system",alpha=0.01,group=FALSE)
diffograph(x,cex.axis=0.8,xlab="Maize Grain Yield [kg/ha]",ylab="")
# Example 2
x<- REGW.test(mod,"season_since",alpha=0.01,group=FALSE)
diffograph(x,cex.axis=0.6,xlab="Maize Grain Yield [kg/ha]",ylab="",color1="brown",color2="green")
mod<-lm(maize_grain~cropping_system+rainfall+season_since+cropping_system+rainfall:season_since+cropping_system:rainfall,data=refined_rainfall)
summary(mod)
##
## Call:
## lm(formula = maize_grain ~ cropping_system + rainfall + season_since +
## cropping_system + rainfall:season_since + cropping_system:rainfall,
## data = refined_rainfall)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3360.1 -1318.0 -350.0 934.7 8821.5
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 3462.8890 495.0114 6.996
## cropping_systemCA_rotation 1108.6767 606.3873 1.828
## cropping_systemCA_sole 1006.5353 543.1850 1.853
## cropping_systemConv_intercrop -3022.4398 1151.7974 -2.624
## cropping_systemConv_sole 1289.8175 600.3191 2.149
## rainfall -0.2740 0.5046 -0.543
## season_since -595.6446 117.9873 -5.048
## rainfall:season_since 0.5704 0.1205 4.735
## cropping_systemCA_rotation:rainfall -0.8868 0.6221 -1.426
## cropping_systemCA_sole:rainfall -0.9438 0.5627 -1.677
## cropping_systemConv_intercrop:rainfall 3.1379 1.5874 1.977
## cropping_systemConv_sole:rainfall -1.4117 0.6215 -2.271
## Pr(>|t|)
## (Intercept) 3.44e-12 ***
## cropping_systemCA_rotation 0.06763 .
## cropping_systemCA_sole 0.06400 .
## cropping_systemConv_intercrop 0.00874 **
## cropping_systemConv_sole 0.03177 *
## rainfall 0.58719
## season_since 4.80e-07 ***
## rainfall:season_since 2.32e-06 ***
## cropping_systemCA_rotation:rainfall 0.15410
## cropping_systemCA_sole:rainfall 0.09364 .
## cropping_systemConv_intercrop:rainfall 0.04819 *
## cropping_systemConv_sole:rainfall 0.02322 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1820 on 2334 degrees of freedom
## Multiple R-squared: 0.03289, Adjusted R-squared: 0.02833
## F-statistic: 7.215 on 11 and 2334 DF, p-value: 3.29e-12
anova(mod)
## Analysis of Variance Table
##
## Response: maize_grain
## Df Sum Sq Mean Sq F value Pr(>F)
## cropping_system 4 116334822 29083705 8.7767 4.962e-07 ***
## rainfall 1 34371126 34371126 10.3723 0.001297 **
## season_since 1 11434054 11434054 3.4505 0.063358 .
## rainfall:season_since 1 61683921 61683921 18.6147 1.666e-05 ***
## cropping_system:rainfall 4 39174496 9793624 2.9555 0.018926 *
## Residuals 2334 7734236540 3313726
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(refined_rainfall, aes(x =season_since,y = maize_grain, color = cropping_system))+ geom_smooth(method=lm, position = "jitter",aes(fill=cropping_system), level = 0.95)+ xlim(0,6)+ylab("Grain yield [kg/ha]") +xlab("Season")+theme(legend.position = c(0.75,0.85))
## Warning: Removed 4 rows containing missing values (geom_smooth).
ggplot(refined_rainfall, aes(x =season_since,y = maize_grain, color = cropping_system))+ geom_smooth(method=lm, position = "jitter",aes(fill=cropping_system), level = 0.95)+ xlim(0,6)+ylab("Grain yield [kg/ha]") +xlab("Season")+theme(legend.position = c(0.75,0.85))+facet_wrap(.~agroecoregion)
## Warning: Removed 6 rows containing missing values (geom_smooth).
ggplot(refined_rainfall, aes(x = season_since, y = maize_grain, color = cropping_system)) +geom_point(data = refined_rainfall, size=3, shape=20, position = "jitter")+ geom_smooth(method="lm", aes(fill=cropping_system))+ylab("Maize grain yield [kg/ha]") + xlab("Years since intervention")+ geom_vline(xintercept = 0, linetype=1, color = "black", size=1)+theme(legend.position = c(0.90, 0.90))
########## agro ecological regions
ggplot(refined_rainfall, aes(x = season_since, y = maize_grain, color = cropping_system)) + geom_smooth(method="lm", aes(fill=cropping_system))+ylab("Maize grain yield [kg/ha]") + xlab("Years since intervention")+ geom_vline(xintercept = 0, linetype=1, color = "black", size=1)+facet_wrap(.~agroecoregion)
#### by treatment
ggplot(refined_rainfall, aes(x = season_since, y = maize_grain, color = cropping_system)) + geom_smooth(method="lm", aes(fill=cropping_system))+ylab("Maize grain yield [kg/ha]") + xlab("Years since intervention")+ geom_vline(xintercept = 0, linetype=1, color = "black", size=1)+facet_wrap(.~treatment)
# Enhanced Scatterplot of MPG vs. Weight
# by Number of Car Cylinders
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:psych':
##
## logit
scatterplot(maize_grain~ rainfall | cropping_system, data=refined_rainfall,
xlab="Seasonal Rainfall", ylab="Maize grain yield [kg/ha]",
main="Maize grain yield [kg/ha]")
#labels=row.names(rf))
# High Density Scatterplot with Binning
library(hexbin)
## Warning: package 'hexbin' was built under R version 3.5.3
bin<-hexbin(rainfall, maize_grain, xbins=50)
plot(bin, main="Hexagonal Binning")
###################
library(psych)
scatter.hist(x=rainfall, y=maize_grain, density=TRUE, ellipse=TRUE,
xlab="Seasonal Rainfall", ylab="Maize grain yield [kg/ha]")
##overal
ggplot(refined_rainfall, aes(x = rainfall, y = maize_grain, color = cropping_system)) +geom_point(data = refined_rainfall, size=3, shape=20, position = "jitter")+ geom_smooth(method="lm", aes(fill=cropping_system))+ylab("Maize grain yield [kg/ha]") + xlab("Seasonal Ranfall")+ geom_vline(xintercept = 250, linetype=1, color = "black", size=1)+theme(legend.position = c(0.85, 0.90))
###agroecology
ggplot(refined_rainfall, aes(x = rainfall, y = maize_grain, color = cropping_system)) +geom_point(data = refined_rainfall, size=3, shape=20, position = "jitter")+ geom_smooth(method="lm", aes(fill=cropping_system))+ylab("Maize grain yield [kg/ha]") + xlab("Seasonal Ranfall")+theme(legend.position = c(0.60, 0.90))+facet_wrap(.~agroecoregion)
###treatment
ggplot(refined_rainfall, aes(x = rainfall, y = maize_grain, color = cropping_system)) +geom_point(data = refined_rainfall, size=3, shape=20, position = "jitter")+ geom_smooth(method="lm", aes(fill=cropping_system))+ylab("Maize grain yield [kg/ha]") + xlab("Seasonal Ranfall")+theme(legend.position = c(0.60, 0.90))+facet_wrap(.~treatment)
###
ggplot(refined_rainfall, aes(y = rainfall, x = season_since, color = cropping_system)) +geom_point(data = refined_rainfall, size=2, shape=cropping_system, position = "jitter")+ geom_smooth(method="lm", aes(fill=cropping_system))+ylab("Season Rainfall") + xlab("Season Since")+theme(legend.position = c(0.60, 0.90))
######
###
ggplot(refined_rainfall, aes(y = rainfall, x = season_since, color = cropping_system)) +geom_point(data = refined_rainfall, size=3, shape=20, position = "jitter")+ geom_smooth(method="lm", aes(fill=cropping_system))+ylab("Season Rainfall") + xlab("Season Since")+theme(legend.position = c(0.60, 0.90))+facet_wrap(.~agroecoregion)
d <- c(maize_grain,rainfall)
df <- data.frame(d, cropping_system)
ggplot(df) +
stat_density(aes(x = d, group = cropping_system), position = "stack", geom = "line", show.legend = F, color = "red") +
stat_density(aes(x = d, linetype = cropping_system,color=cropping_system), position = "identity", geom = "line")+xlab("Total seasonal rainfall(mm) and Maize grain yield [kg/ha] ")+ theme(legend.position = c(0.80, 0.80))
##### overal
d <- c(maize_grain,rainfall)
df <- data.frame(d, cropping_system)
ggplot(df) +
stat_density(aes(x = d), position = "stack", geom = "line", show.legend = F, color = "red") +
stat_density(aes(x = d), position = "identity", geom = "line")+xlab("Total seasonal rainfall(mm) and Maize grain yield [kg/ha] ")
