library(readr)
conso <- read_csv("~/Dr Nyagumbo/conso.csv")
## Parsed with column specification:
## cols(
##   country = col_character(),
##   site = col_character(),
##   agroecoregion = col_character(),
##   farmer = col_character(),
##   treatment = col_character(),
##   cropping_system = col_character(),
##   tillage_practice = col_character(),
##   season = col_character(),
##   season_since = col_double(),
##   rainfall = col_double(),
##   textural_class = col_character(),
##   Revised_class = col_character(),
##   sand_0_20_cm = col_double(),
##   silt_0_20_cm = col_double(),
##   `Clay 0_20_cm` = col_double(),
##   drainage_class = col_character(),
##   drainage_scale = col_double(),
##   `Mean slope%` = col_double(),
##   maize_grain = col_double()
## )
View(conso)
attach(conso)

defining string variables

country=as.factor(country)
site=as.factor(site)
agroecoregion=as.factor(agroecoregion)
treatment=as.factor(treatment)
Revised_class=as.factor(Revised_class)
cropping_system=c(cropping_system)
cropping_systems=factor(cropping_system,levels = c("Conv_sole","CA_sole","CA_intercrop","CA_rotation"),ordered = TRUE)
drainage_class=c(drainage_class)
drainage_clas=factor(drainage_class,levels = c("Well drained",
"Moderately well drained","Somewhat poorly drained","Poorly drained"),ordered=TRUE)
tillage_practice=as.factor(tillage_practice)
textural_class=as.factor(textural_class)
#sand_0_20_cm = col_double(),
#silt_0_20_cm = col_double(),
#`Clay 0_20_cm` = col_double(),
#drainage_class = col_character(),
#drainage_scale = col_double(),
#maize_grain = col_double()

Loading the required packages

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)

Performing the regression analysis : Combined analysis

full<-lm(maize_grain~agroecoregion+cropping_systems+drainage_class+rainfall+Revised_class+site+country,data = conso)
anova(full)
## Analysis of Variance Table
## 
## Response: maize_grain
##                    Df     Sum Sq   Mean Sq  F value    Pr(>F)    
## agroecoregion       1  976563658 976563658 440.8854 < 2.2e-16 ***
## cropping_systems    3   63071719  21023906   9.4916 3.138e-06 ***
## drainage_class      3  153338655  51112885  23.0757 1.038e-14 ***
## rainfall            1  297394288 297394288 134.2635 < 2.2e-16 ***
## Revised_class       4   35399292   8849823   3.9954  0.003112 ** 
## site                9 1118270234 124252248  56.0957 < 2.2e-16 ***
## Residuals        2248 4979332339   2215005                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(full)
## 
## Call:
## lm(formula = maize_grain ~ agroecoregion + cropping_systems + 
##     drainage_class + rainfall + Revised_class + site + country, 
##     data = conso)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4903.8  -943.6  -144.9   733.9  6836.3 
## 
## Coefficients: (2 not defined because of singularities)
##                                         Estimate Std. Error t value
## (Intercept)                            4886.6609   408.8256  11.953
## agroecoregionLowland                  -2157.6575   238.2118  -9.058
## cropping_systems.L                      272.5631    67.3746   4.045
## cropping_systems.Q                        1.3929    65.3313   0.021
## cropping_systems.C                       77.2160    64.8169   1.191
## drainage_classPoorly drained          -1962.2688   252.8911  -7.759
## drainage_classSomewhat poorly drained  -264.3601    98.3056  -2.689
## drainage_classWell drained               27.8028   114.7489   0.242
## rainfall                                 -0.6454     0.2791  -2.312
## Revised_classClay Loam                 -267.3666   316.2716  -0.845
## Revised_classLoamy Sand                -189.3863   318.0836  -0.595
## Revised_classSand                      -431.6701   336.4233  -1.283
## Revised_classSandy Loam                 -81.6148   300.5202  -0.272
## siteCabango                            -124.9340   207.4306  -0.602
## siteChipole                            -259.2159   208.1397  -1.245
## siteGorongosa                           -67.8990   183.1426  -0.371
## siteKasungu                             188.2829   206.1454   0.913
## siteManica                              213.3218   224.1811   0.952
## siteMchinji                           -1652.5231   228.9433  -7.218
## siteMitundu                                   NA         NA      NA
## siteNtcheu                             2890.7018   210.2231  13.751
## siteSalima                             1561.9971   204.8154   7.626
## siteSussundenga                           5.2997   206.5678   0.026
## countryMozambique                             NA         NA      NA
##                                       Pr(>|t|)    
## (Intercept)                            < 2e-16 ***
## agroecoregionLowland                   < 2e-16 ***
## cropping_systems.L                    5.40e-05 ***
## cropping_systems.Q                     0.98299    
## cropping_systems.C                     0.23366    
## drainage_classPoorly drained          1.29e-14 ***
## drainage_classSomewhat poorly drained  0.00722 ** 
## drainage_classWell drained             0.80858    
## rainfall                               0.02085 *  
## Revised_classClay Loam                 0.39799    
## Revised_classLoamy Sand                0.55164    
## Revised_classSand                      0.19958    
## Revised_classSandy Loam                0.78597    
## siteCabango                            0.54704    
## siteChipole                            0.21312    
## siteGorongosa                          0.71086    
## siteKasungu                            0.36116    
## siteManica                             0.34142    
## siteMchinji                           7.18e-13 ***
## siteMitundu                                 NA    
## siteNtcheu                             < 2e-16 ***
## siteSalima                            3.54e-14 ***
## siteSussundenga                        0.97953    
## countryMozambique                           NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1488 on 2248 degrees of freedom
## Multiple R-squared:  0.3468, Adjusted R-squared:  0.3407 
## F-statistic: 56.84 on 21 and 2248 DF,  p-value: < 2.2e-16
mod<-aov(maize_grain~agroecoregion+cropping_systems+drainage_class+rainfall+Revised_class,data = conso)

y<-HSD.test(mod,"drainage_class",group = TRUE)
y
## $statistics
##   MSerror   Df     Mean       CV
##   2701640 2257 3019.922 54.42745
## 
## $parameters
##    test         name.t ntr StudentizedRange alpha
##   Tukey drainage_class   4         3.635867  0.05
## 
## $means
##                         maize_grain       std   r       Min      Max
## Moderately well drained    3128.197 1971.2584 971  273.0000 10923.00
## Poorly drained             2055.385  675.0205  44 1192.0800  3875.78
## Somewhat poorly drained    3028.129 1653.2203 529  225.3206  8773.00
## Well drained               2927.586 1793.8029 726  131.7970 11870.00
##                              Q25      Q50      Q75
## Moderately well drained 1671.670 2659.704 4095.298
## Poorly drained          1547.532 1923.218 2354.411
## Somewhat poorly drained 1675.254 2856.000 4013.047
## Well drained            1731.861 2583.231 3774.895
## 
## $comparison
## NULL
## 
## $groups
##                         maize_grain groups
## Moderately well drained    3128.197      a
## Somewhat poorly drained    3028.129      a
## Well drained               2927.586      a
## Poorly drained             2055.385      b
## 
## attr(,"class")
## [1] "group"

Model with interactions

ful<-lm(maize_grain~agroecoregion+cropping_systems+drainage_class+rainfall+Revised_class+agroecoregion*cropping_systems+agroecoregion*drainage_class+agroecoregion*rainfall+agroecoregion*Revised_class+cropping_systems*drainage_class+cropping_systems*rainfall+cropping_systems*Revised_class+drainage_class*rainfall+drainage_class*Revised_class+rainfall*Revised_class+site*cropping_systems+site*drainage_class+site*rainfall+site*Revised_class ,data = conso)
anova(ful)
## Analysis of Variance Table
## 
## Response: maize_grain
##                                   Df     Sum Sq   Mean Sq  F value
## agroecoregion                      1  976563658 976563658 456.8820
## cropping_systems                   3   63071719  21023906   9.8360
## drainage_class                     3  153338655  51112885  23.9130
## rainfall                           1  297394288 297394288 139.1349
## Revised_class                      4   35399292   8849823   4.1404
## site                               9 1118270234 124252248  58.1310
## agroecoregion:cropping_systems     3    3271710   1090570   0.5102
## agroecoregion:drainage_class       2    2187082   1093541   0.5116
## agroecoregion:rainfall             1     135348    135348   0.0633
## agroecoregion:Revised_class        4   20432986   5108247   2.3899
## cropping_systems:drainage_class    9   15769099   1752122   0.8197
## cropping_systems:rainfall          3    1223537    407846   0.1908
## cropping_systems:Revised_class    12    9007325    750610   0.3512
## drainage_class:rainfall            3   35500445  11833482   5.5363
## drainage_class:Revised_class       6   15407260   2567877   1.2014
## rainfall:Revised_class             3   25438065   8479355   3.9670
## cropping_systems:site             24   82184580   3424357   1.6021
## drainage_class:site                7   64226225   9175175   4.2926
## rainfall:site                      9   78472076   8719120   4.0792
## Revised_class:site                13   32690866   2514682   1.1765
## Residuals                       2149 4593385733   2137453         
##                                    Pr(>F)    
## agroecoregion                   < 2.2e-16 ***
## cropping_systems                1.924e-06 ***
## drainage_class                  3.200e-15 ***
## rainfall                        < 2.2e-16 ***
## Revised_class                   0.0024106 ** 
## site                            < 2.2e-16 ***
## agroecoregion:cropping_systems  0.6752560    
## agroecoregion:drainage_class    0.5996028    
## agroecoregion:rainfall          0.8013442    
## agroecoregion:Revised_class     0.0488764 *  
## cropping_systems:drainage_class 0.5979442    
## cropping_systems:rainfall       0.9027014    
## cropping_systems:Revised_class  0.9791374    
## drainage_class:rainfall         0.0008750 ***
## drainage_class:Revised_class    0.3024932    
## rainfall:Revised_class          0.0078382 ** 
## cropping_systems:site           0.0322620 *  
## drainage_class:site             0.0001001 ***
## rainfall:site                   3.282e-05 ***
## Revised_class:site              0.2902545    
## Residuals                                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

agroecoregioncropping_systems agroecoregiondrainage_class agroecoregionrainfall agroecoregiontextural_class cropping_systemsdrainage_class cropping_systemsrainfall cropping_systemstextural_class drainage_classrainfall drainage_classtextural_class rainfalltextural_class

Gneralized mixed effect model

fm2 <- lm(maize_grain~agroecoregion+cropping_systems+drainage_class+Revised_class+site+rainfall,data =conso)
summary(fm2)
## 
## Call:
## lm(formula = maize_grain ~ agroecoregion + cropping_systems + 
##     drainage_class + Revised_class + site + rainfall, data = conso)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4903.8  -943.6  -144.9   733.9  6836.3 
## 
## Coefficients: (1 not defined because of singularities)
##                                         Estimate Std. Error t value
## (Intercept)                            4886.6609   408.8256  11.953
## agroecoregionLowland                  -2157.6575   238.2118  -9.058
## cropping_systems.L                      272.5631    67.3746   4.045
## cropping_systems.Q                        1.3929    65.3313   0.021
## cropping_systems.C                       77.2160    64.8169   1.191
## drainage_classPoorly drained          -1962.2688   252.8911  -7.759
## drainage_classSomewhat poorly drained  -264.3601    98.3056  -2.689
## drainage_classWell drained               27.8028   114.7489   0.242
## Revised_classClay Loam                 -267.3666   316.2716  -0.845
## Revised_classLoamy Sand                -189.3863   318.0836  -0.595
## Revised_classSand                      -431.6701   336.4233  -1.283
## Revised_classSandy Loam                 -81.6148   300.5202  -0.272
## siteCabango                            -124.9340   207.4306  -0.602
## siteChipole                            -259.2159   208.1397  -1.245
## siteGorongosa                           -67.8990   183.1426  -0.371
## siteKasungu                             188.2829   206.1454   0.913
## siteManica                              213.3218   224.1811   0.952
## siteMchinji                           -1652.5231   228.9433  -7.218
## siteMitundu                                   NA         NA      NA
## siteNtcheu                             2890.7018   210.2231  13.751
## siteSalima                             1561.9971   204.8154   7.626
## siteSussundenga                           5.2997   206.5678   0.026
## rainfall                                 -0.6454     0.2791  -2.312
##                                       Pr(>|t|)    
## (Intercept)                            < 2e-16 ***
## agroecoregionLowland                   < 2e-16 ***
## cropping_systems.L                    5.40e-05 ***
## cropping_systems.Q                     0.98299    
## cropping_systems.C                     0.23366    
## drainage_classPoorly drained          1.29e-14 ***
## drainage_classSomewhat poorly drained  0.00722 ** 
## drainage_classWell drained             0.80858    
## Revised_classClay Loam                 0.39799    
## Revised_classLoamy Sand                0.55164    
## Revised_classSand                      0.19958    
## Revised_classSandy Loam                0.78597    
## siteCabango                            0.54704    
## siteChipole                            0.21312    
## siteGorongosa                          0.71086    
## siteKasungu                            0.36116    
## siteManica                             0.34142    
## siteMchinji                           7.18e-13 ***
## siteMitundu                                 NA    
## siteNtcheu                             < 2e-16 ***
## siteSalima                            3.54e-14 ***
## siteSussundenga                        0.97953    
## rainfall                               0.02085 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1488 on 2248 degrees of freedom
## Multiple R-squared:  0.3468, Adjusted R-squared:  0.3407 
## F-statistic: 56.84 on 21 and 2248 DF,  p-value: < 2.2e-16
##vcov(fm2)
print(fm2)
## 
## Call:
## lm(formula = maize_grain ~ agroecoregion + cropping_systems + 
##     drainage_class + Revised_class + site + rainfall, data = conso)
## 
## Coefficients:
##                           (Intercept)  
##                             4886.6609  
##                  agroecoregionLowland  
##                            -2157.6575  
##                    cropping_systems.L  
##                              272.5631  
##                    cropping_systems.Q  
##                                1.3929  
##                    cropping_systems.C  
##                               77.2160  
##          drainage_classPoorly drained  
##                            -1962.2688  
## drainage_classSomewhat poorly drained  
##                             -264.3601  
##            drainage_classWell drained  
##                               27.8028  
##                Revised_classClay Loam  
##                             -267.3666  
##               Revised_classLoamy Sand  
##                             -189.3863  
##                     Revised_classSand  
##                             -431.6701  
##               Revised_classSandy Loam  
##                              -81.6148  
##                           siteCabango  
##                             -124.9340  
##                           siteChipole  
##                             -259.2159  
##                         siteGorongosa  
##                              -67.8990  
##                           siteKasungu  
##                              188.2829  
##                            siteManica  
##                              213.3218  
##                           siteMchinji  
##                            -1652.5231  
##                           siteMitundu  
##                                    NA  
##                            siteNtcheu  
##                             2890.7018  
##                            siteSalima  
##                             1561.9971  
##                       siteSussundenga  
##                                5.2997  
##                              rainfall  
##                               -0.6454

Graphics

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

#interaction.plot(maize_grain, rainfall, cropping_systems, fixed = TRUE)
theme_set(theme_gray(base_size =9))
ggplot(conso, aes(x = cropping_systems, y = maize_grain, colour = agroecoregion)) + 
    geom_boxplot(size=1.2,varwidth = TRUE) + 
    geom_point(data = conso, aes(y = mean(maize_grain))) +
    geom_line(data = conso, aes(y =mean(maize_grain), group = agroecoregion))+ylab("Maize Grain yield [kg/ha]") + 
  xlab("Cropping Systems")+theme(legend.position = c(0.2, 0.85))

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

theme_set(theme_gray(base_size =8))
ggplot(conso, aes(x = cropping_systems, y = maize_grain, colour = agroecoregion)) + geom_boxplot(size=1.2,varwidth = TRUE) + geom_point(data = conso, aes(y = mean(maize_grain))) +geom_line(data = conso, aes(y =mean(maize_grain), group = agroecoregion)) +ylab("Maize Grain yield [kg/ha]") + xlab("Cropping Systems")+facet_wrap(.~ drainage_class)+theme(legend.position = c(0.85, 0.85))

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

theme_set(theme_gray(base_size =9))
ggplot(conso, aes(x = drainage_clas, y = maize_grain, colour = agroecoregion)) + 
    geom_boxplot(size=1.2,varwidth = TRUE) + 
    geom_point(data = conso, aes(y = mean(maize_grain))) +
    geom_line(data = conso, aes(y =mean(maize_grain), group = agroecoregion))+ylab("Maize Grain yield [kg/ha]") + 
  xlab("Drainage Class") +theme(legend.position = c(0.85, 0.85))

############################################################################
theme_set(theme_gray(base_size =7))
y<-ggplot(conso, aes(x = drainage_clas, y = maize_grain, colour = agroecoregion)) + 
    geom_boxplot(size=1.4,varwidth = TRUE) + 
    geom_point(data = conso, aes(y = mean(maize_grain))) +
    geom_line(data = conso, aes(y =mean(maize_grain), group = agroecoregion))+ylab("Maize Grain yield [kg/ha]") + 
  xlab("Drainage Class")+facet_wrap(.~ country)+theme(legend.position = c(0.85, 0.85))
y+theme_set(theme_gray(base_size =6))+theme(legend.position = c(0.85, 0.85))

######################################################################
######################################################################
###drainage class
theme_set(theme_gray(base_size = 10)) ###This sets the font sizes of anything writt
m<-ggplot(conso, aes(x = drainage_clas, y = maize_grain))+ geom_boxplot(size=1.5,varwidth = TRUE,outlier.colour = "red",outlier.shape = 1, shape=6,fill=c("grey","grey","grey","grey")) + geom_smooth(method=lm)+ ylab("Maize Grain yield [kg/ha]") + xlab("Drainage Class")
m

###### By country##################

theme_set(theme_gray(base_size =6)) ###This sets the font sizes of anything writt
m<-ggplot(conso, aes(x = drainage_clas, y = maize_grain))+ geom_boxplot(size=1.2,varwidth = TRUE,outlier.colour = "red",outlier.shape = 1, shape=6) + geom_smooth(method=lm)+ ylab("Maize Grain yield [kg/ha]") + xlab("Drainage Class")+facet_wrap(.~country)
m

##### wrapping by country and region
theme_set(theme_gray(base_size =6)) ###This sets the font sizes of anything writt
m<-ggplot(conso, aes(x = drainage_clas, y = maize_grain))+ geom_boxplot(size=1.2,varwidth = TRUE,outlier.colour = "red",outlier.shape = 1, shape=6) + geom_smooth(method=lm)+ ylab("Maize Grain yield [kg/ha]") + xlab("Drainage Class")+facet_wrap(country~agroecoregion)
m

############################### country and cropping system
theme_set(theme_gray(base_size =4))
m<-ggplot(conso, aes(x = drainage_clas, y = maize_grain))+ geom_boxplot(size=1,varwidth = TRUE,outlier.colour = "red",outlier.shape = 1, shape=6) + geom_smooth(method=lm)+ ylab("Maize Grain yield [kg/ha]") + xlab("Drainage Class")+facet_wrap(country~cropping_system)
m

########################## agroecology and cropping system
theme_set(theme_gray(base_size =4))
m<-ggplot(conso, aes(x = drainage_clas, y = maize_grain))+ geom_boxplot(size=1,varwidth = TRUE,outlier.colour = "red",outlier.shape = 1, shape=6) + geom_smooth(method=lm)+ ylab("Maize Grain yield [kg/ha]") + xlab("Drainage Class")+facet_wrap(agroecoregion~cropping_system)
m

##################### country by site 
theme_set(theme_gray(base_size =4))
m<-ggplot(conso, aes(x = drainage_clas, y = maize_grain))+ geom_boxplot(size=1,varwidth = TRUE,outlier.colour = "red",outlier.shape = 1, shape=6) + geom_smooth(method=lm)+ ylab("Maize Grain yield [kg/ha]") + xlab("Drainage Class")+facet_wrap(country~site)
m

################### site by cropping system 
theme_set(theme_gray(base_size =4))
m<-ggplot(conso, aes(x = drainage_clas, y = maize_grain))+ geom_boxplot(size=1,varwidth = TRUE,outlier.colour = "red",outlier.shape = 1, shape=6) + geom_smooth(method=lm)+ ylab("Maize Grain yield [kg/ha]") + xlab("Drainage Class")+facet_wrap(site~cropping_system)
m

Violin plots embeded with box plots

## violin plots of Country maize grain yields 
theme_set(theme_gray(base_size =10))
m <- ggplot(data=conso,aes(x=drainage_clas, y=maize_grain))
m + geom_violin(size=1.3,shape=8) + geom_boxplot(width=.3, outlier.size=0,fill=c("red","yellow","grey","green"))+ylab("Maize Grain yield [kg/ha]") + xlab("Drainage Class")
## Warning: Ignoring unknown parameters: shape

########################################################################
theme_set(theme_gray(base_size =10))
m <- ggplot(data=conso,aes(x=Revised_class, y=maize_grain))
m + geom_violin(size=1.1,shape=8) + geom_boxplot(width=.3, outlier.size=0,fill=c("red","yellow","grey","green","black"))+ylab("Maize Grain yield [kg/ha]") + xlab("Textural Class")
## Warning: Ignoring unknown parameters: shape

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

### by country
theme_set(theme_gray(base_size =10))
m <- ggplot(data=conso,aes(x=Revised_class, y=maize_grain))
m + geom_violin(size=1,shape=8) + geom_boxplot(width=.2, outlier.size=0,fill=c("red","yellow","grey","green","black","red","yellow","grey","green"))+ylab("Maize Grain yield [kg/ha]") + xlab("Textural Class")+facet_wrap(.~country)
## Warning: Ignoring unknown parameters: shape

### by agro ecology

theme_set(theme_gray(base_size =10))
m <- ggplot(data=conso,aes(x=Revised_class, y=maize_grain))
m + geom_violin(size=1,shape=8) + geom_boxplot(width=.2, outlier.size=0,fill=c("red","yellow","grey","green","black","red","yellow","grey","green","black"))+ylab("Maize Grain yield [kg/ha]") + xlab("Textural Class")+facet_wrap(.~agroecoregion)
## Warning: Ignoring unknown parameters: shape

#### by agroecology and country

theme_set(theme_gray(base_size =10))
m <- ggplot(data=conso,aes(x=Revised_class, y=maize_grain))
m + geom_violin(size=1,shape=8) + geom_boxplot(width=.2, outlier.size=0,fill=c("red","yellow","grey","green","black","red","yellow","grey","green","black","red","yellow","grey","green","black","red","yellow","grey"))+ylab("Maize Grain yield [kg/ha]") + xlab("Textural Class")+facet_wrap(country~agroecoregion)
## Warning: Ignoring unknown parameters: shape

#### country and site
theme_set(theme_gray(base_size =6))
m <- ggplot(data=conso,aes(x=Revised_class, y=maize_grain))
m + geom_violin(size=1,shape=8) + geom_boxplot(width=.2, outlier.size=0)+ylab("Maize Grain yield [kg/ha]") + xlab("Textural Class")+facet_wrap(country~site)
## Warning: Ignoring unknown parameters: shape

rainfall and maize grain yield by drainage classes

#rainfall and maize grain yield by drainage classes

ggplot(conso, aes(x = rainfall , y = maize_grain, color = drainage_clas)) +  
  geom_point(size=3,  aes(shape=drainage_clas)) + 
  geom_smooth(method=lm, position = "jitter", aes(fill=drainage_clas), level = 0.95)+ylab("Maize grain yield [kg/ha]") + 
  xlab("Total seasonal rainfall[mm]")+
  ylim(0,12000) + xlim(0,1400)+
  geom_abline(xintercept = 0, linetype=2, color = "red", size=1)+theme(legend.position = c(0.85, 0.85))
## Warning: Ignoring unknown parameters: xintercept
## Warning: Using shapes for an ordinal variable is not advised
## Warning: Removed 14 rows containing non-finite values (stat_smooth).
## Warning: Removed 14 rows containing missing values (geom_point).

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

one to one plots

######################################################################
######################################################################
##grain yield and drainage_scale
theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =drainage_scale, y =maize_grain ))+ geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Drainage Scale")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)
## Warning: Ignoring unknown parameters: xintercept

#######################################
## faceting by country and agroecology
theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =drainage_scale, y =maize_grain ))+ geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Drainage Scale")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)+facet_wrap(country~agroecoregion)
## Warning: Ignoring unknown parameters: xintercept

######################################
## facetting by country and cropping systems
theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =drainage_scale, y =maize_grain ))+ geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Drainage Scale")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)+theme(legend.position =c(0.45,0.85)) +facet_wrap(country~cropping_system)
## Warning: Ignoring unknown parameters: xintercept

############################################################################
############################################################################
##maize grain yield and sand_0_20_cm

theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =sand_0_20_cm, y =maize_grain )) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Sand %")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)
## Warning: Ignoring unknown parameters: xintercept

##################################
theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =sand_0_20_cm, y =maize_grain )) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Sand %")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)+facet_wrap(country~agroecoregion)
## Warning: Ignoring unknown parameters: xintercept

####################################
theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =sand_0_20_cm, y =maize_grain )) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Sand %")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)+facet_wrap(country~cropping_system)
## Warning: Ignoring unknown parameters: xintercept

##########################################################################
##########################################################################
### silt and maize grain yield 
theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =silt_0_20_cm, y =maize_grain )) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Silt %")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)
## Warning: Ignoring unknown parameters: xintercept

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

theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =silt_0_20_cm, y =maize_grain )) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Silt %")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)+facet_wrap(country~agroecoregion)
## Warning: Ignoring unknown parameters: xintercept

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

theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =silt_0_20_cm, y =maize_grain )) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Silt %")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)+facet_wrap(country~cropping_system)
## Warning: Ignoring unknown parameters: xintercept

#########################################################################
#########################################################################
## clay and grain yield accross countries 
theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =`Clay 0_20_cm`, y =maize_grain )) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Clay %")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)
## Warning: Ignoring unknown parameters: xintercept

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

theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =`Clay 0_20_cm`, y =maize_grain )) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Clay %")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)+facet_wrap(country~agroecoregion)
## Warning: Ignoring unknown parameters: xintercept

###############################
theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =`Clay 0_20_cm`, y =maize_grain )) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Maize grain yield [kg/ha]")+ xlab("Clay %")+geom_abline(xintercept = 0, linetype=2, color = "black", size=1)+facet_wrap(country~cropping_system)
## Warning: Ignoring unknown parameters: xintercept

data1=data.frame(rainfall, maize_grain)
nlsfit(data1, model=1)
## $Model
## [1] "y~a+b*x"
## 
## $Parameters
##                      maize_grain
## coefficient a         4436.31232
## coefficient b           -1.75137
## p-value t.test for a     0.00000
## p-value t.test for b     0.00000
## r-squared                0.02080
## adjusted r-squared       0.02030
## AIC                  40511.55303
## BIC                  40528.73563
nlsplot(data1, model=1)

nlsfit(data1, model=4)
## $Model
## [1] "y ~ (a + b * x + c * I(x^2)) * (x <= -0.5 * b/c) + (a + I(-b^2/(4 * c))) * (x > -0.5 * b/c)"
## 
## $Parameters
##                                  maize_grain
## coefficient a                   7.718575e+03
## coefficient b                  -1.012523e+01
## coefficient c                   5.185710e-03
## p-value t.test for a            0.000000e+00
## p-value t.test for b            6.683400e-04
## p-value t.test for c            1.119106e-02
## r-squared                       3.300000e-02
## adjusted r-squared              3.210000e-02
## AIC                             4.048498e+04
## BIC                             4.050789e+04
## maximum or minimum value for y  2.776136e+03
## critical point in x             9.762619e+02
nlsplot(data1, model=4)

nlsfit(data1, model=2)
## $Model
## [1] "y~a+b*x+c*x^2"
## 
## $Parameters
##                                  maize_grain
## coefficient a                   6.506163e+03
## coefficient b                  -6.809407e+00
## coefficient c                   2.985870e-03
## p-value t.test for a            0.000000e+00
## p-value t.test for b            2.000000e-08
## p-value t.test for c            1.918000e-05
## r-squared                       2.860000e-02
## adjusted r-squared              2.780000e-02
## AIC                             4.049526e+04
## BIC                             4.051817e+04
## maximum or minimum value for y  2.623882e+03
## critical point in x             1.140270e+03
nlsplot(data1, model=2)

ggplot(conso, aes(x = rainfall , y = maize_grain, color = agroecoregion)) +  
  geom_point(size=3,  aes(shape=agroecoregion)) + 
  geom_smooth(method=lm, position = "jitter", aes(fill=agroecoregion), level = 0.95)+ylab("Maize grain yield [kg/ha]") + 
  xlab("Total seasonal rainfall[mm]")+
  geom_abline(xintercept = 0, linetype=2, color = "red", size=1)+theme(legend.position = c(0.85, 0.85))+facet_wrap(.~drainage_class)
## Warning: Ignoring unknown parameters: xintercept

Multivariate Gaussan Distribution

set <- data.frame(rainfall, maize_grain)
result <- mvn(set, mvnTest = "hz", multivariatePlot = "persp")

result <- mvn(set, mvnTest = "hz", multivariatePlot = "contour")

result <- mvn(set,multivariatePlot = "persp")

bvn<-gibbs(10000,0.98) par(mfrow=c(3,2)) plot(bvn,col=1:10000,main=“bivariate normal distribution”,xlab=“X”,ylab=“Y”) plot(bvn,type=“l”,main=“bivariate normal distribution”,xlab=“X”,ylab=“Y”)

hist(bvn[,1],40,main=“bivariate normal distribution”,xlab=“X”,ylab=“”) hist(bvn[,2],40,main=“bivariate normal distribution”,xlab=“Y”,ylab=“”) par(mfrow=c(1,1))`

Maize grain yied and Rainfall

###
theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =rainfall, y =maize_grain))+ geom_smooth()+ylab("Maize grain yield [kg/ha]")+ xlab("Total Seasonal Rainfall [mm]")+xlim(400,1600)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 4 rows containing non-finite values (stat_smooth).

#######Maize grain yied and Rainfall by cropping system

theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =rainfall, y =maize_grain,color=cropping_systems))+ geom_smooth()+ylab("Maize grain yield [kg/ha]")+ xlab("Total Seasonal Rainfall [mm]")+xlim(400,1600)+theme(legend.position = c(0.85, 0.85))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 4 rows containing non-finite values (stat_smooth).

#####

theme_set(theme_gray(base_size =12))
ggplot(conso, aes(x =rainfall, y =maize_grain,color=drainage_clas))+ geom_smooth()+ylab("Maize grain yield [kg/ha]")+ xlab("Total Seasonal Rainfall [mm]")+xlim(400,1600)+theme(legend.position = c(0.70, 0.85))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 699.62
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 159.38
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 77498
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 699.62
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 159.38
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 77498

season Since
ggplot(conso, aes(x =season_since,y = maize_grain, color = cropping_system))+ geom_smooth(method=lm, position = "jitter",aes(fill=cropping_system), level = 0.95)+ ylim(0,10000) + xlim(0,6)+ylab("Grain yield  [kg/ha]") + 
xlab("Season")+theme(legend.position = c(0.2,0.85))+facet_wrap(.~agroecoregion)
## Warning: Removed 7 rows containing non-finite values (stat_smooth).
## Warning: Removed 3 rows containing missing values (geom_smooth).

ggplot(conso, aes(season_since, maize_grain, color=cropping_system, shape=cropping_system)) + geom_point(data = conso, size=1.5, aes(color=cropping_system, shape=cropping_system), position = "jitter") + geom_smooth(method="lm", level = 0.85,)+  ylim(0,7000) + 
 xlim(0,6)+ylab("Maize grain yield [kg/ha]") + xlab("Season Since")+theme(legend.position = c(0.5,0.85))+facet_wrap(.~agroecoregion)
## Warning: Removed 88 rows containing non-finite values (stat_smooth).
## Warning: Removed 291 rows containing missing values (geom_point).