library(ggplot2)
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")
## Loading required package: usethis
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.8
## 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)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
## sROC 0.1-2 loaded
library(lme4)
## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car
library(ggsignif)
library(ggpubr)
## Loading required package: magrittr
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
library(multcompView)
library(Rmisc)
## Loading required package: lattice
## Loading required package: plyr
## 
## Attaching package: 'plyr'
## The following object is masked from 'package:ggpubr':
## 
##     mutate
library(readr)
yld <- read_csv("~/Dr Nyagumbo/Chemical analysis/final yield with socio-analysis/yld.csv")
## Warning: Duplicated column names deduplicated: 'biomass' => 'biomass_1' [33]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   age = col_double(),
##   hhsize = col_double(),
##   experience_CA = col_double(),
##   weeding = col_double(),
##   row_spacing = col_double(),
##   plant_spacing = col_double(),
##   number_plants = col_double(),
##   plant_population = col_double(),
##   days_seeding = col_double(),
##   biomass = col_double(),
##   grain = col_double(),
##   biomass_1 = col_double()
## )
## See spec(...) for full column specifications.
View(yld)
attach(yld)
names(yld)
##  [1] "farmer_name"      "tech"             "season"           "district"        
##  [5] "village"          "site"             "agroecology"      "cropping_system" 
##  [9] "class"            "age"              "sex"              "hhsize"          
## [13] "experience_CA"    "local_soil_name"  "soil_colour"      "texture"         
## [17] "field_treatment"  "weeding"          "crop_name"        "variety"         
## [21] "row_spacing"      "plant_spacing"    "number_plants"    "tillage"         
## [25] "residue_applied"  "residue_type"     "plant_population" "planting_rains"  
## [29] "seeding"          "days_seeding"     "biomass"          "grain"           
## [33] "biomass_1"
tech=as.factor(tech)
season=as.factor(season)
district=as.factor(district)
agroecology=as.factor(agroecology)
site=as.factor(site)
#age
sex=as.factor(sex)
## hhsze
local_soil_name=as.factor(local_soil_name)
texture=as.factor(texture)
soil_colour=as.factor(soil_colour)
variet=factor(variety,levels = c("DKC 8033","MH33","SC627","Kanyani","Kanthochi","SC503","Mapasa","Syngenta","ZM523","MH18","Mphangala","Peacock","DKC 9089","Local","Masika","MH44","Nkango","Pannar","QPM","SC403"),ordered = TRUE)

###############################################
cropping_systems=factor(class,levels = c("CA","CONV","Intercrop maize","Monocrop maize","Mbeya fert","Normal fert"))
residue_applied=as.factor(residue_applied)
residue_type=as.factor(residue_type)
class=as.factor(class)
##############################################
biod<-lm(grain~tech+cropping_system+agroecology+residue_applied+plant_population+weeding+number_plants,data = yld)

summary(biod)
## 
## Call:
## lm(formula = grain ~ tech + cropping_system + agroecology + residue_applied + 
##     plant_population + weeding + number_plants, data = yld)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3213.5 -1415.4  -344.4   873.5  7142.2 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               2.840e+03  1.278e+03   2.222 0.028302 *  
## techNon-CSA              -6.962e+02  4.614e+02  -1.509 0.134135    
## cropping_systemIntercrop -1.332e+02  4.779e+02  -0.279 0.780949    
## cropping_systemMbeya      1.265e+03  6.444e+02   1.963 0.052106 .  
## agroecologyLow rainfall  -8.937e+01  4.194e+02  -0.213 0.831636    
## residue_appliedYes       -5.809e+01  5.280e+02  -0.110 0.912592    
## plant_population          6.963e-02  2.167e-02   3.213 0.001709 ** 
## weeding                  -7.879e+01  3.798e+02  -0.207 0.836030    
## number_plants            -1.171e+03  3.358e+02  -3.485 0.000701 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2128 on 113 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.2582, Adjusted R-squared:  0.2057 
## F-statistic: 4.917 on 8 and 113 DF,  p-value: 3.131e-05
anova(biod)
## Analysis of Variance Table
## 
## Response: grain
##                   Df    Sum Sq  Mean Sq F value    Pr(>F)    
## tech               1  25687973 25687973  5.6748 0.0188808 *  
## cropping_system    2  20505345 10252672  2.2649 0.1085321    
## agroecology        1   5852766  5852766  1.2929 0.2579121    
## residue_applied    1   6322492  6322492  1.3967 0.2397565    
## plant_population   1  64571763 64571763 14.2646 0.0002550 ***
## weeding            1    146504   146504  0.0324 0.8575527    
## number_plants      1  54990576 54990576 12.1480 0.0007005 ***
## Residuals        113 511517993  4526708                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model_1<-aov(grain~tech+agroecology+site+texture+soil_colour+variety+cropping_system+tech*agroecology+tech*site+tech*texture+tech*soil_colour+tech*cropping_system+tech*residue_applied+tech*cropping_system*agroecology+tech*cropping_system*site+tech*cropping_system*soil_colour+cropping_system*soil_colour*residue_applied+cropping_system*soil_colour*variety+residue_applied+residue_type+plant_population+weeding+number_plants+age+sex+hhsize+experience_CA,data = yld)

#############
out<-LSD.test(model_1,"tech",p.adj="hommel",console=TRUE)
## 
## Study: model_1 ~ "tech"
## 
## LSD t Test for grain 
## P value adjustment method: hommel 
## 
## Mean Square Error:  1623275 
## 
## tech,  means and individual ( 95 %) CI
## 
##            grain      std  r      LCL      UCL     Min      Max
## CSA     3903.330 2662.174 51 3534.267 4272.392 539.844 11091.00
## Non-CSA 2869.872 2127.760 51 2500.809 3238.934   0.000  8502.05
## 
## Alpha: 0.05 ; DF Error: 23
## Critical Value of t: 2.068658 
## 
## Minimum Significant Difference: 521.9328 
## 
## Treatments with the same letter are not significantly different.
## 
##            grain groups
## CSA     3903.330      a
## Non-CSA 2869.872      b
out
## $statistics
##   MSerror Df     Mean       CV  t.value      MSD
##   1623275 23 3386.601 37.62115 2.068658 521.9328
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD    hommel   tech   2  0.05
## 
## $means
##            grain      std  r      LCL      UCL     Min      Max      Q25
## CSA     3903.330 2662.174 51 3534.267 4272.392 539.844 11091.00 1817.968
## Non-CSA 2869.872 2127.760 51 2500.809 3238.934   0.000  8502.05 1421.559
##              Q50      Q75
## CSA     3103.922 5272.583
## Non-CSA 2028.926 4131.753
## 
## $comparison
## NULL
## 
## $groups
##            grain groups
## CSA     3903.330      a
## Non-CSA 2869.872      b
## 
## attr(,"class")
## [1] "group"
###################################################

out<-LSD.test(model_1,"cropping_system",p.adj="hommel",console=TRUE)
## 
## Study: model_1 ~ "cropping_system"
## 
## LSD t Test for grain 
## P value adjustment method: hommel 
## 
## Mean Square Error:  1623275 
## 
## cropping_system,  means and individual ( 95 %) CI
## 
##              grain      std  r      LCL      UCL     Min      Max
## CA        3236.625 2224.705 56 2884.424 3588.826   0.000 10082.23
## Intercrop 3179.784 2519.803 28 2681.696 3677.871 539.844 10186.92
## Mbeya     4174.907 2977.093 18 3553.682 4796.131 813.145 11091.00
## 
## Alpha: 0.05 ; DF Error: 23
## Critical Value of t: 2.068658 
## 
## Groups according to probability of means differences and alpha level( 0.05 )
## 
## Treatments with the same letter are not significantly different.
## 
##              grain groups
## Mbeya     4174.907      a
## CA        3236.625      b
## Intercrop 3179.784      b
out
## $statistics
##   MSerror Df     Mean       CV
##   1623275 23 3386.601 37.62115
## 
## $parameters
##         test p.ajusted          name.t ntr alpha
##   Fisher-LSD    hommel cropping_system   3  0.05
## 
## $means
##              grain      std  r      LCL      UCL     Min      Max      Q25
## CA        3236.625 2224.705 56 2884.424 3588.826   0.000 10082.23 1566.081
## Intercrop 3179.784 2519.803 28 2681.696 3677.871 539.844 10186.92 1681.100
## Mbeya     4174.907 2977.093 18 3553.682 4796.131 813.145 11091.00 1483.844
##                Q50      Q75
## CA        2791.701 4219.068
## Intercrop 2216.329 4155.100
## Mbeya     4294.150 6186.985
## 
## $comparison
## NULL
## 
## $groups
##              grain groups
## Mbeya     4174.907      a
## CA        3236.625      b
## Intercrop 3179.784      b
## 
## attr(,"class")
## [1] "group"
##################################################
md_social<-lm(grain~age+sex+hhsize+experience_CA,data = yld)
summary(md_social)
## 
## Call:
## lm(formula = grain ~ age + sex + hhsize + experience_CA, data = yld)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3341.4 -1680.6  -675.3  1074.8  6872.6 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     179.21    1495.05   0.120   0.9048  
## age              70.89      30.51   2.324   0.0222 *
## sexMale         414.13     568.48   0.728   0.4680  
## hhsize           13.13     123.84   0.106   0.9158  
## experience_CA   -78.41      31.49  -2.490   0.0144 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2349 on 99 degrees of freedom
##   (20 observations deleted due to missingness)
## Multiple R-squared:  0.1312, Adjusted R-squared:  0.09614 
## F-statistic: 3.739 on 4 and 99 DF,  p-value: 0.007086
anova(md_social)
## Analysis of Variance Table
## 
## Response: grain
##               Df    Sum Sq  Mean Sq F value  Pr(>F)  
## age            1  32457026 32457026  5.8811 0.01712 *
## sex            1  15804384 15804384  2.8637 0.09374 .
## hhsize         1     50810    50810  0.0092 0.92375  
## experience_CA  1  34224125 34224125  6.2013 0.01443 *
## Residuals     99 546363913  5518827                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### mean and sample size function

stat_box_data <- function(grain, upper_limit = max(grain)*1.5) {
  return( 
    data.frame(
      y =11000,
      label = paste('N =', length(grain), '\n',
                    'Mean =', round(mean(grain), 0), '\n')
    )
  )
}
#######################################################################
##tech in general

theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =tech, y = grain))+ geom_boxplot(size=0.7,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey")) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=3,
    color="blue"
  )+ theme_classic(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=12),
          axis.text.y = element_text(face="bold", color="black", 
                           size=11))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+
  geom_text(data = yld, x = 1, y = 4500, label = "a",color="blue",size=4)+
  geom_text(data = yld, x = 2, y = 3800, label = "b",color="blue",size=4)

agric<-subset(yld,agroecology=="Low rainfall")

ag<-aov(grain~tech,data=agric)

h<-HSD.test(model_1,"tech",group =TRUE)
h
## $statistics
##   MSerror Df     Mean       CV      MSD
##   1623275 23 3386.601 37.62115 521.9327
## 
## $parameters
##    test name.t ntr StudentizedRange alpha
##   Tukey   tech   2         2.925523  0.05
## 
## $means
##            grain      std  r     Min      Max      Q25      Q50      Q75
## CSA     3903.330 2662.174 51 539.844 11091.00 1817.968 3103.922 5272.583
## Non-CSA 2869.872 2127.760 51   0.000  8502.05 1421.559 2028.926 4131.753
## 
## $comparison
## NULL
## 
## $groups
##            grain groups
## CSA     3903.330      a
## Non-CSA 2869.872      b
## 
## attr(,"class")
## [1] "group"
out<-LSD.test(ag,"tech",p.adj="hommel",console=TRUE)
## 
## Study: ag ~ "tech"
## 
## LSD t Test for grain 
## P value adjustment method: hommel 
## 
## Mean Square Error:  6475313 
## 
## tech,  means and individual ( 95 %) CI
## 
##            grain      std  r      LCL      UCL     Min       Max
## CSA     3594.786 2950.114 22 2499.928 4689.644 539.844 11091.004
## Non-CSA 2674.928 2060.935 22 1580.071 3769.786 675.177  7786.251
## 
## Alpha: 0.05 ; DF Error: 42
## Critical Value of t: 2.018082 
## 
## Minimum Significant Difference: 1548.363 
## 
## Treatments with the same letter are not significantly different.
## 
##            grain groups
## CSA     3594.786      a
## Non-CSA 2674.928      a
out
## $statistics
##   MSerror Df     Mean       CV  t.value      MSD
##   6475313 42 3134.857 81.17319 2.018082 1548.363
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD    hommel   tech   2  0.05
## 
## $means
##            grain      std  r      LCL      UCL     Min       Max      Q25
## CSA     3594.786 2950.114 22 2499.928 4689.644 539.844 11091.004 1542.783
## Non-CSA 2674.928 2060.935 22 1580.071 3769.786 675.177  7786.251 1213.626
##              Q50      Q75
## CSA     2317.118 4529.520
## Non-CSA 1732.519 4191.446
## 
## $comparison
## NULL
## 
## $groups
##            grain groups
## CSA     3594.786      a
## Non-CSA 2674.928      a
## 
## attr(,"class")
## [1] "group"
theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =tech, y = grain))+ geom_boxplot(size=0.7,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey","grey","grey")) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=3,
    color="blue"
  )+ theme_classic(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=12),
          axis.text.y = element_text(face="bold", color="black", 
                           size=11))+facet_wrap(.~agroecology)+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+
  geom_text(data = yld, x = 1, y = 4400, label = "a",color="blue",size=4)+
  geom_text(data = yld, x = 2, y = 3800, label = "b",color="blue",size=4)

theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =cropping_system, y = grain))+ geom_boxplot(size=0.7,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey","grey")) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=3,
    color="blue"
  )+ theme_classic(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=12),
          axis.text.y = element_text(face="bold", color="black", 
                           size=11))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+
  geom_text(data = yld, x = 1, y = 3800, label = "b",color="blue",size=4)+
  geom_text(data = yld, x = 2, y = 3500, label = "b",color="blue",size=4)+
  geom_text(data = yld, x = 3, y = 4800, label = "a",color="blue",size=4)

ca<-subset(yld,cropping_system=="CA")

a1<-aov(grain~class,data=ca)
out<-LSD.test(a1,"class",p.adj="hommel",console=TRUE)
## 
## Study: a1 ~ "class"
## 
## LSD t Test for grain 
## P value adjustment method: hommel 
## 
## Mean Square Error:  5037591 
## 
## class,  means and individual ( 95 %) CI
## 
##         grain      std  r      LCL      UCL     Min       Max
## CA   4106.249 2617.562 32 3313.121 4899.376 624.155 10082.226
## CONV 2713.485 1795.425 32 1920.358 3506.613   0.000  6777.729
## 
## Alpha: 0.05 ; DF Error: 62
## Critical Value of t: 1.998972 
## 
## Minimum Significant Difference: 1121.652 
## 
## Treatments with the same letter are not significantly different.
## 
##         grain groups
## CA   4106.249      a
## CONV 2713.485      b
out
## $statistics
##   MSerror Df     Mean       CV  t.value      MSD
##   5037591 62 3409.867 65.82245 1.998972 1121.652
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD    hommel  class   2  0.05
## 
## $means
##         grain      std  r      LCL      UCL     Min       Max      Q25      Q50
## CA   4106.249 2617.562 32 3313.121 4899.376 624.155 10082.226 2073.808 3392.720
## CONV 2713.485 1795.425 32 1920.358 3506.613   0.000  6777.729 1439.203 2196.515
##           Q75
## CA   5472.482
## CONV 3842.888
## 
## $comparison
## NULL
## 
## $groups
##         grain groups
## CA   4106.249      a
## CONV 2713.485      b
## 
## attr(,"class")
## [1] "group"
####
biod<-lm(grain~tech+class+agroecology+residue_applied+plant_population+weeding+number_plants,data = ca)

summary(biod)
## 
## Call:
## lm(formula = grain ~ tech + class + agroecology + residue_applied + 
##     plant_population + weeding + number_plants, data = ca)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3461.4 -1134.1  -185.9   804.2  6000.7 
## 
## Coefficients: (1 not defined because of singularities)
##                           Estimate Std. Error t value Pr(>|t|)   
## (Intercept)             -995.74131 1980.45835  -0.503  0.61712   
## techNon-CSA             -658.82627 1188.67884  -0.554  0.58165   
## classCONV                       NA         NA      NA       NA   
## agroecologyLow rainfall  -55.37537  571.52509  -0.097  0.92317   
## residue_appliedYes       547.00610 1201.05082   0.455  0.65059   
## plant_population           0.08740    0.02927   2.986  0.00421 **
## weeding                  917.61604  415.54828   2.208  0.03142 * 
## number_plants           -293.43188  549.25733  -0.534  0.59533   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1960 on 55 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.3434, Adjusted R-squared:  0.2718 
## F-statistic: 4.794 on 6 and 55 DF,  p-value: 0.0005361
anova(biod)
## Analysis of Variance Table
## 
## Response: grain
##                  Df    Sum Sq  Mean Sq F value    Pr(>F)    
## tech              1  30738005 30738005  8.0044 0.0065003 ** 
## agroecology       1     51500    51500  0.0134 0.9082284    
## residue_applied   1   3647145  3647145  0.9497 0.3340504    
## plant_population  1  53565167 53565167 13.9488 0.0004478 ***
## weeding           1  21356511 21356511  5.5614 0.0219428 *  
## number_plants     1   1095992  1095992  0.2854 0.5953327    
## Residuals        55 211206536  3840119                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##################################

mbeya<-subset(yld,cropping_system=="Mbeya")

m1<-aov(grain~class,data=mbeya)
out<-LSD.test(m1,"class",p.adj="hommel",console=TRUE)
## 
## Study: m1 ~ "class"
## 
## LSD t Test for grain 
## P value adjustment method: hommel 
## 
## Mean Square Error:  8266295 
## 
## class,  means and individual ( 95 %) CI
## 
##                grain      std  r      LCL      UCL      Min       Max
## Mbeya fert  4535.499 3188.240 10 2625.358 6445.640 1339.578 11091.004
## Normal fert 3792.518 2523.434 10 1882.377 5702.659  813.145  7786.251
## 
## Alpha: 0.05 ; DF Error: 18
## Critical Value of t: 2.100922 
## 
## Minimum Significant Difference: 2701.347 
## 
## Treatments with the same letter are not significantly different.
## 
##                grain groups
## Mbeya fert  4535.499      a
## Normal fert 3792.518      a
out
## $statistics
##   MSerror Df     Mean       CV  t.value      MSD
##   8266295 18 4164.008 69.04685 2.100922 2701.347
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD    hommel  class   2  0.05
## 
## $means
##                grain      std  r      LCL      UCL      Min       Max      Q25
## Mbeya fert  4535.499 3188.240 10 2625.358 6445.640 1339.578 11091.004 1569.714
## Normal fert 3792.518 2523.434 10 1882.377 5702.659  813.145  7786.251 1502.000
##                  Q50      Q75
## Mbeya fert  4761.506 5797.551
## Normal fert 3800.357 5789.862
## 
## $comparison
## NULL
## 
## $groups
##                grain groups
## Mbeya fert  4535.499      a
## Normal fert 3792.518      a
## 
## attr(,"class")
## [1] "group"
######
biod<-lm(grain~tech+class+agroecology+residue_applied+plant_population+weeding+number_plants,data = mbeya)

summary(biod)
## 
## Call:
## lm(formula = grain ~ tech + class + agroecology + residue_applied + 
##     plant_population + weeding + number_plants, data = mbeya)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3393.1 -1293.6  -234.2   827.5  4896.5 
## 
## Coefficients: (1 not defined because of singularities)
##                           Estimate Std. Error t value Pr(>|t|)  
## (Intercept)              1.100e+04  4.803e+03   2.289   0.0394 *
## techNon-CSA             -7.127e+02  1.181e+03  -0.604   0.5565  
## classNormal fert                NA         NA      NA       NA  
## agroecologyLow rainfall  7.874e+01  1.585e+03   0.050   0.9611  
## residue_appliedYes      -1.822e+03  1.429e+03  -1.275   0.2247  
## plant_population         7.451e-02  6.007e-02   1.240   0.2368  
## weeding                 -2.004e+03  1.372e+03  -1.461   0.1677  
## number_plants           -2.698e+03  1.168e+03  -2.310   0.0379 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2582 on 13 degrees of freedom
## Multiple R-squared:  0.428,  Adjusted R-squared:  0.164 
## F-statistic: 1.621 on 6 and 13 DF,  p-value: 0.2185
anova(biod)
## Analysis of Variance Table
## 
## Response: grain
##                  Df   Sum Sq  Mean Sq F value  Pr(>F)  
## tech              1  2760105  2760105  0.4139 0.53119  
## agroecology       1  1143890  1143890  0.1715 0.68551  
## residue_applied   1    72085    72085  0.0108 0.91878  
## plant_population  1  5446893  5446893  0.8168 0.38256  
## weeding           1 19846756 19846756  2.9761 0.10818  
## number_plants     1 35591268 35591268  5.3371 0.03793 *
## Residuals        13 86692425  6668648                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
######################################
intercrop<-subset(yld,cropping_system=="Intercrop")

int<-aov(grain~class,data=intercrop)
out<-LSD.test(int,"class",p.adj="hommel",console=TRUE)
## 
## Study: int ~ "class"
## 
## LSD t Test for grain 
## P value adjustment method: hommel 
## 
## Mean Square Error:  5138851 
## 
## class,  means and individual ( 95 %) CI
## 
##                    grain      std  r      LCL      UCL     Min      Max
## Intercrop maize 3047.227 2329.012 20 2021.071 4073.382 539.844 10186.92
## Monocrop maize  2802.389 2203.044 20 1776.234 3828.544 884.808  8502.05
## 
## Alpha: 0.05 ; DF Error: 38
## Critical Value of t: 2.024394 
## 
## Minimum Significant Difference: 1451.203 
## 
## Treatments with the same letter are not significantly different.
## 
##                    grain groups
## Intercrop maize 3047.227      a
## Monocrop maize  2802.389      a
out
## $statistics
##   MSerror Df     Mean       CV  t.value      MSD
##   5138851 38 2924.808 77.50606 2.024394 1451.203
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD    hommel  class   2  0.05
## 
## $means
##                    grain      std  r      LCL      UCL     Min      Max
## Intercrop maize 3047.227 2329.012 20 2021.071 4073.382 539.844 10186.92
## Monocrop maize  2802.389 2203.044 20 1776.234 3828.544 884.808  8502.05
##                      Q25      Q50      Q75
## Intercrop maize 1726.993 2322.294 3911.369
## Monocrop maize  1629.169 1931.782 3184.842
## 
## $comparison
## NULL
## 
## $groups
##                    grain groups
## Intercrop maize 3047.227      a
## Monocrop maize  2802.389      a
## 
## attr(,"class")
## [1] "group"
######
biod<-lm(grain~tech+class+agroecology+residue_applied+plant_population+weeding+number_plants,data = intercrop)

summary(biod)
## 
## Call:
## lm(formula = grain ~ tech + class + agroecology + residue_applied + 
##     plant_population + weeding + number_plants, data = intercrop)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3176.3 -1247.0   -43.6   783.3  4743.7 
## 
## Coefficients: (1 not defined because of singularities)
##                           Estimate Std. Error t value Pr(>|t|)   
## (Intercept)              1.044e+04  3.028e+03   3.448  0.00156 **
## techNon-CSA             -3.705e+02  5.485e+02  -0.675  0.50407   
## classMonocrop maize             NA         NA      NA       NA   
## agroecologyLow rainfall -1.017e+03  5.998e+02  -1.695  0.09942 . 
## residue_appliedYes      -8.617e+01  8.326e+02  -0.103  0.91820   
## plant_population         5.548e-02  3.974e-02   1.396  0.17208   
## weeding                 -3.806e+03  1.087e+03  -3.502  0.00135 **
## number_plants           -1.099e+03  4.003e+02  -2.746  0.00970 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1703 on 33 degrees of freedom
## Multiple R-squared:  0.5116, Adjusted R-squared:  0.4228 
## F-statistic: 5.762 on 6 and 33 DF,  p-value: 0.0003457
anova(biod)
## Analysis of Variance Table
## 
## Response: grain
##                  Df   Sum Sq  Mean Sq F value   Pr(>F)   
## tech              1   599456   599456  0.2068 0.652270   
## agroecology       1 21621126 21621126  7.4586 0.010055 * 
## residue_applied   1  1604797  1604797  0.5536 0.462115   
## plant_population  1 23890176 23890176  8.2413 0.007097 **
## weeding           1 30644062 30644062 10.5712 0.002648 **
## number_plants     1 21854922 21854922  7.5392 0.009696 **
## Residuals        33 95661247  2898826                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
################################# facet wrap by cropping system
theme_set(theme_gray(base_size = 10))
ggplot(ca, aes(x =class, y = grain))+ geom_boxplot(size=0.7,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey")) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=4,
    color="blue"
  )+ theme_bw(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=12),
          axis.text.y = element_text(face="bold", color="black", 
                           size=12))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000)) +
  geom_text(data = ca, x = 1, y = 5000, label = "a",color="blue",size=5)+
  geom_text(data = ca, x = 2, y = 3500, label = "b",color="blue",size=5)

#+facet_wrap(.~agroecology)


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

theme_set(theme_gray(base_size = 10))
ggplot(mbeya, aes(x =class, y = grain))+ geom_boxplot(size=0.7,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey")) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=4,
    color="blue"
  )+ theme_bw(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=12),
          axis.text.y = element_text(face="bold", color="black", 
                           size=12))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000)) +
  geom_text(data = mbeya, x = 1, y = 3500, label = "a",color="blue",size=5)+
  geom_text(data = mbeya, x = 2, y = 5000, label = "a",color="blue",size=5)

#+facet_wrap(.~agroecology)

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

theme_set(theme_gray(base_size = 10))
ggplot(intercrop, aes(x =class, y = grain))+ geom_boxplot(size=0.7,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey")) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=4,
    color="blue"
  )+ theme_bw(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=12),
          axis.text.y = element_text(face="bold", color="black", 
                           size=12))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000)) +
  geom_text(data = intercrop, x = 1, y = 7000, label = "a",color="blue",size=6)+
  geom_text(data = intercrop, x = 2, y = 7000, label = "a",color="blue",size=6)

#+facet_wrap(.~agroecology)

########################## suggested error bar plots 


ggplot(intercrop, aes(x=class, y=grain)) + 
  stat_summary(fun.y=mean, 
               geom="bar",position=position_dodge(),colour="black",width=.5,size=.5) + 
  stat_summary(fun.ymin=min,fun.ymax=max,geom="errorbar",
               color="blue",position=position_dodge(.7), width=.1) +
  scale_fill_manual("Legend")+ 
  xlab("")+
  ylab("Grain Yield [kg/ha] ") +
  theme(panel.background = element_rect(fill = 'white', colour = 'black'))+theme_bw()+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=12),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))

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

ggplot(ca, aes(x=class, y=grain)) + 
  stat_summary(fun.y=mean, 
               geom="bar",position=position_dodge(),colour="black",width=.5,size=.5) + 
  stat_summary(fun.ymin=min,fun.ymax=max,geom="errorbar",
               color="blue",position=position_dodge(.7), width=.1) +
  scale_fill_manual("Legend")+ 
  xlab("")+
  ylab("Grain Yield [kg/ha] ") +
  theme(panel.background = element_rect(fill = 'white', colour = 'black'))+theme_bw()+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=12),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))

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

ggplot(mbeya, aes(x=class, y=grain)) + 
  stat_summary(fun.y=mean, 
               geom="bar",position=position_dodge(),colour="black",width=.5,size=.5) + 
  stat_summary(fun.ymin=min,fun.ymax=max,geom="errorbar",
               color="blue",position=position_dodge(.7), width=.1) +
  scale_fill_manual("Legend")+ 
  xlab("")+
  ylab("Grain Yield [kg/ha] ") +
  theme(panel.background = element_rect(fill = 'white', colour = 'black'))+theme_bw()+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=12),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))

######################### alll




##################
tgc <- summarySE(mbeya, measurevar="grain", groupvars=c("class","cropping_system"))
tgc
##         class cropping_system  N    grain       sd        se       ci
## 1  Mbeya fert           Mbeya 10 4535.499 3188.240 1008.2099 2280.729
## 2 Normal fert           Mbeya 10 3792.518 2523.434  797.9799 1805.156
# Error bars represent standard error of the mean
ggplot(tgc, aes(x=class, y=grain)) + 
    geom_bar(position=position_dodge(), stat="identity",width = 0.3) +
    geom_errorbar(aes(ymin=grain-se, ymax=grain+se),
                  width=.1,color="blue" ,                  # Width of the error bars
                  position=position_dodge(.9))+ 
  ylab("Grain Yield [kg/ha] ")+theme_bw()+
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=13),
          axis.text.y = element_text(face="bold",color="black", 
                           size=10))+
  geom_text(data =ca, x = 1, y = 3300, label = "a",color="blue",size=6)+
  geom_text(data = ca, x = 2, y = 2500, label = "a",color="blue",size=6)

stat_summary(geom = ‘text’, fun.y = max, position = position_dodge(.7), label = c(“a”,“a”))+

out<-LSD.test(model_1,"soil_colour",p.adj="hommel",console=TRUE)
## 
## Study: model_1 ~ "soil_colour"
## 
## LSD t Test for grain 
## P value adjustment method: hommel 
## 
## Mean Square Error:  1623275 
## 
## soil_colour,  means and individual ( 95 %) CI
## 
##          grain       std  r       LCL      UCL     Min       Max
## Black 3285.136 2469.0849 56  2932.935 3637.337 622.803 11091.004
## Brown 3888.476 2593.3389 34  3436.469 4340.483   0.000 10186.923
## Grey  2831.054 1544.7883 10  1997.594 3664.514 539.844  4591.708
## Red    473.463  213.1107  2 -1390.210 2337.136 322.771   624.155
## 
## Alpha: 0.05 ; DF Error: 23
## Critical Value of t: 2.068658 
## 
## Groups according to probability of means differences and alpha level( 0.05 )
## 
## Treatments with the same letter are not significantly different.
## 
##          grain groups
## Brown 3888.476      a
## Black 3285.136      a
## Grey  2831.054     ab
## Red    473.463      b
out
## $statistics
##   MSerror Df     Mean       CV
##   1623275 23 3386.601 37.62115
## 
## $parameters
##         test p.ajusted      name.t ntr alpha
##   Fisher-LSD    hommel soil_colour   4  0.05
## 
## $means
##          grain       std  r       LCL      UCL     Min       Max      Q25
## Black 3285.136 2469.0849 56  2932.935 3637.337 622.803 11091.004 1504.778
## Brown 3888.476 2593.3389 34  3436.469 4340.483   0.000 10186.923 1915.004
## Grey  2831.054 1544.7883 10  1997.594 3664.514 539.844  4591.708 1390.146
## Red    473.463  213.1107  2 -1390.210 2337.136 322.771   624.155  398.117
##            Q50      Q75
## Black 2302.666 4794.970
## Brown 3397.566 4700.572
## Grey  3417.270 4084.907
## Red    473.463  548.809
## 
## $comparison
## NULL
## 
## $groups
##          grain groups
## Brown 3888.476      a
## Black 3285.136      a
## Grey  2831.054     ab
## Red    473.463      b
## 
## attr(,"class")
## [1] "group"
### soil color
theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =soil_colour, y = grain))+ geom_boxplot(size=0.7,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey","grey","grey")) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=3,
    color="blue"
  )+ theme_bw(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=8))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+
  geom_text(data = yld, x = 1, y = 4000, label = "b",color="blue",size=4)+
  geom_text(data = yld, x = 2, y = 4000, label = "a",color="blue",size=4)+
  geom_text(data = yld, x = 3, y = 3800, label = "b",color="blue",size=4)+
  geom_text(data = yld, x = 4, y = 3000, label = "c",color="blue",size=4)

############## facet wrapping soil color by agroecology
theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =soil_colour, y = grain))+ geom_boxplot(size=0.7,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=3,
    color="blue"
  )+ theme_bw(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=8))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+facet_wrap(.~agroecology)

################## wrapping by soil color
theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =tech, y = grain))+ geom_boxplot(size=0.7,width = 0.5,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey","grey","grey","grey","grey","grey","grey")) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=3,
    color="blue"
  )+ theme_bw(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=8))+ 
  scale_y_continuous(breaks = c(0, 1000,3000,5000,7000,9000))+facet_wrap(.~soil_colour)+
  geom_text(data = yld, x = 1, y = 8000, label = "a",color="blue",size=4)+
  geom_text(data = yld, x = 2, y = 8000, label = "a",color="blue",size=4)

#### grain yield and texture

theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =tech, y = grain))+ geom_boxplot(size=0.7,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey","grey","grey","grey","grey")) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=3,
    color="blue"
  )+ theme_bw(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=7),
          axis.text.y = element_text(face="bold", color="black", 
                           size=8))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+
  facet_wrap(.~texture)+
  geom_text(data = yld, x = 1, y = 9000, label = "a",color="blue",size=6)+
  geom_text(data = yld, x = 2, y = 9000, label = "b",color="blue",size=6)

p <- ggplot(yld, aes(x=variet))
p + geom_bar(aes(y = (..count..)/sum(..count..)))+ theme_bw(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=8,angle = 90),
          axis.text.y = element_text(face="bold", color="black", 
                           size=8))+ 
  scale_y_continuous() +
    geom_text(aes( label = scales::percent(..count../sum(..count..)),
                   y= ..count../sum(..count..)), stat= "count", vjust = -.5,size=2.5)+ ylab("Proportion of Varieties Grown")+xlab("Varieties Grown")

#### anova model

resiapp<-aov(grain~residue_applied+agroecology,data = yld)

out<-LSD.test(resiapp,"agroecology",p.adj="hommel",console=TRUE)
## 
## Study: resiapp ~ "agroecology"
## 
## LSD t Test for grain 
## P value adjustment method: hommel 
## 
## Mean Square Error:  5588037 
## 
## agroecology,  means and individual ( 95 %) CI
## 
##                  grain      std  r      LCL      UCL     Min      Max
## High rainfall 3507.128 2322.774 80 2983.892 4030.364   0.000 10186.92
## Low rainfall  3134.857 2557.573 44 2429.325 3840.389 539.844 11091.00
## 
## Alpha: 0.05 ; DF Error: 121
## Critical Value of t: 1.979764 
## 
## Groups according to probability of means differences and alpha level( 0.05 )
## 
## Treatments with the same letter are not significantly different.
## 
##                  grain groups
## High rainfall 3507.128      a
## Low rainfall  3134.857      a
out
## $statistics
##   MSerror  Df     Mean      CV
##   5588037 121 3375.032 70.0409
## 
## $parameters
##         test p.ajusted      name.t ntr alpha
##   Fisher-LSD    hommel agroecology   2  0.05
## 
## $means
##                  grain      std  r      LCL      UCL     Min      Max      Q25
## High rainfall 3507.128 2322.774 80 2983.892 4030.364   0.000 10186.92 1732.171
## Low rainfall  3134.857 2557.573 44 2429.325 3840.389 539.844 11091.00 1332.043
##                    Q50      Q75
## High rainfall 2986.055 4741.526
## Low rainfall  2067.860 4300.314
## 
## $comparison
## NULL
## 
## $groups
##                  grain groups
## High rainfall 3507.128      a
## Low rainfall  3134.857      a
## 
## attr(,"class")
## [1] "group"
### effect of residue application

theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =residue_applied, y = grain))+ geom_boxplot(size=0.7,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey")) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("Residue Application")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=3.3,
    color="blue"
  )+ theme_bw(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=8))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+
  geom_text(data = yld, x = 1, y = 3500, label = "b",color="blue",size=5)+
  geom_text(data = yld, x = 2, y = 4500, label = "a",color="blue",size=5)

################## gender by technology
theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =sex, y = grain))+ geom_boxplot(size=0.7,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey","grey","grey")) + geom_smooth(method=lm)+ ylab("Grain Yield [kg/ha]") + xlab("")+stat_summary(fun.y=mean, geom="point", shape=10, size=2, color="blue", fill="red")+stat_summary(
    fun.data = stat_box_data , 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=3,
    color="blue"
  )+ theme_bw(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black", 
                           size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=8))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+facet_wrap(.~tech)+
  geom_text(data = yld, x = 1, y = 9000, label = "b",color="blue",size=5)+
  geom_text(data = yld, x = 2, y = 9000, label = "a",color="blue",size=5)

socio

################################## experience by sex
ggplot(yld, aes(x =experience_CA ,y =grain)) +  
  geom_point(size=3) + geom_smooth(method=lm, position = "jitter", level = 0.95, formula = y ~ x + I(x^2))+ylab("Grain Yield [kg/ha]") + 
  xlab("Experience in CA [yrs]")+theme_bw()+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+ theme( axis.line = element_line(colour = "black", 
                      size = 0.8, linetype = "solid"))+xlim(0,10)+
  stat_cor(label.x = 6)  + 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))
## Warning: Removed 20 rows containing non-finite values (stat_smooth).
## Warning: Removed 20 rows containing non-finite values (stat_cor).
## Warning: Removed 20 rows containing missing values (geom_point).

####################################age in general
ggplot(yld, aes(x =age ,y =grain)) +  
  geom_point(size=3,shape=2,color="black") + geom_smooth(method=lm, position = "jitter", level = 0.95, formula = y ~ x + I(x^2)+I(x^3))+ylab("Grain Yield [kg/ha]") + 
  xlab("Age[yrs]")+theme_bw()+ theme(axis.text.x = element_text(color="black", size=10),
          axis.text.y = element_text(color="black", 
                           size=10))+ theme( axis.line = element_line(colour = "black", 
                      size = 0.8, linetype = "solid"))+ stat_cor(method = "pearson", label.x = 55, label.y = 9000)+xlim(30,70)+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))
## Warning: Removed 8 rows containing non-finite values (stat_smooth).
## Warning: Removed 8 rows containing non-finite values (stat_cor).
## Warning: Removed 8 rows containing missing values (geom_point).

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

################################################# age by sex

ggplot(yld, aes(x =age ,y =grain,color=sex,shape=sex)) +  
  geom_point(size=3) + geom_smooth(method=lm, position = "jitter", level = 0.95, formula = y ~ x + I(x^2))+ylab("Grain Yield [kg/ha]") + 
  xlab("Age[yrs]")+theme_bw()+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+ theme( axis.line = element_line(colour = "black", 
                      size = 0.8, linetype = "solid"))+
  stat_cor(aes(color = sex), label.x = 45)  + 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))
## Warning: Removed 8 rows containing non-finite values (stat_smooth).
## Warning: Removed 8 rows containing non-finite values (stat_cor).
## Warning: Removed 8 rows containing missing values (geom_point).

plant population

## effect of plant population
ggplot(yld, aes(x =plant_population,y =grain)) +  
  geom_point(size=3) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Grain Yield [kg/ha]") + xlab("Plant population")+theme_bw()+ stat_cor(method = "pearson", label.x = 5000, label.y = 9000)+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+ 
  scale_x_continuous(breaks = c(0, 10000, 20000,30000,40000,50000,60000,70000))

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

#### facet wrap by technology
ggplot(yld, aes(x =plant_population,y =grain,color=tech,shape=tech)) +  
  geom_point(size=3) + geom_smooth(method=lm, position = "jitter", level = 0.95, formula = y ~ x + I(x^2))+ylab("Grain Yield [kg/ha]") + xlab("Plant population")+theme_bw()+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+ 
  scale_x_continuous(breaks = c(0, 10000, 20000,30000,40000,50000,60000,70000))+
  stat_cor(aes(color = tech), label.x = 5000)

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

### using the facet wrap approach
ggplot(yld, aes(x =plant_population,y =grain,color=tech,shape=tech)) +  
  geom_point(size=3) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Grain Yield [kg/ha]") + xlab("Plant population")+theme_bw()+ theme(axis.text.x = element_text(color="black", size=7),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+ 
  scale_x_continuous(breaks = c(0, 10000,30000,50000,70000))+facet_wrap(.~cropping_system)+
  stat_cor(aes(color = tech), label.x = 5000)

Effects of weeding

## effects of weeding 

ggplot(ca, aes(x =weeding,y =grain)) +  
  geom_point(size=3,shape=2) + geom_smooth(method=lm, position = "jitter", level = 0.95)+ylab("Grain Yield [kg/ha]") + xlab("Weeding")+theme_bw()+ stat_cor(method = "pearson", label.x = 2, label.y = 9000)+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))
## Warning: Removed 2 rows containing non-finite values (stat_smooth).
## Warning: Removed 2 rows containing non-finite values (stat_cor).
## Warning: Removed 2 rows containing missing values (geom_point).

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


#######################################
ggplot(ca, aes(x =weeding,y =grain,color=tech,shape=tech)) +  
  geom_point(size=3) + geom_smooth(method=lm, position = "jitter", level = 0.95, formula = y ~ x + I(x^2))+ylab("Grain Yield [kg/ha]") + xlab("Weeding")+theme_bw()+ stat_cor(method = "pearson", label.x = 1, label.y = 9000)+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+ 
  scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))+facet_wrap(.~tech) + theme(legend.position="none")
## Warning: Removed 2 rows containing non-finite values (stat_smooth).
## Warning: Removed 2 rows containing non-finite values (stat_cor).
## Warning: Removed 2 rows containing missing values (geom_point).

Number of plants

###### number of plants
## number_plants
ggplot(yld, aes(x =number_plants,y =grain)) +  
  geom_bar(stat="identity")+ylab("Grain Yield [kg/ha]") + xlab("Number of plants")+theme_bw()+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))

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

ggplot(yld, aes(x =number_plants,y =grain))+
  geom_bar(stat="identity")+ylab("Grain Yield [kg/ha]") + xlab("Number of plants")+theme_bw()+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+facet_wrap(.~tech)+ scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))

 library(readr)
rainfall <- read_csv("~/Dr Nyagumbo/Chemical analysis/final yield with socio-analysis/rainfall.csv")
## Parsed with column specification:
## cols(
##   month = col_character(),
##   label = col_character(),
##   rfall = col_double(),
##   cumm = col_double()
## )
attach(rainfall)


month=as.factor(month)
label=as.factor(label)
mo=factor(month,levels =c("Sep","Oct","Nov","Dec","Jan","Feb","Mar","Apr","May","Jun","Jul"),ordered = TRUE)
ggplot(rainfall, aes(x =mo,y =rfall,color=label,shape=label)) + geom_smooth(method=lm, position = "jitter", level = 0.95, formula = y ~ x + I(x^2))+ylab("Average Rainfall [mm]") + xlab("Month")+theme_bw()+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+ theme(legend.position = c(0.8, 0.75))

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

ggplot(data=rainfall, aes(x=mo, y=rfall, group=label,color=label,shape=label)) +
  geom_line(aes(linetype=label),size=1)+
  geom_point(aes(shape=label))+theme_bw()+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+ theme(legend.position = c(0.8, 0.75))+ylab(" Average Rainfall [mm]") + xlab("Month")

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

ggplot(data=rainfall, aes(x=mo, y=cumm, group=label,color=label,shape=label)) +
  geom_line(aes(linetype=label),size=1)+
  geom_point(aes(shape=label),size=2)+theme_bw()+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+ theme(legend.position = c(0.8, 0.45))+ 
  scale_y_continuous(breaks = c(0,200,400,600,800,1000,1200,1400,1600))+ylab("Commulative Average Rainfall [mm]") + xlab("Month")+
  geom_text(data = rainfall, x = 7, y = 1200, label = "Cyclone IDAI (March 14 and 15, 2019) ",color="blue",size=3)

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

ggplot(data=rainfall, aes(x=mo, y=cumm, group=label,color=label,shape=label))+ 
  geom_smooth(method = "loess",level = 0.95)+theme_bw()+ theme(axis.text.x = element_text(face="bold", color="black", size=10),
          axis.text.y = element_text(face="bold", color="black", 
                           size=10))+ theme(legend.position = c(0.8, 0.45))+ 
  scale_y_continuous(breaks = c(0,200,400,600,800,1000,1200,1400,1600))+ylab("Commulative Average Rainfall [mm]") + xlab("Month")