Yiel ,Socio-economic Analysis
Lovemore Chipindu [lovemore.datascience@gmail.com], [0778796212]
6/19/2019
library(ggplot2)## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlanglibrary(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 Zeileislibrary(psych)##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(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.6
## 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':
##
## diamondslibrary(agricolae)
library(easynls)
library(MVN)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2## sROC 0.1-2 loadedlibrary(lme4)## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod carlibrary(ggsignif)
library(ggpubr)## Loading required package: magrittrlibrary(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':
##
## geyserlibrary(multcompView)library(readr)
yld <- read_csv("~/Dr Nyagumbo/Chemical analysis/final yield with socio-analysis/yld.csv")## Warning: Duplicated column names deduplicated: 'biomass' =>
## 'biomass_1' [32]## 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"
## [4] "district" "village" "site"
## [7] "agroecology" "cropping_system" "age"
## [10] "sex" "hhsize" "experience_CA"
## [13] "local_soil_name" "soil_colour" "texture"
## [16] "field_treatment" "weeding" "crop_name"
## [19] "variety" "row_spacing" "plant_spacing"
## [22] "number_plants" "tillage" "residue_applied"
## [25] "residue_type" "plant_population" "planting_rains"
## [28] "seeding" "days_seeding" "biomass"
## [31] "grain" "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)
variety=as.factor(variety)
cropping_system=as.factor(cropping_system)
residue_applied=as.factor(residue_applied)
residue_type=as.factor(residue_type)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)
summary(model_1)## Df Sum Sq Mean Sq F value Pr(>F)
## tech 1 11533316 11533316 5.878 0.024930
## agroecology 1 10436099 10436099 5.319 0.031932
## site 3 35263821 11754607 5.991 0.004385
## texture 2 22237709 11118855 5.667 0.011226
## soil_colour 3 33668583 11222861 5.720 0.005390
## variety 26 128172204 4929700 2.512 0.019106
## cropping_system 2 1086174 543087 0.277 0.761076
## residue_applied 1 15355236 15355236 7.826 0.011119
## residue_type 1 1407347 1407347 0.717 0.407067
## plant_population 1 32576819 32576819 16.603 0.000591
## weeding 1 12321671 12321671 6.280 0.020966
## number_plants 1 12863370 12863370 6.556 0.018653
## age 1 3091007 3091007 1.575 0.223903
## sex 1 10099033 10099033 5.147 0.034510
## hhsize 1 1417649 1417649 0.723 0.405385
## experience_CA 1 4432519 4432519 2.259 0.148464
## tech:agroecology 1 294215 294215 0.150 0.702677
## tech:site 3 6979335 2326445 1.186 0.340349
## tech:texture 2 1142206 571103 0.291 0.750582
## tech:soil_colour 3 3454022 1151341 0.587 0.630692
## tech:cropping_system 2 1169753 584877 0.298 0.745481
## tech:residue_applied 1 2601136 2601136 1.326 0.263157
## agroecology:cropping_system 2 10023736 5011868 2.554 0.102820
## site:cropping_system 2 42429434 21214717 10.812 0.000656
## soil_colour:cropping_system 2 17553125 8776562 4.473 0.024799
## cropping_system:residue_applied 2 320241 160121 0.082 0.921941
## soil_colour:residue_applied 2 25069906 12534953 6.388 0.007155
## variety:cropping_system 2 11573929 5786965 2.949 0.075428
## tech:agroecology:cropping_system 2 162277 81138 0.041 0.959573
## tech:site:cropping_system 4 6193279 1548320 0.789 0.545844
## tech:soil_colour:cropping_system 2 2090649 1045324 0.533 0.595090
## Residuals 20 39242723 1962136
##
## tech *
## agroecology *
## site **
## texture *
## soil_colour **
## variety *
## cropping_system
## residue_applied *
## residue_type
## plant_population ***
## weeding *
## number_plants *
## age
## sex *
## hhsize
## experience_CA
## tech:agroecology
## tech:site
## tech:texture
## tech:soil_colour
## tech:cropping_system
## tech:residue_applied
## agroecology:cropping_system
## site:cropping_system ***
## soil_colour:cropping_system *
## cropping_system:residue_applied
## soil_colour:residue_applied **
## variety:cropping_system .
## tech:agroecology:cropping_system
## tech:site:cropping_system
## tech:soil_colour:cropping_system
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 24 observations deleted due to missingness#############
h<-HSD.test(model_1,"cropping_system",group =TRUE,unbalanced = TRUE)
h## $statistics
## MSerror Df Mean CV
## 1962136 20 3230.09 43.36605
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey cropping_system 3 3.577935 0.05
##
## $means
## grain std r Min Max Q25 Q50
## CA 3075.911 2038.471 56 0.000 9755.272 1504.446 2579.320
## Intercrop 3292.682 2582.729 26 539.844 10186.923 1719.497 2433.195
## Mbeya 3619.351 2506.402 18 813.145 7786.251 1389.090 3010.270
## Q75
## CA 4102.766
## Intercrop 4155.100
## Mbeya 5797.551
##
## $comparison
## NULL
##
## $groups
## grain groups
## Mbeya 3619.351 a
## Intercrop 3292.682 a
## CA 3075.911 a
##
## attr(,"class")
## [1] "group"##################################################
md_social<-aov(grain~age+sex+hhsize+experience_CA,data = yld)
summary(md_social)## Df Sum Sq Mean Sq F value Pr(>F)
## age 1 23818624 23818624 5.096 0.0262 *
## sex 1 24593605 24593605 5.262 0.0239 *
## hhsize 1 631642 631642 0.135 0.7140
## experience_CA 1 22400615 22400615 4.792 0.0309 *
## Residuals 99 462738959 4674131
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 20 observations deleted due to missingnesstechagroecology+techsite+techtexture+techsoil_colour+techcropping_system+techresidue_applied+techcropping_systemagroecology+techcropping_systemsite+techcropping_systemsoil_colour+cropping_systemsoil_colourresidue_applied+cropping_systemsoil_colourvariety
### mean and sample size function
stat_box_data <- function(grain, upper_limit = max(grain)*1.5) {
return(
data.frame(
y =9000,
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))
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))
##############################
theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =tech, y = grain))+geom_violin(width=1,size=1) +geom_boxplot(fill=c("grey","grey","grey","grey"),width=0.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=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))
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))
################################# facet wrap by cropping system
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_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=8))+facet_wrap(.~cropping_system)+
scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000)) 
###################### cropping system and agroecology
theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =tech, y = grain))+ geom_boxplot(size=0.5,varwidth = FALSE,outlier.colour = "red",outlier.shape = 2, shape=2,fill=c("grey","grey","grey","grey","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=9),
axis.text.y = element_text(face="bold", color="black",
size=8))+facet_wrap(cropping_system~agroecology)+
scale_y_continuous(breaks = c(0, 1000,3000,5000,7000,9000))
##########################################
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","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=8))+ facet_wrap(.~agroecology)+
scale_y_continuous(breaks = c(0, 1000, 2000,3000,4000,5000,6000,7000,8000,9000))
#### grain yield and texture
theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =texture, 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=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))
### 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_classic(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 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)
theme_set(theme_gray(base_size = 10))
ggplot(yld, aes(x =variety, y = grain))+ geom_boxplot(size=0.7,width = 0.5,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")+ theme_bw(base_size = 12)+ theme(axis.text.x = element_text(face="bold", color="black",
size=3),
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))
### 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,
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 technology
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","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,
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)
####################################################################socio economic and yied (yield data only)
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) + 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))
################## 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)
###
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"))+ stat_cor(method = "pearson", label.x = 7, label.y = 9000)+xlim(0,10)## 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).
################################## experience by sex
ggplot(yld, aes(x =experience_CA ,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("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)## Warning: Removed 20 rows containing non-finite values (stat_smooth).
## 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="red") + 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)## 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"))## Warning: Removed 8 rows containing non-finite values (stat_smooth).
## Warning: Removed 8 rows containing missing values (geom_point).
## 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()+ 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))
###################################################################
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()+ stat_cor(method = "pearson", label.x = 5000, label.y = 9000)+ 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))+facet_wrap(.~tech) + theme(legend.position="none")+
scale_x_continuous(breaks = c(0, 10000, 20000,30000,40000,50000,60000,70000))
### 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()+ stat_cor(method = "pearson", label.x = 5000, label.y = 9000)+ 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)
## effects of weeding
ggplot(yld, aes(x =weeding,y =grain)) +
geom_point(size=3) + 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(yld, 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_plants
ggplot(yld, aes(x =number_plants,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("Number of plants")+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))
#####################################
ggplot(yld, aes(x =number_plants,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("Number of plants")+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")
seeding
ggplot(yld, aes(x =days_seeding,y =grain)) + geom_line()+
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("Seeding date")+theme_bw()+ stat_cor(method = "pearson", label.x = 2, label.y = 9000)+ theme(axis.text.x = element_text(face="bold", color="black", size=8),
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))
Yield and soil attributes
yield_chem <- read_csv("~/Dr Nyagumbo/Chemical analysis/final yield with socio-analysis/yield_chem.csv")## Parsed with column specification:
## cols(
## tech = col_character(),
## treated = col_character(),
## bulk_density = col_double(),
## compaction = col_double(),
## om = col_double(),
## oc = col_double(),
## n = col_double(),
## pppm = col_double(),
## k = col_double(),
## ca = col_double(),
## mg = col_double(),
## pomp = col_double(),
## pomr = col_double(),
## grain = col_double()
## )View(yield_chem)
attach(yield_chem)## The following object is masked _by_ .GlobalEnv:
##
## tech## The following objects are masked from yld:
##
## grain, tech### linear model
md<-lm(grain~bulk_density+compaction+om+oc+n+pppm+k+ca+mg+pomp+pomr,data = yield_chem)
summary(md)##
## Call:
## lm(formula = grain ~ bulk_density + compaction + om + oc + n +
## pppm + k + ca + mg + pomp + pomr, data = yield_chem)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3161.9 -1531.8 -366.6 1098.0 6328.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.702e+03 5.779e+03 0.987 0.3306
## bulk_density 1.473e+02 3.883e+03 0.038 0.9700
## compaction -6.482e+02 4.255e+02 -1.524 0.1366
## om -1.382e+07 5.966e+06 -2.317 0.0265 *
## oc 2.383e+07 1.029e+07 2.317 0.0265 *
## n 2.601e+04 4.334e+04 0.600 0.5522
## pppm 1.295e+01 1.318e+01 0.983 0.3326
## k -1.012e+04 1.227e+04 -0.825 0.4152
## ca -9.681e+01 3.355e+02 -0.289 0.7746
## mg -2.465e+03 2.084e+03 -1.183 0.2449
## pomp -3.689e+02 1.804e+02 -2.045 0.0484 *
## pomr 2.713e+02 1.274e+02 2.130 0.0403 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2596 on 35 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.3413, Adjusted R-squared: 0.1342
## F-statistic: 1.648 on 11 and 35 DF, p-value: 0.1279#### poisson distribution
md<-glm(grain~bulk_density+compaction+om+oc+n+pppm+k+ca+mg+pomp+pomr,family = poisson,data = yield_chem)## Warning in dpois(y, mu, log = TRUE): non-integer x = 1839.309003## Warning in dpois(y, mu, log = TRUE): non-integer x = 1584.922213## Warning in dpois(y, mu, log = TRUE): non-integer x = 624.155031## Warning in dpois(y, mu, log = TRUE): non-integer x = 322.771429## Warning in dpois(y, mu, log = TRUE): non-integer x = 7668.866486## Warning in dpois(y, mu, log = TRUE): non-integer x = 4747.527833## Warning in dpois(y, mu, log = TRUE): non-integer x = 6164.363432## Warning in dpois(y, mu, log = TRUE): non-integer x = 7077.823389## Warning in dpois(y, mu, log = TRUE): non-integer x = 813.145432## Warning in dpois(y, mu, log = TRUE): non-integer x = 1339.577952## Warning in dpois(y, mu, log = TRUE): non-integer x = 3120.638703## Warning in dpois(y, mu, log = TRUE): non-integer x = 2028.925575## Warning in dpois(y, mu, log = TRUE): non-integer x = 923.400000## Warning in dpois(y, mu, log = TRUE): non-integer x = 771.474738## Warning in dpois(y, mu, log = TRUE): non-integer x = 1584.500342## Warning in dpois(y, mu, log = TRUE): non-integer x = 2393.380454## Warning in dpois(y, mu, log = TRUE): non-integer x = 539.844455
## Warning in dpois(y, mu, log = TRUE): non-integer x = 539.844455## Warning in dpois(y, mu, log = TRUE): non-integer x = 1200.243777
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1200.243777## Warning in dpois(y, mu, log = TRUE): non-integer x = 3410.300409## Warning in dpois(y, mu, log = TRUE): non-integer x = 4721.546365## Warning in dpois(y, mu, log = TRUE): non-integer x = 2824.164682## Warning in dpois(y, mu, log = TRUE): non-integer x = 4155.100006## Warning in dpois(y, mu, log = TRUE): non-integer x = 4300.904036## Warning in dpois(y, mu, log = TRUE): non-integer x = 6554.179060## Warning in dpois(y, mu, log = TRUE): non-integer x = 3581.388507## Warning in dpois(y, mu, log = TRUE): non-integer x = 5314.094848## Warning in dpois(y, mu, log = TRUE): non-integer x = 2485.885490## Warning in dpois(y, mu, log = TRUE): non-integer x = 2403.733009## Warning in dpois(y, mu, log = TRUE): non-integer x = 8502.050440## Warning in dpois(y, mu, log = TRUE): non-integer x = 10186.922800
## Warning in dpois(y, mu, log = TRUE): non-integer x = 10186.922800## Warning in dpois(y, mu, log = TRUE): non-integer x = 1067.821530## Warning in dpois(y, mu, log = TRUE): non-integer x = 2393.380454## Warning in dpois(y, mu, log = TRUE): non-integer x = 1584.500342## Warning in dpois(y, mu, log = TRUE): non-integer x = 2025.065422## Warning in dpois(y, mu, log = TRUE): non-integer x = 4266.872066## Warning in dpois(y, mu, log = TRUE): non-integer x = 1584.922213## Warning in dpois(y, mu, log = TRUE): non-integer x = 1839.309003## Warning in dpois(y, mu, log = TRUE): non-integer x = 2089.496170## Warning in dpois(y, mu, log = TRUE): non-integer x = 1903.711399## Warning in dpois(y, mu, log = TRUE): non-integer x = 10186.922800## Warning in dpois(y, mu, log = TRUE): non-integer x = 8502.050440## Warning in dpois(y, mu, log = TRUE): non-integer x = 1456.917402## Warning in dpois(y, mu, log = TRUE): non-integer x = 4801.466659## Warning in dpois(y, mu, log = TRUE): non-integer x = 4397.886855summary(md)##
## Call:
## glm(formula = grain ~ bulk_density + compaction + om + oc + n +
## pppm + k + ca + mg + pomp + pomr, family = poisson, data = yield_chem)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -62.22 -33.57 -6.25 18.31 86.63
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.602e+00 4.026e-02 213.641 < 2e-16 ***
## bulk_density 2.135e-01 2.754e-02 7.753 8.94e-15 ***
## compaction -2.222e-01 3.149e-03 -70.558 < 2e-16 ***
## om -2.275e+03 2.554e+01 -89.089 < 2e-16 ***
## oc 3.921e+03 4.402e+01 89.073 < 2e-16 ***
## n 5.928e+00 3.471e-01 17.082 < 2e-16 ***
## pppm 3.802e-03 8.885e-05 42.793 < 2e-16 ***
## k -2.979e+00 8.647e-02 -34.450 < 2e-16 ***
## ca -2.374e-02 2.229e-03 -10.649 < 2e-16 ***
## mg -7.430e-01 1.461e-02 -50.840 < 2e-16 ***
## pomp -1.076e-01 1.257e-03 -85.556 < 2e-16 ***
## pomr 7.714e-02 8.053e-04 95.798 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 94879 on 46 degrees of freedom
## Residual deviance: 63215 on 35 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: Inf
##
## Number of Fisher Scoring iterations: 5##### correlation
co<-data.frame(grain,bulk_density,compaction,om,oc,n,pppm,k,ca,mg,pomp,pomr)
library(PerformanceAnalytics)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo##
## Attaching package: 'PerformanceAnalytics'## The following objects are masked from 'package:agricolae':
##
## kurtosis, skewness## The following object is masked from 'package:graphics':
##
## legendchart.Correlation(co)
######
flattenCorrMatrix <- function(cormat, pmat) {
ut <- upper.tri(cormat)
data.frame(
row = rownames(cormat)[row(cormat)[ut]],
column = rownames(cormat)[col(cormat)[ut]],
cor =(cormat)[ut],
p = pmat[ut]
)
}
library(Hmisc)## Loading required package: lattice## Loading required package: Formula##
## Attaching package: 'Hmisc'## The following object is masked from 'package:psych':
##
## describe## The following objects are masked from 'package:base':
##
## format.pval, units#res2<-rcorr(as.matrix(co))
#flattenCorrMatrix(res2$r, res2$P)