library(readr)
rf <- read_csv("~/Dr Nyagumbo/Chemical analysis/rfall.csv")
## Parsed with column specification:
## cols(
##   country = col_character(),
##   site = col_character(),
##   agroecoregion = col_character(),
##   treatment = col_character(),
##   association = col_character(),
##   cropping_system = col_character(),
##   tillage_practice = col_character(),
##   season = col_character(),
##   season_since = col_double(),
##   slope = col_double(),
##   rainfall = col_double(),
##   textureclass = col_character(),
##   maize_grain = col_double(),
##   rainfall_class = col_character(),
##   season_class = col_character()
## )
## Warning: 38 parsing failures.
##  row   col expected actual                                        file
## 3929 slope a double      - '~/Dr Nyagumbo/Chemical analysis/rfall.csv'
## 3930 slope a double      - '~/Dr Nyagumbo/Chemical analysis/rfall.csv'
## 3932 slope a double      - '~/Dr Nyagumbo/Chemical analysis/rfall.csv'
## 3933 slope a double      - '~/Dr Nyagumbo/Chemical analysis/rfall.csv'
## 3938 slope a double      - '~/Dr Nyagumbo/Chemical analysis/rfall.csv'
## .... ..... ........ ...... ...........................................
## See problems(...) for more details.
View(rf)
attach(rf)
names(rf)
##  [1] "country"          "site"             "agroecoregion"   
##  [4] "treatment"        "association"      "cropping_system" 
##  [7] "tillage_practice" "season"           "season_since"    
## [10] "slope"            "rainfall"         "textureclass"    
## [13] "maize_grain"      "rainfall_class"   "season_class"
country=as.factor(country)
site=as.factor(site)
agroecoregion=as.factor(agroecoregion)
treatment=as.factor(treatment)
#treatment=c(treatment)
system=factor(treatment,levels = c("Conv","CA",order=TRUE))
cropping_system=as.factor(cropping_system)
#cropping_system=c(cropping_system)
cropping_systems=factor(cropping_system,levels = c("Conv","CA_sole","CA_intercrop","CA_rotation"),ordered = TRUE)
#season=as.factor(season)
#season_since=as.factor(season_since)
season_class=as.factor(season_class)
#season_class=c(season_class)
season_clas=factor(season_class,levels = c("(0-2)","(3-5)","(>5)"),ordered = TRUE)
rainfall_class=as.factor(rainfall_class)
#rainfall_class=c(rainfall_class)
rainfall_clas=factor(rainfall_class,levels = c("Less_than_700","Between_700-1300","Greater_than_1300"),ordered = TRUE)
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)
### linear model
mod<-lm(maize_grain~cropping_system+rainfall_class+agroecoregion+season_class+cropping_system:rainfall_class+cropping_system*season_class+rainfall_class*season_class,data=rf)
summary(mod)
## 
## Call:
## lm(formula = maize_grain ~ cropping_system + rainfall_class + 
##     agroecoregion + season_class + cropping_system:rainfall_class + 
##     cropping_system * season_class + rainfall_class * season_class, 
##     data = rf)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3782.1 -1339.1  -322.0   947.1  9350.9 
## 
## Coefficients: (1 not defined because of singularities)
##                                                            Estimate
## (Intercept)                                                 3186.04
## cropping_systemCA_rotation                                  -573.00
## cropping_systemCA_sole                                      -363.35
## cropping_systemConv                                         -377.43
## rainfall_classGreater_than_1300                             3198.69
## rainfall_classLess_than_700                                -1217.25
## agroecoregionLowland                                        -535.94
## season_class(0-2)                                            100.39
## season_class(3-5)                                            256.63
## cropping_systemCA_rotation:rainfall_classGreater_than_1300  -256.62
## cropping_systemCA_sole:rainfall_classGreater_than_1300      -809.87
## cropping_systemConv:rainfall_classGreater_than_1300           45.43
## cropping_systemCA_rotation:rainfall_classLess_than_700      2030.31
## cropping_systemCA_sole:rainfall_classLess_than_700          1624.49
## cropping_systemConv:rainfall_classLess_than_700              269.53
## cropping_systemCA_rotation:season_class(0-2)                 780.21
## cropping_systemCA_sole:season_class(0-2)                     469.90
## cropping_systemConv:season_class(0-2)                        -89.96
## cropping_systemCA_rotation:season_class(3-5)                 930.78
## cropping_systemCA_sole:season_class(3-5)                     286.09
## cropping_systemConv:season_class(3-5)                        264.15
## rainfall_classGreater_than_1300:season_class(0-2)          -3325.74
## rainfall_classLess_than_700:season_class(0-2)                327.42
## rainfall_classGreater_than_1300:season_class(3-5)          -2897.71
## rainfall_classLess_than_700:season_class(3-5)                    NA
##                                                            Std. Error
## (Intercept)                                                    314.47
## cropping_systemCA_rotation                                     457.66
## cropping_systemCA_sole                                         384.56
## cropping_systemConv                                            433.65
## rainfall_classGreater_than_1300                               1912.15
## rainfall_classLess_than_700                                    210.57
## agroecoregionLowland                                            57.45
## season_class(0-2)                                              319.68
## season_class(3-5)                                              328.62
## cropping_systemCA_rotation:rainfall_classGreater_than_1300     406.56
## cropping_systemCA_sole:rainfall_classGreater_than_1300         417.89
## cropping_systemConv:rainfall_classGreater_than_1300            425.27
## cropping_systemCA_rotation:rainfall_classLess_than_700         328.56
## cropping_systemCA_sole:rainfall_classLess_than_700             224.99
## cropping_systemConv:rainfall_classLess_than_700                167.35
## cropping_systemCA_rotation:season_class(0-2)                   476.97
## cropping_systemCA_sole:season_class(0-2)                       403.81
## cropping_systemConv:season_class(0-2)                          444.47
## cropping_systemCA_rotation:season_class(3-5)                   483.83
## cropping_systemCA_sole:season_class(3-5)                       406.85
## cropping_systemConv:season_class(3-5)                          458.58
## rainfall_classGreater_than_1300:season_class(0-2)             1888.91
## rainfall_classLess_than_700:season_class(0-2)                  207.57
## rainfall_classGreater_than_1300:season_class(3-5)             1905.64
## rainfall_classLess_than_700:season_class(3-5)                      NA
##                                                            t value
## (Intercept)                                                 10.132
## cropping_systemCA_rotation                                  -1.252
## cropping_systemCA_sole                                      -0.945
## cropping_systemConv                                         -0.870
## rainfall_classGreater_than_1300                              1.673
## rainfall_classLess_than_700                                 -5.781
## agroecoregionLowland                                        -9.329
## season_class(0-2)                                            0.314
## season_class(3-5)                                            0.781
## cropping_systemCA_rotation:rainfall_classGreater_than_1300  -0.631
## cropping_systemCA_sole:rainfall_classGreater_than_1300      -1.938
## cropping_systemConv:rainfall_classGreater_than_1300          0.107
## cropping_systemCA_rotation:rainfall_classLess_than_700       6.179
## cropping_systemCA_sole:rainfall_classLess_than_700           7.220
## cropping_systemConv:rainfall_classLess_than_700              1.611
## cropping_systemCA_rotation:season_class(0-2)                 1.636
## cropping_systemCA_sole:season_class(0-2)                     1.164
## cropping_systemConv:season_class(0-2)                       -0.202
## cropping_systemCA_rotation:season_class(3-5)                 1.924
## cropping_systemCA_sole:season_class(3-5)                     0.703
## cropping_systemConv:season_class(3-5)                        0.576
## rainfall_classGreater_than_1300:season_class(0-2)           -1.761
## rainfall_classLess_than_700:season_class(0-2)                1.577
## rainfall_classGreater_than_1300:season_class(3-5)           -1.521
## rainfall_classLess_than_700:season_class(3-5)                   NA
##                                                            Pr(>|t|)    
## (Intercept)                                                 < 2e-16 ***
## cropping_systemCA_rotation                                   0.2106    
## cropping_systemCA_sole                                       0.3448    
## cropping_systemConv                                          0.3842    
## rainfall_classGreater_than_1300                              0.0944 .  
## rainfall_classLess_than_700                                7.96e-09 ***
## agroecoregionLowland                                        < 2e-16 ***
## season_class(0-2)                                            0.7535    
## season_class(3-5)                                            0.4349    
## cropping_systemCA_rotation:rainfall_classGreater_than_1300   0.5279    
## cropping_systemCA_sole:rainfall_classGreater_than_1300       0.0527 .  
## cropping_systemConv:rainfall_classGreater_than_1300          0.9149    
## cropping_systemCA_rotation:rainfall_classLess_than_700     7.02e-10 ***
## cropping_systemCA_sole:rainfall_classLess_than_700         6.08e-13 ***
## cropping_systemConv:rainfall_classLess_than_700              0.1073    
## cropping_systemCA_rotation:season_class(0-2)                 0.1020    
## cropping_systemCA_sole:season_class(0-2)                     0.2446    
## cropping_systemConv:season_class(0-2)                        0.8396    
## cropping_systemCA_rotation:season_class(3-5)                 0.0544 .  
## cropping_systemCA_sole:season_class(3-5)                     0.4820    
## cropping_systemConv:season_class(3-5)                        0.5646    
## rainfall_classGreater_than_1300:season_class(0-2)            0.0784 .  
## rainfall_classLess_than_700:season_class(0-2)                0.1148    
## rainfall_classGreater_than_1300:season_class(3-5)            0.1284    
## rainfall_classLess_than_700:season_class(3-5)                    NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1839 on 4360 degrees of freedom
## Multiple R-squared:  0.08431,    Adjusted R-squared:  0.07948 
## F-statistic: 17.45 on 23 and 4360 DF,  p-value: < 2.2e-16
######################################################################
###### general anova table
mod<-aov(maize_grain~cropping_system+rainfall_class+agroecoregion+season_class+cropping_system:rainfall_class+cropping_system*season_class+rainfall_class*season_class,data=rf)
anova(mod)
## Analysis of Variance Table
## 
## Response: maize_grain
##                                  Df     Sum Sq   Mean Sq F value    Pr(>F)
## cropping_system                   3 4.1994e+08 139979335 41.3914 < 2.2e-16
## rainfall_class                    2 1.7799e+08  88995900 26.3158 4.362e-12
## agroecoregion                     1 2.8280e+08 282795368 83.6215 < 2.2e-16
## season_class                      2 7.7543e+07  38771315 11.4645 1.082e-05
## cropping_system:rainfall_class    6 3.2127e+08  53544541 15.8329 < 2.2e-16
## cropping_system:season_class      6 5.3411e+07   8901912  2.6323   0.01503
## rainfall_class:season_class       3 2.4720e+07   8240091  2.4366   0.06280
## Residuals                      4360 1.4745e+10   3381849                  
##                                   
## cropping_system                ***
## rainfall_class                 ***
## agroecoregion                  ***
## season_class                   ***
## cropping_system:rainfall_class ***
## cropping_system:season_class   *  
## rainfall_class:season_class    .  
## Residuals                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##################################################################### end
###cropping systems
x<- LSD.test(mod,"cropping_system",alpha=0.01,group=FALSE)
diffograph(x,cex.axis=0.8,xlab="Maize Grain Yield [kg/ha]",ylab="Maize Grain Yield [kg/ha]")

h<-HSD.test(mod,"cropping_system",group =TRUE)
h
## $statistics
##   MSerror   Df     Mean       CV
##   3381849 4360 2811.201 65.41619
## 
## $parameters
##    test          name.t ntr StudentizedRange alpha
##   Tukey cropping_system   4         3.634561  0.05
## 
## $means
##              maize_grain      std    r Min   Max      Q25      Q50
## CA_intercrop    2810.941 1972.895 1292   3 11870 1302.145 2430.212
## CA_rotation     3285.613 2009.936  601 199 10923 1735.650 2906.688
## CA_sole         3087.620 1758.100  993 148 10822 1821.905 2766.608
## Conv            2437.857 1857.933 1498   4 11634 1050.000 2078.500
##                   Q75
## CA_intercrop 3830.126
## CA_rotation  4479.570
## CA_sole      4023.144
## Conv         3328.873
## 
## $comparison
## NULL
## 
## $groups
##              maize_grain groups
## CA_rotation     3285.613      a
## CA_sole         3087.620      a
## CA_intercrop    2810.941      b
## Conv            2437.857      c
## 
## attr(,"class")
## [1] "group"
bar.group(h$groups,ylim=c(0,4000),density=4,border="blue",ylab="Maize Grain Yield [kg/ha]")
######################################################################

############ rainfall
x<- LSD.test(mod,"rainfall_class",alpha=0.01,group=FALSE)
diffograph(x,cex.axis=0.8,xlab="Maize Grain Yield [kg/ha]",ylab="Maize Grain Yield [kg/ha]")

t<-HSD.test(mod,"rainfall_class",group =TRUE)
t
## $statistics
##   MSerror   Df     Mean       CV
##   3381849 4360 2811.201 65.41619
## 
## $parameters
##    test         name.t ntr StudentizedRange alpha
##   Tukey rainfall_class   3         3.315628  0.05
## 
## $means
##                   maize_grain      std    r      Min   Max      Q25
## Between_700-1300     2939.605 1968.668 3332  24.0000 11870 1460.573
## Greater_than_1300    2843.480 1579.183  180 304.4322  8187 1760.533
## Less_than_700        2313.891 1685.878  872   3.0000  9090  883.750
##                        Q50      Q75
## Between_700-1300  2553.702 3964.756
## Greater_than_1300 2476.800 3587.050
## Less_than_700     2065.097 3310.846
## 
## $comparison
## NULL
## 
## $groups
##                   maize_grain groups
## Between_700-1300     2939.605      a
## Greater_than_1300    2843.480      a
## Less_than_700        2313.891      b
## 
## attr(,"class")
## [1] "group"
bar.group(t$groups,ylim=c(0,4000),density=4,border="blue",ylab="Maize Grain Yield [kg/ha]")
#bar.err(t$means,variation="range",ylim=c(0,4000),bar=FALSE,col=0)

######################################################################
##### season 
x<- LSD.test(mod,"season_class",alpha=0.01,group=FALSE)
diffograph(x,cex.axis=0.8,xlab="Maize Grain Yield [kg/ha]",ylab="Maize Grain Yield [kg/ha]")

s<-HSD.test(mod,"season_class",group =TRUE)
s
## $statistics
##   MSerror   Df     Mean       CV
##   3381849 4360 2811.201 65.41619
## 
## $parameters
##    test       name.t ntr StudentizedRange alpha
##   Tukey season_class   3         3.315628  0.05
## 
## $means
##       maize_grain      std    r      Min       Max      Q25      Q50
## (>5)     2391.549 1770.472  169 465.8601  9792.077 1421.432 1859.060
## (0-2)    2686.666 1989.615 2821   3.0000 11870.000 1160.000 2249.636
## (3-5)    3114.094 1738.744 1394 100.0000 10277.999 1775.218 2832.032
##            Q75
## (>5)  2699.963
## (0-2) 3769.660
## (3-5) 4087.795
## 
## $comparison
## NULL
## 
## $groups
##       maize_grain groups
## (3-5)    3114.094      a
## (0-2)    2686.666      b
## (>5)     2391.549      b
## 
## attr(,"class")
## [1] "group"
bar.group(s$groups,ylim=c(0,4000),density=4,border="blue",ylab="Maize Grain Yield [kg/ha]")
#bar.err(t$means,variation="range",ylim=c(0,4000),bar=FALSE,col=0)

######################################################################

###### 
theme_set(theme_gray(base_size =10))
ggplot(rf, aes(x = season_clas, y = maize_grain, color = system)) +  
  geom_boxplot(size=1) + ylab("Maize Grain yield [kg/ha]") + xlab("Years")+theme(legend.position = c(0.75,0.85))

##############

theme_set(theme_gray(base_size =10))
ggplot(rf, aes(x = season_clas, y = maize_grain, color = system)) + geom_violin(scale = "area") +ylab("Maize Grain yield [kg/ha]") + xlab("Cropping Season [Years]")+theme(legend.position = c(0.75,0.85))

########################
mod<-aov(maize_grain~cropping_system+season_class+cropping_system*season_class,data=rf)
anova(mod)
## Analysis of Variance Table
## 
## Response: maize_grain
##                                Df     Sum Sq   Mean Sq F value    Pr(>F)
## cropping_system                 3 4.1994e+08 139979335 39.5834 < 2.2e-16
## season_class                    2 1.1496e+08  57480170 16.2543 9.268e-08
## cropping_system:season_class    6 1.0688e+08  17813335  5.0373 3.708e-05
## Residuals                    4372 1.5461e+10   3536310                  
##                                 
## cropping_system              ***
## season_class                 ***
## cropping_system:season_class ***
## Residuals                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#############

mod<-lm(maize_grain~cropping_system+season_class+cropping_system*season_class,data=rf)
summary(mod)
## 
## Call:
## lm(formula = maize_grain ~ cropping_system + season_class + cropping_system * 
##     season_class, data = rf)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3198.5 -1416.1  -363.4   956.9  9332.4 
## 
## Coefficients:
##                                              Estimate Std. Error t value
## (Intercept)                                  2741.972    317.864   8.626
## cropping_systemCA_rotation                   -664.879    467.883  -1.421
## cropping_systemCA_sole                       -406.505    393.215  -1.034
## cropping_systemConv                          -327.522    440.566  -0.743
## season_class(0-2)                               5.487    323.834   0.017
## season_class(3-5)                             251.626    334.104   0.753
## cropping_systemCA_rotation:season_class(0-2) 1152.108    483.561   2.383
## cropping_systemCA_sole:season_class(0-2)      869.413    408.620   2.128
## cropping_systemConv:season_class(0-2)        -118.379    448.294  -0.264
## cropping_systemCA_rotation:season_class(3-5) 1165.231    493.493   2.361
## cropping_systemCA_sole:season_class(3-5)      497.991    414.864   1.200
## cropping_systemConv:season_class(3-5)         310.769    465.438   0.668
##                                              Pr(>|t|)    
## (Intercept)                                    <2e-16 ***
## cropping_systemCA_rotation                     0.1554    
## cropping_systemCA_sole                         0.3013    
## cropping_systemConv                            0.4573    
## season_class(0-2)                              0.9865    
## season_class(3-5)                              0.4514    
## cropping_systemCA_rotation:season_class(0-2)   0.0172 *  
## cropping_systemCA_sole:season_class(0-2)       0.0334 *  
## cropping_systemConv:season_class(0-2)          0.7917    
## cropping_systemCA_rotation:season_class(3-5)   0.0183 *  
## cropping_systemCA_sole:season_class(3-5)       0.2301    
## cropping_systemConv:season_class(3-5)          0.5044    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1881 on 4372 degrees of freedom
## Multiple R-squared:  0.03986,    Adjusted R-squared:  0.03744 
## F-statistic:  16.5 on 11 and 4372 DF,  p-value: < 2.2e-16
######################

theme_set(theme_gray(base_size =10))
ggplot(rf, aes(x = season_clas, y = maize_grain, color = cropping_systems)) +  
  geom_boxplot(size=1) + ylab("Maize Grain yield [kg/ha]") + xlab("Years")+theme(legend.position = c(0.75,0.85))

###########################
theme_set(theme_gray(base_size =7))
ggplot(rf, aes(x = season_clas, y = maize_grain, color = system)) +  
  geom_boxplot(size=1) + ylab("Maize Grain yield [kg/ha]") + xlab("Cropping Season [Years]")+theme(legend.position = c(0.85,0.85))+facet_wrap(.~rainfall_class)

theme_set(theme_gray(base_size =10))
ggplot(rf, aes(x = rainfall_clas, y = maize_grain, color = system)) +  
  geom_boxplot(size=1) + ylab("Maize Grain yield [kg/ha]") + xlab("Seasonal Rainfall [mm]")+theme(legend.position = c(0.75,0.85))

#######################
theme_set(theme_gray(base_size =10))
ggplot(rf, aes(x = rainfall_clas, y = maize_grain, color = cropping_systems)) +  
  geom_boxplot(size=1) + ylab("Maize Grain yield [kg/ha]") + xlab("Seasonal Rainfall [mm]")+theme(legend.position = c(0.75,0.85))

ggplot(rf, aes(x = season_since , y = maize_grain, color = cropping_systems))+geom_point(method=glm,size=3,  aes(shape=cropping_systems))  + geom_smooth(position = "jitter", aes(fill=cropping_systems), level = 0.95)+ylab("Maize Grain yield [t/ha]") + xlab("Years")
## Warning: Ignoring unknown parameters: method
## Warning: Using shapes for an ordinal variable is not advised
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Computation failed in `stat_smooth()`:
## x has insufficient unique values to support 10 knots: reduce k.
## Warning: Removed 1 rows containing missing values (geom_point).

ggplot(rf, aes(season_since, maize_grain, shape=cropping_systems, color=cropping_systems)) + 
  geom_point(data =rf, size=1.5, aes(color=cropping_system, shape=cropping_system), position = "jitter") + 
  geom_smooth(method="lm", level = 0.95)+  ylim(0,8000) + 
  xlim(0,6)+ylab("Maize grain yield [kg/ha]") + xlab("Time since intervention  [years]")+
  theme(legend.position = c(0.90, 0.82))+
  geom_hline(yintercept =  mean(rf$maize_grain), linetype=2, color = "black", size=1)#+
## Warning: Removed 87 rows containing non-finite values (stat_smooth).
## Warning: Removed 665 rows containing missing values (geom_point).

  #geom_hline(yintercept = median(sim_summary$maize_grain), linetype=1, color = "black", size=1)

#+facet_wrap(agroecoregion ~ ., nc=5)
ggplot(rf, aes(x = season_since, y = maize_grain,fill=cropping_systems)) +  
  geom_point(data = rf, size=3, shape=20, position = "jitter") + 
  geom_smooth(method = "loess")+  ylim(0,8000) + xlim(0,6)+ylab("Maize grain yield [kg/ha]") + xlab("Cropping Season")+geom_vline(xintercept = 0, linetype=1, color = "black", size=1)
## Warning: Removed 87 rows containing non-finite values (stat_smooth).
## Warning: Removed 668 rows containing missing values (geom_point).

ggplot(rf, aes(x = season_since, y = maize_grain, color=cropping_systems)) +  
  geom_point(data = rf, size=2, shape=cropping_systems, position = "jitter") + 
  geom_smooth(method = "loess")+  ylim(0,8000) + xlim(0,6)+ylab("Maize grain yield [kg/ha]") + xlab("Cropping Season")+geom_vline(xintercept = 0, linetype=1, color = "black", size=1)
## Warning: Removed 87 rows containing non-finite values (stat_smooth).
## Warning: Removed 658 rows containing missing values (geom_point).

ggplot(rf, aes(x = season_since, y = maize_grain, color=cropping_systems)) +  
  geom_point(data = rf, size=2, shape=cropping_systems, position = "jitter") + 
  geom_smooth(method = "loess")+  ylim(0,8000) + xlim(0,6)+ylab("Maize grain yield [kg/ha]") + xlab("Cropping Season")+geom_vline(xintercept = 0, linetype=1, color = "black", size=1)+facet_wrap(.~agroecoregion)
## Warning: Removed 87 rows containing non-finite values (stat_smooth).
## Warning: Removed 678 rows containing missing values (geom_point).

ggplot(rf, aes(x = rainfall, y = maize_grain)) +  
  geom_point(data = rf, size=2, position = "jitter") + 
  geom_smooth(method = "loess")+  ylim(0,8000) + xlim(200,1800)+ylab("Maize grain yield [kg/ha]") + xlab("Seasonal Rainfall")+geom_vline(xintercept = 1, linetype=300, color = "black", size=1)
## Warning: Removed 86 rows containing non-finite values (stat_smooth).
## Warning: Removed 88 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_vline).

p <- ggplot(rf, aes(rainfall, maize_grain)) +
  geom_point()
# Add regression line
p +  geom_smooth(method = "loess")

p + geom_line(data = rf, aes(x = rainfall, y = maize_grain))

polynomial

library(tidyr)
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:Matrix':
## 
##     expand
ggplot(rf, aes(rainfall,maize_grain)) +
      geom_point() +stat_smooth(method = "lm", formula = y ~ x + I(x^2)+I(x^3)+I(x^4+I(x^5))+I(x^6), size = 1)

theme_set(theme_gray(base_size = 10)) ###This sets the font sizes of anything writt
m<-ggplot(rf, aes(x = rainfall_clas, y = maize_grain,fill=system))+ geom_boxplot(size=1.5,varwidth = TRUE,outlier.colour = "red",outlier.shape = 1, shape=6,fill=c("grey","grey","grey","grey","grey","grey")) + geom_smooth(method=lm)+ ylab("Maize Grain yield [kg/ha]") + xlab("Seasonal Rainfall [mm]")+facet_wrap(.~agroecoregion)
m

theme_set(theme_gray(base_size = 10)) ###This sets the font sizes of anything writt
m<-ggplot(rf, aes(x = rainfall_clas, y = maize_grain,color=system))+ geom_boxplot(size=1,varwidth = TRUE,outlier.colour = "red",outlier.shape = 1, shape=6,fill=c("grey","grey","grey","grey","grey","grey","grey","grey","grey","grey","grey","grey")) + geom_smooth(method=lm)+ ylab("Maize Grain yield [kg/ha]") + xlab("Seasonal Rainfall [mm]")+theme(legend.position = c(0.75,0.85))+facet_grid(.~agroecoregion)
m

p + geom_smooth(method = “loess”)

+geom_point(method=lm,size=3, aes(shape=cropping_systems))