library(readr)
prt_yld <- read_csv("~/CIMMYT-Zimbabwe/Dr Nyagumbo/Energy and protein/prt_yld.csv")
##
## -- Column specification --------------------------------------------------------
## cols(
## season = col_character(),
## system = col_character(),
## cum = col_double(),
## protein = col_double()
## )
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.2
library(patchwork)
## Warning: package 'patchwork' was built under R version 4.1.1
library(ggsignif)
## Warning: package 'ggsignif' was built under R version 4.1.2
library(ggpubr)
library(agricolae)
library(easynls)
library(MVN)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
## sROC 0.1-2 loaded
library(Rmisc)
## Loading required package: lattice
## Loading required package: plyr
##
## Attaching package: 'plyr'
## The following object is masked from 'package:ggpubr':
##
## mutate
library(extrafontdb)
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:patchwork':
##
## area
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(gridExtra)
library(repr) ### adjusting the length and width of your plot
library(beanplot)
library("devtools")
## Loading required package: usethis
library(SSM)
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
##
## Attaching package: 'factoextra'
## The following object is masked from 'package:agricolae':
##
## hcut
library(FactoMineR)
library(ggplot2)
library(nFactors)
##
## Attaching package: 'nFactors'
## The following object is masked from 'package:lattice':
##
## parallel
Protein
a<-ggplot(prt_yld, aes(x=season, y=protein, fill=system)) +
geom_bar(stat="identity", color="black", position=position_dodge())+
xlab("")+
ylab("Protein yield (kg/ha)")+
theme_bw()+
theme(plot.title=element_text(size=8,face="bold"),
axis.text.x=element_text(size=9,face="bold",color="black"),
axis.text.y=element_text(size=9,face="bold",color="black"))+
theme(legend.background = element_rect(fill="white", size=0.5, linetype="solid",
colour ="black"))+
theme(legend.position = c(.15, .7))
a
prt=a+plot_layout(ncol = 1)
ggsave("protein_yield.pdf",prt, width = 6, height =4, dpi=700)
###################################
ma<-ggplot(prt_yld, aes(x=season, y=protein,fill=system)) +
geom_bar(stat="identity", color="black", position=position_dodge())+
xlab("")+
ylab("Protein yield (kg/ha)")+
theme_bw()+
theme(plot.title=element_text(size=8,face="bold"),
axis.text.x=element_text(angle=45, hjust=1,size=7,face="bold",color="black"),
axis.text.y=element_text(size=9,face="bold",color="black"))+
theme(legend.background = element_rect(fill="white", size=0.5, linetype="solid",
colour ="black"))+
theme(legend.position = "none")+
facet_wrap(~system,ncol=5)
ma
zb=ma+plot_layout(ncol = 1)
ggsave("Protein revised.pdf",zb, width = 6, height =4, dpi=700)
Energy
engy<- read_csv("~/CIMMYT-Zimbabwe/Dr Nyagumbo/Energy and protein/engy.csv")
##
## -- Column specification --------------------------------------------------------
## cols(
## season = col_character(),
## system = col_character(),
## cum = col_double(),
## energy = col_double()
## )
b<-ggplot(engy, aes(x=season, y=energy, fill=system)) +
geom_bar(stat="identity", color="black", position=position_dodge())+
xlab("")+
ylab("Energy (GJ/ha)")+
theme_bw()+
theme(plot.title=element_text(size=8,face="bold"),
axis.text.x=element_text(size=9,face="bold",color="black"),
axis.text.y=element_text(size=9,face="bold",color="black"))+
theme(legend.background = element_rect(fill="white", size=0.5, linetype="solid",
colour ="black"))+
theme(legend.position = c(.2, .7))
b
en=b+plot_layout(ncol = 1)
ggsave("Energy.pdf",en, width = 6, height =4, dpi=700)
##############################################
mb<-ggplot(engy, aes(x=season, y=energy,fill=system)) +
geom_bar(stat="identity", color="black", position=position_dodge())+
xlab("")+
ylab("Energy (GJ/ha)")+
theme_bw()+
theme(plot.title=element_text(size=8,face="bold"),
axis.text.x=element_text(angle=45, hjust=1,size=7,face="bold",color="black"),
axis.text.y=element_text(size=9,face="bold",color="black"))+
theme(legend.background = element_rect(fill="white", size=0.5, linetype="solid",
colour ="black"))+
theme(legend.position = "none")+
facet_wrap(~system,ncol=5)
mb
za=mb+plot_layout(ncol = 1)
ggsave("Energy revised.pdf",za, width = 6, height =4, dpi=700)
en_prt=a+b+plot_layout(ncol = 2)
ggsave("Energy_Protein.pdf",en_prt, width = 9, height =4, dpi=700)
Rotation and inter-crops
rot_int <- read_csv("~/CIMMYT-Zimbabwe/Dr Nyagumbo/Energy and protein/rot_int.csv")
##
## -- Column specification --------------------------------------------------------
## cols(
## season = col_character(),
## cropping_system = col_character(),
## system = col_character(),
## protein = col_double(),
## energy = col_double()
## )
rot=subset(rot_int,cropping_system=="CA_Rotation")
## Rotation
rtb<-ggplot(rot, aes(x=season, y=protein,fill=system)) +
geom_bar(stat="identity", color="black", position=position_dodge())+
xlab("")+
ylab("Protein (kg/ha)")+
theme_bw()+
theme(plot.title=element_text(size=8,face="bold"),
axis.text.x=element_text(size=7,face="bold",color="black"),
axis.text.y=element_text(size=9,face="bold",color="black"))+
theme(legend.background = element_rect(fill="white", size=0.5, linetype="solid",
colour ="black"))+
theme(legend.position = c(.8, .8))
rtb
##############################
rta<-ggplot(rot, aes(x=season, y=energy,fill=system)) +
geom_bar(stat="identity", color="black", position=position_dodge())+
xlab("")+
ylab("Energy (GJ/ha)")+
theme_bw()+
theme(plot.title=element_text(size=8,face="bold"),
axis.text.x=element_text(size=7,face="bold",color="black"),
axis.text.y=element_text(size=9,face="bold",color="black"))+
theme(legend.background = element_rect(fill="white", size=0.5, linetype="solid",
colour ="black"))+
theme(legend.position = "none")
rta
## combining rta and rtb
rt=rta+rtb+plot_layout(ncol = 2)
ggsave("Rotation legumes.pdf",rt, width = 8, height =4, dpi=700)
# intercrop
int=subset(rot_int,cropping_system=="CA_Intercrop")
inta<-ggplot(int, aes(x=season, y=protein,fill=system)) +
geom_bar(stat="identity", color="black", position=position_dodge())+
xlab("")+
ylab("Protein (kg/ha)")+
theme_bw()+
theme(plot.title=element_text(size=8,face="bold"),
axis.text.x=element_text(size=7,face="bold",color="black"),
axis.text.y=element_text(size=9,face="bold",color="black"))+
theme(legend.background = element_rect(fill="white", size=0.5, linetype="solid",
colour ="black"))+
theme(legend.position = "none")
inta
##########################
intb<-ggplot(int, aes(x=season, y=energy,fill=system)) +
geom_bar(stat="identity", color="black", position=position_dodge())+
xlab("")+
ylab("Energy (GJ/ha)")+
theme_bw()+
theme(plot.title=element_text(size=8,face="bold"),
axis.text.x=element_text(size=7,face="bold",color="black"),
axis.text.y=element_text(size=9,face="bold",color="black"))+
theme(legend.background = element_rect(fill="white", size=0.5, linetype="solid",
colour ="black"))+
theme(legend.position = c(.4, .8))
intb
## combining inta and intbint
int=intb+inta+plot_layout(ncol = 2)
ggsave("Intercrop legumes.pdf",int, width = 8, height =4, dpi=700)
cow=subset(rot_int,system=="Cowpea")
## cowpea
cowa<-ggplot(cow, aes(x=season, y=protein,fill=cropping_system)) +
geom_bar(stat="identity", color="black", position=position_dodge())+
xlab("")+
ylab("Protein (kg/ha)")+
theme_bw()+
theme(plot.title=element_text(size=8,face="bold"),
axis.text.x=element_text(size=7,face="bold",color="black"),
axis.text.y=element_text(size=9,face="bold",color="black"))+
theme(legend.background = element_rect(fill="white", size=0.5, linetype="solid",
colour ="black"))+
theme(legend.position = "none")
cowa
# energy
cowb<-ggplot(cow, aes(x=season, y=energy,fill=cropping_system)) +
geom_bar(stat="identity", color="black", position=position_dodge())+
xlab("")+
ylab("Energy (GJ/ha))")+
theme_bw()+
theme(plot.title=element_text(size=8,face="bold"),
axis.text.x=element_text(size=7,face="bold",color="black"),
axis.text.y=element_text(size=9,face="bold",color="black"))+
theme(legend.background = element_rect(fill="white", size=0.5, linetype="solid",
colour ="black"))+
theme(legend.position = c(.3, .8))
cowb
## combining inta and intbint
cw=cowb+cowa+plot_layout(ncol = 2)
ggsave("cowpea in rot and int.pdf",cw, width = 8, height =4, dpi=700)
Protein and Energy full-data
pe<- read_csv("~/CIMMYT-Zimbabwe/Dr Nyagumbo/Energy and protein/energy_protein.csv")
##
## -- Column specification --------------------------------------------------------
## cols(
## season = col_character(),
## treatment = col_character(),
## cropping_system = col_character(),
## legume_species = col_character(),
## association_type = col_character(),
## rainfall = col_double(),
## maize_grain = col_double(),
## legume_grain = col_double(),
## protein = col_double(),
## energy = col_double()
## )
xa <-ggplot(data =pe,
aes(x = cropping_system,
y = protein,
fill = cropping_system,
colour =cropping_system)) +
geom_violin(alpha = 95,
position = position_nudge(x = .2, y = 0)) +
stat_ydensity(geom = "bar", colour = "white", fill = "white", width = .5,
position = position_nudge(x = -.05, y = 0)) +
geom_boxplot(width = 0.2,
outlier.shape = NA,
color = "black",
notch=FALSE)+
geom_jitter( width = .15,
size = .3,
alpha = 0.5,
color="black") +
theme_bw() +
theme(legend.position = "none") +
scale_color_brewer(palette = "Spectral") +
scale_fill_brewer(palette = "Spectral")+ ylim (0,750)+
theme(axis.text.x = element_text(face="bold", color="black",
size=9),
axis.text.y = element_text(face="bold", color="black",
size=9))+
ylab("Protein (kg/ha)")+
xlab("")
xa
## Warning: Removed 49 rows containing non-finite values (stat_ydensity).
## Warning: Removed 49 rows containing non-finite values (stat_ydensity).
## Warning: Removed 49 rows containing non-finite values (stat_boxplot).
## Warning: Removed 49 rows containing missing values (geom_point).
## Legumes associations
xx<-ggplot(data =pe,
aes(x = legume_species,
y = protein,
fill = legume_species,
colour =legume_species)) +
geom_violin(alpha = 95,
position = position_nudge(x = .2, y = 0)) +
stat_ydensity(geom = "bar", colour = "white", fill = "white", width = .5,
position = position_nudge(x = -.05, y = 0)) +
geom_boxplot(width = 0.2,
outlier.shape = NA,
color = "black",
notch=FALSE)+
geom_jitter( width = .15,
size = .3,
alpha = 0.5,
color="black") +
theme_bw() +
theme(legend.position = "none") +
scale_color_brewer(palette = "Spectral") +
scale_fill_brewer(palette = "Spectral")+ ylim (0,750)+
theme(axis.text.x = element_text(face="bold", color="black",
size=9),
axis.text.y = element_text(face="bold", color="black",
size=9))+
ylab("Protein (kg/ha)")+
xlab("")
xx
## Warning: Removed 49 rows containing non-finite values (stat_ydensity).
## Warning: Removed 49 rows containing non-finite values (stat_ydensity).
## Warning: Removed 49 rows containing non-finite values (stat_boxplot).
## Warning: Removed 49 rows containing missing values (geom_point).
###
ppx=xa+xx+plot_layout(ncol = 2)
ggsave("Systems_and legune_species.pdf",ppx, width = 10, height =5, dpi=700)
## Warning: Removed 49 rows containing non-finite values (stat_ydensity).
## Warning: Removed 49 rows containing non-finite values (stat_ydensity).
## Warning: Removed 49 rows containing non-finite values (stat_boxplot).
## Warning: Removed 49 rows containing missing values (geom_point).
## Warning: Removed 49 rows containing non-finite values (stat_ydensity).
## Warning: Removed 49 rows containing non-finite values (stat_ydensity).
## Warning: Removed 49 rows containing non-finite values (stat_boxplot).
## Warning: Removed 49 rows containing missing values (geom_point).
xb<-ggplot(data =pe,
aes(x = cropping_system,
y = energy,
fill = cropping_system,
colour =cropping_system)) +
geom_violin(alpha = 95,
position = position_nudge(x = .2, y = 0)) +
stat_ydensity(geom = "bar", colour = "white", fill = "white", width = .5,
position = position_nudge(x = -.05, y = 0)) +
geom_boxplot(width = 0.2,
outlier.shape = NA,
color = "black",
notch=FALSE) +
geom_jitter(width = .05,
size = .3,
alpha = 0.5,
color="black") +
theme_bw() +
theme(legend.position = "none") +
scale_color_brewer(palette = "Spectral") +
scale_fill_brewer(palette = "Spectral")+
theme(axis.text.x = element_text(face="bold", color="black",
size=9),
axis.text.y = element_text(face="bold", color="black",
size=9))+
ylab("Energy")+
xlab("")
xb
####### legume_association
xbb<-ggplot(data =pe,
aes(x = legume_species,
y = energy,
fill = legume_species,
colour =legume_species)) +
geom_violin(alpha = 95,
position = position_nudge(x = .2, y = 0)) +
stat_ydensity(geom = "bar", colour = "white", fill = "white", width = .5,
position = position_nudge(x = -.05, y = 0)) +
geom_boxplot(width = 0.2,
outlier.shape = NA,
color = "black",
notch=FALSE)+
geom_jitter( width = .15,
size = .3,
alpha = 0.5,
color="black") +
theme_bw() +
theme(legend.position = "none") +
scale_color_brewer(palette = "Spectral") +
scale_fill_brewer(palette = "Spectral")+ ylim (0,150)+
theme(axis.text.x = element_text(face="bold", color="black",
size=9),
axis.text.y = element_text(face="bold", color="black",
size=9))+
ylab("Energy (GJ/ha)")+
xlab("")
xbb
## Warning: Removed 1 rows containing non-finite values (stat_ydensity).
## Warning: Removed 1 rows containing non-finite values (stat_ydensity).
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1 rows containing missing values (geom_point).
####################################
## combining xa and xb
ep=xb+xbb+plot_layout(ncol = 2)
ggsave("Systems_legumes_association_energy.pdf",ep, width = 10, height =5, dpi=700)
## Warning: Removed 1 rows containing non-finite values (stat_ydensity).
## Warning: Removed 1 rows containing non-finite values (stat_ydensity).
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1 rows containing missing values (geom_point).
ya<-ggplot(data =pe,
aes(x = association_type,
y = energy,
fill = association_type,
colour =association_type)) +
geom_violin(alpha = 95,
position = position_nudge(x = .2, y = 0)) +
stat_ydensity(geom = "bar", colour = "white", fill = "white", width = .5,
position = position_nudge(x = -.05, y = 0)) +
geom_boxplot(width = 0.3,
outlier.shape = NA,
color = "black",
notch=FALSE) +
geom_jitter(width = .1,
size = .2,
alpha = 0.5,
color="black") +
theme_bw() +
theme(legend.position = "none") +
scale_color_brewer(palette = "Spectral") +
scale_fill_brewer(palette = "Spectral")+
theme(axis.text.x = element_text(angle=45,hjust=1,face="bold", color="black",
size=8),
axis.text.y = element_text(face="bold", color="black",
size=8))+
ylab("Energy (GJ/ha)")+
xlab("")+
ylim(-1,120)
ya
## Warning: Removed 16 rows containing non-finite values (stat_ydensity).
## Warning: Removed 16 rows containing non-finite values (stat_ydensity).
## Warning: Removed 16 rows containing non-finite values (stat_boxplot).
## Warning: Removed 16 rows containing missing values (geom_point).
##########
## protein
yb<-ggplot(data =pe,
aes(x = association_type,
y = protein,
fill = association_type,
colour =association_type)) +
geom_violin(alpha = 95,
position = position_nudge(x = .2, y = 0)) +
stat_ydensity(geom = "bar", colour = "white", fill = "white", width = .5,
position = position_nudge(x = -.05, y = 0)) +
geom_boxplot(width = 0.3,
outlier.shape = NA,
color = "black",
notch=FALSE) +
geom_jitter(width = .1,
size = .2,
alpha = 0.5,
color="black") +
theme_bw() +
theme(legend.position = "none") +
scale_color_brewer(palette = "Spectral") +
scale_fill_brewer(palette = "Spectral")+
theme(axis.text.x = element_text(angle=45,hjust=1,face="bold", color="black",
size=8),
axis.text.y = element_text(face="bold", color="black",
size=8))+
ylab("Protein (kg/ha)")+
xlab("")+
ylim(-1,500)
yb
## Warning: Removed 174 rows containing non-finite values (stat_ydensity).
## Warning: Removed 174 rows containing non-finite values (stat_ydensity).
## Warning: Removed 174 rows containing non-finite values (stat_boxplot).
## Warning: Removed 174 rows containing missing values (geom_point).
ey=ya+yb+plot_layout(ncol = 2)
ggsave("Energy_Protein_associations.pdf",ey, width = 8, height =4, dpi=700)
## Warning: Removed 16 rows containing non-finite values (stat_ydensity).
## Warning: Removed 16 rows containing non-finite values (stat_ydensity).
## Warning: Removed 16 rows containing non-finite values (stat_boxplot).
## Warning: Removed 16 rows containing missing values (geom_point).
## Warning: Removed 174 rows containing non-finite values (stat_ydensity).
## Warning: Removed 174 rows containing non-finite values (stat_ydensity).
## Warning: Removed 174 rows containing non-finite values (stat_boxplot).
## Warning: Removed 174 rows containing missing values (geom_point).
Association type plots
protein
Error Bars
Protein error bars
# summarySE provides the standard deviation, standard error of the mean, and a (default 95%) confidence interval
tgc <- summarySE(pe, measurevar="protein", groupvars=c("season","cropping_system"))
tgc
## season cropping_system N protein sd se ci
## 1 2011/12 CA_Intercrop 36 243.5233 77.14442 12.857404 26.10192
## 2 2011/12 CA_Rotation 48 173.3340 100.28010 14.474186 29.11831
## 3 2011/12 CA_Sole 36 258.3272 71.88131 11.980219 24.32114
## 4 2011/12 Conv_Intercrop 12 234.6267 52.59364 15.182477 33.41641
## 5 2011/12 Conv_Sole 11 253.5927 66.80594 20.142750 44.88084
## 6 2012/13 CA_Intercrop 32 230.0672 86.85258 15.353512 31.31369
## 7 2012/13 CA_Rotation 48 150.7783 48.04374 6.934517 13.95045
## 8 2012/13 CA_Sole 36 245.0628 74.27230 12.378716 25.13013
## 9 2012/13 Conv_Intercrop 12 324.6608 55.75380 16.094735 35.42427
## 10 2012/13 Conv_Sole 12 337.1267 59.71615 17.238568 37.94183
## 11 2013/14 CA_Intercrop 32 213.2319 196.52713 34.741416 70.85559
## 12 2013/14 CA_Rotation 24 687.3737 689.03812 140.649317 290.95528
## 13 2013/14 CA_Sole 36 143.9153 49.35845 8.226408 16.70050
## 14 2013/14 Conv_Intercrop 12 181.7933 38.34992 11.070670 24.36638
## 15 2013/14 Conv_Sole 12 194.2492 49.57051 14.309775 31.49560
## 16 2014/15 CA_Intercrop 54 323.3039 88.42317 12.032870 24.13488
## 17 2014/15 CA_Rotation 72 168.5799 72.57897 8.553514 17.05522
## 18 2014/15 CA_Sole 54 347.6246 107.30675 14.602599 29.28910
## 19 2014/15 Conv_Intercrop 18 360.2550 83.88167 19.771098 41.71337
## 20 2014/15 Conv_Sole 18 382.4961 150.18827 35.399715 74.68687
## 21 2015/16 CA_Intercrop 54 737.9096 382.35995 52.032598 104.36418
## 22 2015/16 CA_Rotation 36 392.4472 363.42973 60.571622 122.96693
## 23 2015/16 CA_Sole 36 389.6750 109.68070 18.280116 37.11061
## 24 2015/16 Conv_Intercrop 9 601.3589 144.84688 48.282295 111.33917
## 25 2015/16 Conv_Sole 9 421.8500 59.46024 19.820079 45.70518
## 26 2016/17 CA_Intercrop 54 400.6015 148.05342 20.147518 40.41080
## 27 2016/17 CA_Rotation 48 255.1431 99.96762 14.429084 29.02757
## 28 2016/17 CA_Sole 54 604.1102 153.50200 20.888976 41.89798
## 29 2016/17 Conv_Intercrop 18 575.1572 148.82138 35.077535 74.00713
## 30 2016/17 Conv_Sole 18 546.6222 122.07669 28.773753 60.70731
## 31 2017/18 CA_Intercrop 48 344.3502 68.12355 9.832787 19.78102
## 32 2017/18 CA_Rotation 72 141.3118 60.69310 7.152750 14.26218
## 33 2017/18 CA_Sole 54 332.5856 85.47299 11.631401 23.32964
## 34 2017/18 Conv_Intercrop 24 283.6675 89.25313 18.218719 37.68829
## 35 2017/18 Conv_Sole 18 283.5744 72.09062 16.991923 35.84982
# Standard error of the mean
ka<-ggplot(tgc, aes(x=season, y=protein, colour=cropping_system)) +
geom_errorbar(aes(ymin=protein-se, ymax=protein+se), width=.7,size=0.6) +
geom_line() +
geom_point()+
facet_wrap(~cropping_system,ncol=5)+
theme_bw()+
theme(legend.position = "none")+
xlab("")+
ylab("Protein")+
theme(axis.text.x = element_text(angle=45,hjust=1,face="bold", color="black",
size=7),
axis.text.y = element_text(face="bold", color="black",
size=8))+
geom_hline(yintercept=500, linetype='dashed', color="red",size=0.7)
ka
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
################legume_species
Energy error bars
# summarySE provides the standard deviation, standard error of the mean, and a (default 95%) confidence interval
tg <- summarySE(pe, measurevar="energy", groupvars=c("season","cropping_system"))
tg
## season cropping_system N energy sd se ci
## 1 2011/12 CA_Intercrop 36 33.24306 10.565575 1.760929 3.574876
## 2 2011/12 CA_Rotation 48 20.58250 11.969927 1.727710 3.475705
## 3 2011/12 CA_Sole 36 38.17944 10.623921 1.770653 3.594618
## 4 2011/12 Conv_Intercrop 12 34.67583 7.772054 2.243599 4.938128
## 5 2011/12 Conv_Sole 11 37.48000 9.871985 2.976515 6.632089
## 6 2012/13 CA_Intercrop 32 27.00688 8.358880 1.477655 3.013697
## 7 2012/13 CA_Rotation 48 15.99500 6.939814 1.001676 2.015112
## 8 2012/13 CA_Sole 36 36.21889 10.976304 1.829384 3.713847
## 9 2012/13 Conv_Intercrop 12 47.98250 8.240868 2.378934 5.235998
## 10 2012/13 Conv_Sole 12 49.82583 8.825780 2.547783 5.607633
## 11 2013/14 CA_Intercrop 32 19.69469 7.902606 1.396997 2.849193
## 12 2013/14 CA_Rotation 24 37.85042 29.122305 5.944566 12.297271
## 13 2013/14 CA_Sole 36 21.27000 7.295254 1.215876 2.468359
## 14 2013/14 Conv_Intercrop 12 26.86917 5.666768 1.635855 3.600493
## 15 2013/14 Conv_Sole 12 28.70917 7.326617 2.115012 4.655110
## 16 2014/15 CA_Intercrop 54 46.87259 13.383933 1.821323 3.653110
## 17 2014/15 CA_Rotation 72 20.31014 13.058588 1.538969 3.068618
## 18 2014/15 CA_Sole 54 51.37667 15.859092 2.158149 4.328699
## 19 2014/15 Conv_Intercrop 18 53.38889 12.434071 2.930739 6.183318
## 20 2014/15 Conv_Sole 18 56.53000 22.196504 5.231766 11.038062
## 21 2015/16 CA_Intercrop 54 85.91593 29.994133 4.081684 8.186822
## 22 2015/16 CA_Rotation 36 21.72556 15.062079 2.510347 5.096274
## 23 2015/16 CA_Sole 36 57.59083 16.210699 2.701783 5.484911
## 24 2015/16 Conv_Intercrop 9 91.28556 22.046201 7.348734 16.946210
## 25 2015/16 Conv_Sole 9 62.34667 8.787855 2.929285 6.754944
## 26 2016/17 CA_Intercrop 54 59.20704 21.881431 2.977686 5.972481
## 27 2016/17 CA_Rotation 48 34.81896 18.229390 2.631186 5.293263
## 28 2016/17 CA_Sole 54 89.28352 22.686513 3.087243 6.192226
## 29 2016/17 Conv_Intercrop 18 85.00500 21.994366 5.184122 10.937541
## 30 2016/17 Conv_Sole 18 80.78833 18.042403 4.252635 8.972276
## 31 2017/18 CA_Intercrop 48 47.48708 9.155624 1.321501 2.658516
## 32 2017/18 CA_Rotation 72 17.36403 11.730171 1.382414 2.756455
## 33 2017/18 CA_Sole 54 49.15481 12.632458 1.719060 3.447997
## 34 2017/18 Conv_Intercrop 24 42.05250 13.212713 2.697034 5.579240
## 35 2017/18 Conv_Sole 18 41.91056 10.654216 2.511223 5.298217
# Standard error of the mean
kb<-ggplot(tg, aes(x=season, y=energy, colour=cropping_system)) +
geom_errorbar(aes(ymin=energy-se, ymax=energy+se), width=.7,size=0.6) +
geom_line() +
geom_point()+
facet_wrap(~cropping_system,ncol=5)+
theme_bw()+
theme(legend.position = "none")+
xlab("")+
ylab("Energy")+
theme(axis.text.x = element_text(angle=45,hjust=1,face="bold", color="black",
size=7),
axis.text.y = element_text(face="bold", color="black",
size=8))+
geom_hline(yintercept=50, linetype='dashed', color="red",size=0.7)
kb
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
#######################################################
## combining the plots
ek=ka+kb+plot_layout(ncol = 1)
ggsave("Energy_Protein_Errorbars.pdf",ek, width = 8, height =5, dpi=700)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
Regression analysis
Protein
cropping systems
## protein
rg<-aov(protein~cropping_system+association_type,data=pe)
summary(rg)
## Df Sum Sq Mean Sq F value Pr(>F)
## cropping_system 4 4414219 1103555 22.34 < 2e-16 ***
## association_type 2 1891649 945825 19.15 6.56e-09 ***
## Residuals 1160 57287900 49386
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
out <- HSD.test(rg,"cropping_system", group=TRUE,console=TRUE)
##
## Study: rg ~ "cropping_system"
##
## HSD Test for protein
##
## Mean Square Error: 49386.12
##
## cropping_system, means
##
## protein std r Min Max
## CA_Intercrop 381.9975 258.4228 310 22.22 2231.84
## CA_Rotation 232.0159 265.8347 348 20.46 3162.21
## CA_Sole 348.6425 170.9581 306 47.56 1040.18
## Conv_Intercrop 361.4347 170.1207 105 57.71 814.05
## Conv_Sole 355.0117 147.0665 98 131.48 756.44
##
## Alpha: 0.05 ; DF Error: 1160
## Critical Value of Studentized Range: 3.863733
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## protein groups
## CA_Intercrop 381.9975 a
## Conv_Intercrop 361.4347 a
## Conv_Sole 355.0117 a
## CA_Sole 348.6425 a
## CA_Rotation 232.0159 b
out
## $statistics
## MSerror Df Mean CV
## 49386.12 1160 324.4106 68.50266
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey cropping_system 5 3.863733 0.05
##
## $means
## protein std r Min Max Q25 Q50 Q75
## CA_Intercrop 381.9975 258.4228 310 22.22 2231.84 235.4200 334.325 451.2550
## CA_Rotation 232.0159 265.8347 348 20.46 3162.21 116.5625 168.635 235.4150
## CA_Sole 348.6425 170.9581 306 47.56 1040.18 230.8325 318.685 431.1550
## Conv_Intercrop 361.4347 170.1207 105 57.71 814.05 238.0100 337.660 422.5400
## Conv_Sole 355.0117 147.0665 98 131.48 756.44 246.8300 323.135 438.2925
##
## $comparison
## NULL
##
## $groups
## protein groups
## CA_Intercrop 381.9975 a
## Conv_Intercrop 361.4347 a
## Conv_Sole 355.0117 a
## CA_Sole 348.6425 a
## CA_Rotation 232.0159 b
##
## attr(,"class")
## [1] "group"
# Variation range: max and min
plot(out)
#endgraph
out<-HSD.test(rg,"cropping_system", group=FALSE)
print(out$comparison)
## difference pvalue signif. LCL UCL
## CA_Intercrop - CA_Rotation 149.981625 0.0000 *** 102.56430 197.39895
## CA_Intercrop - CA_Sole 33.355065 0.3383 -15.57135 82.28148
## CA_Intercrop - Conv_Intercrop 20.562849 0.9247 -47.99278 89.11848
## CA_Intercrop - Conv_Sole 26.985781 0.8329 -43.37502 97.34659
## CA_Rotation - CA_Sole -116.626560 0.0000 *** -164.20751 -69.04561
## CA_Rotation - Conv_Intercrop -129.418776 0.0000 *** -197.02076 -61.81679
## CA_Rotation - Conv_Sole -122.995844 0.0000 *** -192.42780 -53.56389
## CA_Sole - Conv_Intercrop -12.792216 0.9865 -81.46113 55.87669
## CA_Sole - Conv_Sole -6.369284 0.9992 -76.84046 64.10189
## Conv_Intercrop - Conv_Sole 6.422932 0.9996 -78.85466 91.70053
Association type
## enery
ene <- HSD.test(rg,"association_type", group=TRUE,console=TRUE)
##
## Study: rg ~ "association_type"
##
## HSD Test for protein
##
## Mean Square Error: 49386.12
##
## association_type, means
##
## protein std r Min Max
## GntsMzRot 179.1949 88.78152 174 37.54 558.16
## MzCwInt 459.2710 354.22459 104 57.14 2231.84
## MzCwRot 284.8369 358.12032 174 20.46 3162.21
## MzPpInt 349.2144 177.89846 311 22.22 949.12
## MzSole 350.1875 165.32630 404 47.56 1040.18
##
## Alpha: 0.05 ; DF Error: 1160
## Critical Value of Studentized Range: 3.863733
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## protein groups
## MzCwInt 459.2710 a
## MzSole 350.1875 b
## MzPpInt 349.2144 b
## MzCwRot 284.8369 c
## GntsMzRot 179.1949 d
ene
## $statistics
## MSerror Df Mean CV
## 49386.12 1160 324.4106 68.50266
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey association_type 5 3.863733 0.05
##
## $means
## protein std r Min Max Q25 Q50 Q75
## GntsMzRot 179.1949 88.78152 174 37.54 558.16 114.4450 165.040 222.1625
## MzCwInt 459.2710 354.22459 104 57.14 2231.84 267.6150 371.530 462.7925
## MzCwRot 284.8369 358.12032 174 20.46 3162.21 118.6475 174.120 267.7150
## MzPpInt 349.2144 177.89846 311 22.22 949.12 228.9250 327.860 429.3600
## MzSole 350.1875 165.32630 404 47.56 1040.18 237.6275 320.545 432.8100
##
## $comparison
## NULL
##
## $groups
## protein groups
## MzCwInt 459.2710 a
## MzSole 350.1875 b
## MzPpInt 349.2144 b
## MzCwRot 284.8369 c
## GntsMzRot 179.1949 d
##
## attr(,"class")
## [1] "group"
# Variation range: max and min
plot(ene)
#endgraph
outen<-HSD.test(rg,"association_type", group=FALSE)
print(outen$comparison)
## difference pvalue signif. LCL UCL
## GntsMzRot - MzCwInt -280.0760765 0.0000 *** -355.32940 -204.822749
## GntsMzRot - MzCwRot -105.6420115 0.0001 *** -170.73511 -40.548908
## GntsMzRot - MzPpInt -170.0195522 0.0000 *** -227.49874 -112.540360
## GntsMzRot - MzSole -170.9925902 0.0000 *** -226.04714 -115.938040
## MzCwInt - MzCwRot 174.4340650 0.0000 *** 99.18074 249.687392
## MzCwInt - MzPpInt 110.0565242 0.0001 *** 41.28292 178.830129
## MzCwInt - MzSole 109.0834863 0.0001 *** 42.32306 175.843913
## MzCwRot - MzPpInt -64.3775407 0.0192 * -121.85673 -6.898349
## MzCwRot - MzSole -65.3505787 0.0107 * -120.40513 -10.296029
## MzPpInt - MzSole -0.9730379 1.0000 -46.77423 44.828150
Energy
# energy
rge<-aov(energy~cropping_system+association_type,data=pe)
summary(rge)
## Df Sum Sq Mean Sq F value Pr(>F)
## cropping_system 4 202077 50519 94.004 <2e-16 ***
## association_type 2 1886 943 1.755 0.173
## Residuals 1160 623407 537
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
outrge <- HSD.test(rge,"cropping_system", group=TRUE,console=TRUE)
##
## Study: rge ~ "cropping_system"
##
## HSD Test for energy
##
## Mean Square Error: 537.4195
##
## cropping_system, means
##
## energy std r Min Max
## CA_Intercrop 49.47852 27.03815 310 3.28 144.86
## CA_Rotation 22.50029 16.22546 348 0.97 149.42
## CA_Sole 51.52725 25.26640 306 7.03 153.73
## Conv_Intercrop 53.67857 25.47739 105 8.64 120.31
## Conv_Sole 52.46878 21.73550 98 19.43 111.80
##
## Alpha: 0.05 ; DF Error: 1160
## Critical Value of Studentized Range: 3.863733
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## energy groups
## Conv_Intercrop 53.67857 a
## Conv_Sole 52.46878 a
## CA_Sole 51.52725 a
## CA_Intercrop 49.47852 a
## CA_Rotation 22.50029 b
outrge
## $statistics
## MSerror Df Mean CV
## 537.4195 1160 42.5998 54.41882
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey cropping_system 5 3.863733 0.05
##
## $means
## energy std r Min Max Q25 Q50 Q75
## CA_Intercrop 49.47852 27.03815 310 3.28 144.86 29.7175 45.665 61.4500
## CA_Rotation 22.50029 16.22546 348 0.97 149.42 9.8975 20.115 30.7850
## CA_Sole 51.52725 25.26640 306 7.03 153.73 34.1150 47.100 63.7250
## Conv_Intercrop 53.67857 25.47739 105 8.64 120.31 35.2200 50.010 62.4500
## Conv_Sole 52.46878 21.73550 98 19.43 111.80 36.4850 47.760 64.7775
##
## $comparison
## NULL
##
## $groups
## energy groups
## Conv_Intercrop 53.67857 a
## Conv_Sole 52.46878 a
## CA_Sole 51.52725 a
## CA_Intercrop 49.47852 a
## CA_Rotation 22.50029 b
##
## attr(,"class")
## [1] "group"
# Variation range: max and min
plot(outrge)
############################
g <- HSD.test(rge,"association_type", group=TRUE,console=TRUE)
##
## Study: rge ~ "association_type"
##
## HSD Test for energy
##
## Mean Square Error: 537.4195
##
## association_type, means
##
## energy std r Min Max
## GntsMzRot 22.57063 12.96847 174 3.84 82.49
## MzCwInt 46.10865 25.57063 104 7.84 133.90
## MzCwRot 22.42994 18.97002 174 0.97 149.42
## MzPpInt 52.02344 26.92366 311 3.28 144.86
## MzSole 51.75564 24.43409 404 7.03 153.73
##
## Alpha: 0.05 ; DF Error: 1160
## Critical Value of Studentized Range: 3.863733
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## energy groups
## MzPpInt 52.02344 a
## MzSole 51.75564 a
## MzCwInt 46.10865 a
## GntsMzRot 22.57063 b
## MzCwRot 22.42994 b
g
## $statistics
## MSerror Df Mean CV
## 537.4195 1160 42.5998 54.41882
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey association_type 5 3.863733 0.05
##
## $means
## energy std r Min Max Q25 Q50 Q75
## GntsMzRot 22.57063 12.96847 174 3.84 82.49 12.2050 20.070 30.4150
## MzCwInt 46.10865 25.57063 104 7.84 133.90 28.1825 42.530 56.9725
## MzCwRot 22.42994 18.97002 174 0.97 149.42 6.5575 20.455 31.4475
## MzPpInt 52.02344 26.92366 311 3.28 144.86 33.8950 48.730 64.0350
## MzSole 51.75564 24.43409 404 7.03 153.73 35.1175 47.375 63.9675
##
## $comparison
## NULL
##
## $groups
## energy groups
## MzPpInt 52.02344 a
## MzSole 51.75564 a
## MzCwInt 46.10865 a
## GntsMzRot 22.57063 b
## MzCwRot 22.42994 b
##
## attr(,"class")
## [1] "group"
# Variation range: max and min
plot(g)
#endgraph
outg<-HSD.test(rge,"association_type", group=FALSE)
print(outg$comparison)
## difference pvalue signif. LCL UCL
## GntsMzRot - MzCwInt -23.5380217 0.0000 *** -31.388208 -15.687836
## GntsMzRot - MzCwRot 0.1406897 1.0000 -6.649614 6.930994
## GntsMzRot - MzPpInt -29.4528083 0.0000 *** -35.448853 -23.456763
## GntsMzRot - MzSole -29.1850114 0.0000 *** -34.928126 -23.441897
## MzCwInt - MzCwRot 23.6787113 0.0000 *** 15.828525 31.528897
## MzCwInt - MzPpInt -5.9147867 0.1614 -13.089029 1.259455
## MzCwInt - MzSole -5.6469897 0.1747 -12.611224 1.317244
## MzCwRot - MzPpInt -29.5934980 0.0000 *** -35.589543 -23.597453
## MzCwRot - MzSole -29.3257010 0.0000 *** -35.068815 -23.582587
## MzPpInt - MzSole 0.2677970 0.9999 -4.510036 5.045630
Relationship assessment
my.formula <- energy~ protein
ra<- ggplot(pe, aes(x =protein, y = energy, color=legume_species)) +
geom_point() +
geom_smooth(method=lm, position = "jitter", level = 0.95)+
ylab("Energy (GJ/ha") +
xlab("Protein (kg/ha)")+
theme(legend.position = c(0.2, 0.85))+
theme_bw()+
stat_regline_equation(label.x = 5, label.y = 2.5,size=3)+
stat_cor(method="pearson",label.x=5,label.y=2.7,size=3)+
theme(axis.text.x = element_text(face="bold", color="black",
size=8),
axis.text.y = element_text(face="bold", color="black",
size=8))+
theme(legend.position='none')+
facet_wrap(~cropping_system)
ra
## `geom_smooth()` using formula 'y ~ x'
Linear Models
Protein Linear Model
plm<-lm(protein~cropping_system+association_type+season,data=pe)
summary(plm)
##
## Call:
## lm(formula = protein ~ cropping_system + association_type + season,
## data = pe)
##
## Residuals:
## Min 1Q Median 3Q Max
## -505.56 -102.17 -13.33 65.23 2913.57
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 229.93 20.81 11.048 < 2e-16 ***
## cropping_systemCA_Rotation -139.04 20.17 -6.894 8.90e-12 ***
## cropping_systemCA_Sole 18.35 17.56 1.045 0.29613
## cropping_systemConv_Intercrop 36.08 23.23 1.553 0.12074
## cropping_systemConv_Sole 29.96 24.03 1.247 0.21278
## association_typeMzCwInt 114.20 23.28 4.906 1.06e-06 ***
## association_typeMzCwRot 105.64 21.08 5.012 6.24e-07 ***
## association_typeMzPpInt NA NA NA NA
## association_typeMzSole NA NA NA NA
## season2012/13 6.75 23.40 0.288 0.77309
## season2013/14 52.11 24.65 2.114 0.03476 *
## season2014/15 65.02 21.21 3.066 0.00222 **
## season2015/16 304.83 23.28 13.092 < 2e-16 ***
## season2016/17 221.09 21.75 10.163 < 2e-16 ***
## season2017/18 40.12 21.21 1.892 0.05879 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 196.6 on 1154 degrees of freedom
## Multiple R-squared: 0.2985, Adjusted R-squared: 0.2912
## F-statistic: 40.93 on 12 and 1154 DF, p-value: < 2.2e-16
####### Anova Models With Interractions
preg<-aov(protein~cropping_system+association_type+season+cropping_system*association_type*season,data=pe)
summary(preg)
## Df Sum Sq Mean Sq F value Pr(>F)
## cropping_system 4 4414219 1103555 61.179 <2e-16 ***
## association_type 2 1891649 945825 52.435 <2e-16 ***
## season 6 12679019 2113170 117.151 <2e-16 ***
## cropping_system:association_type 1 77980 77980 4.323 0.0378 *
## cropping_system:season 24 13075258 544802 30.203 <2e-16 ***
## association_type:season 12 11307169 942264 52.238 <2e-16 ***
## Residuals 1117 20148474 18038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Energy Linear Model
elm<-lm(energy~cropping_system+association_type+season,data=pe)
summary(elm)
##
## Call:
## lm(formula = energy ~ cropping_system + association_type + season,
## data = pe)
##
## Residuals:
## Min 1Q Median 3Q Max
## -63.643 -11.704 -1.336 9.417 145.850
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.0916 1.9583 19.962 < 2e-16 ***
## cropping_systemCA_Rotation -26.3784 1.8978 -13.899 < 2e-16 ***
## cropping_systemCA_Sole 1.7768 1.6522 1.075 0.28242
## cropping_systemConv_Intercrop 4.9486 2.1863 2.263 0.02379 *
## cropping_systemConv_Sole 3.1689 2.2612 1.401 0.16136
## association_typeMzCwInt -5.7001 2.1903 -2.602 0.00938 **
## association_typeMzCwRot -0.1407 1.9835 -0.071 0.94347
## association_typeMzPpInt NA NA NA NA
## association_typeMzSole NA NA NA NA
## season2012/13 -1.4964 2.2023 -0.679 0.49698
## season2013/14 -9.0028 2.3199 -3.881 0.00011 ***
## season2014/15 9.5713 1.9955 4.796 1.83e-06 ***
## season2015/16 28.9312 2.1910 13.205 < 2e-16 ***
## season2016/17 32.8051 2.0471 16.025 < 2e-16 ***
## season2017/18 5.7360 1.9959 2.874 0.00413 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.5 on 1154 degrees of freedom
## Multiple R-squared: 0.5226, Adjusted R-squared: 0.5176
## F-statistic: 105.3 on 12 and 1154 DF, p-value: < 2.2e-16
####### Anova Models With Interractions
ereg<-aov(energy~cropping_system+association_type+season+cropping_system*association_type*season,data=pe)
summary(ereg)
## Df Sum Sq Mean Sq F value Pr(>F)
## cropping_system 4 202077 50519 223.144 < 2e-16 ***
## association_type 2 1886 943 4.165 0.0158 *
## season 6 228415 38069 168.152 < 2e-16 ***
## cropping_system:association_type 1 11 11 0.051 0.8220
## cropping_system:season 24 121687 5070 22.396 < 2e-16 ***
## association_type:season 12 20407 1701 7.511 1.82e-13 ***
## Residuals 1117 252886 226
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yield ANOVA Model
# yield
mzge<-aov(maize_grain~cropping_system+association_type,data=pe)
summary(mzge)
## Df Sum Sq Mean Sq F value Pr(>F)
## cropping_system 4 2.023e+08 50571093 21.482 < 2e-16 ***
## association_type 2 2.540e+07 12701069 5.395 0.00468 **
## Residuals 968 2.279e+09 2354088
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 192 observations deleted due to missingness
outmz <- HSD.test(mzge,"cropping_system", group=TRUE,console=TRUE)
##
## Study: mzge ~ "cropping_system"
##
## HSD Test for maize_grain
##
## Mean Square Error: 2354088
##
## cropping_system, means
##
## maize_grain std r Min Max
## CA_Intercrop 2748.667 1353.922 310 46.85 7013.89
## CA_Rotation 4069.375 1607.868 156 1205.83 11163.30
## CA_Sole 3486.424 1709.582 306 475.64 10401.79
## Conv_Intercrop 3352.442 1505.777 105 467.17 8140.50
## Conv_Sole 3550.112 1470.667 98 1314.76 7564.40
##
## Alpha: 0.05 ; DF Error: 968
## Critical Value of Studentized Range: 3.86494
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## maize_grain groups
## CA_Rotation 4069.375 a
## Conv_Sole 3550.112 ab
## CA_Sole 3486.424 b
## Conv_Intercrop 3352.442 b
## CA_Intercrop 2748.667 c
outmz
## $statistics
## MSerror Df Mean CV
## 2354088 968 3337.1 45.97717
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey cropping_system 5 3.86494 0.05
##
## $means
## maize_grain std r Min Max Q25 Q50
## CA_Intercrop 2748.667 1353.922 310 46.85 7013.89 1763.615 2756.200
## CA_Rotation 4069.375 1607.868 156 1205.83 11163.30 2947.648 3683.315
## CA_Sole 3486.424 1709.582 306 475.64 10401.79 2308.300 3186.850
## Conv_Intercrop 3352.442 1505.777 105 467.17 8140.50 2272.330 3195.580
## Conv_Sole 3550.112 1470.667 98 1314.76 7564.40 2468.323 3231.345
## Q75
## CA_Intercrop 3516.113
## CA_Rotation 4875.382
## CA_Sole 4311.558
## Conv_Intercrop 3927.670
## Conv_Sole 4382.905
##
## $comparison
## NULL
##
## $groups
## maize_grain groups
## CA_Rotation 4069.375 a
## Conv_Sole 3550.112 ab
## CA_Sole 3486.424 b
## Conv_Intercrop 3352.442 b
## CA_Intercrop 2748.667 c
##
## attr(,"class")
## [1] "group"
# Variation range: max and min
print(mzge$comparison)
## NULL
############################
gmz <- HSD.test(mzge,"association_type", group=TRUE,console=TRUE)
##
## Study: mzge ~ "association_type"
##
## HSD Test for maize_grain
##
## Mean Square Error: 2354088
##
## association_type, means
##
## maize_grain std r Min Max
## GntsMzRot 4157.833 1741.391 78 1205.83 11163.30
## MzCwInt 2371.822 1407.565 104 46.85 6694.34
## MzCwRot 3980.917 1468.250 78 1783.61 8476.23
## MzPpInt 3078.533 1377.041 311 222.19 8140.50
## MzSole 3501.873 1653.264 404 475.64 10401.79
##
## Alpha: 0.05 ; DF Error: 968
## Critical Value of Studentized Range: 3.86494
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## maize_grain groups
## GntsMzRot 4157.833 a
## MzCwRot 3980.917 ab
## MzSole 3501.873 b
## MzPpInt 3078.533 c
## MzCwInt 2371.822 d
gmz
## $statistics
## MSerror Df Mean CV
## 2354088 968 3337.1 45.97717
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey association_type 5 3.86494 0.05
##
## $means
## maize_grain std r Min Max Q25 Q50 Q75
## GntsMzRot 4157.833 1741.391 78 1205.83 11163.30 2894.615 3791.525 5084.260
## MzCwInt 2371.822 1407.565 104 46.85 6694.34 1266.757 2378.455 3358.970
## MzCwRot 3980.917 1468.250 78 1783.61 8476.23 3028.490 3623.015 4624.150
## MzPpInt 3078.533 1377.041 311 222.19 8140.50 2132.030 2990.190 3755.120
## MzSole 3501.873 1653.264 404 475.64 10401.79 2376.295 3205.450 4328.073
##
## $comparison
## NULL
##
## $groups
## maize_grain groups
## GntsMzRot 4157.833 a
## MzCwRot 3980.917 ab
## MzSole 3501.873 b
## MzPpInt 3078.533 c
## MzCwInt 2371.822 d
##
## attr(,"class")
## [1] "group"
# Variation range: max and min
plot(gmz)
#endgraph
outgmz<-HSD.test(rg,"association_type", group=FALSE)
print(outgmz$comparison)
## difference pvalue signif. LCL UCL
## GntsMzRot - MzCwInt -280.0760765 0.0000 *** -355.32940 -204.822749
## GntsMzRot - MzCwRot -105.6420115 0.0001 *** -170.73511 -40.548908
## GntsMzRot - MzPpInt -170.0195522 0.0000 *** -227.49874 -112.540360
## GntsMzRot - MzSole -170.9925902 0.0000 *** -226.04714 -115.938040
## MzCwInt - MzCwRot 174.4340650 0.0000 *** 99.18074 249.687392
## MzCwInt - MzPpInt 110.0565242 0.0001 *** 41.28292 178.830129
## MzCwInt - MzSole 109.0834863 0.0001 *** 42.32306 175.843913
## MzCwRot - MzPpInt -64.3775407 0.0192 * -121.85673 -6.898349
## MzCwRot - MzSole -65.3505787 0.0107 * -120.40513 -10.296029
## MzPpInt - MzSole -0.9730379 1.0000 -46.77423 44.828150
Principal Component Analysis
library(readr)
descrip <- read_csv("ep_pc.csv")
##
## -- Column specification --------------------------------------------------------
## cols(
## season = col_character(),
## cropping_system = col_character(),
## legume_species = col_character(),
## rainfall_ca = col_character(),
## rainfall = col_double(),
## maize_grain = col_double()
## )
res.pca <- prcomp(descrip[, -1:-4], scale = TRUE)
fviz_eig(res.pca)
fviz_pca_var(res.pca,
col.var = "contrib", # Color by contributions to the PC
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Avoid text overlappin
)
m<-fviz_pca_ind(res.pca, label="none", habillage=descrip$cropping_system)+
theme_bw()+
theme(legend.position = "none")
m
p <- fviz_pca_ind(res.pca, label="none", habillage=descrip$cropping_system,
addEllipses=TRUE, ellipse.level=0.95)+
theme_bw()
print(p)
mp=m+p+plot_layout(ncol = 2)
ggsave("PCA plots.pdf",mp, width =9, height =4, dpi=700)
Mean Differences Plots
library(readr)
mn_df <- read_csv("~/CIMMYT-Zimbabwe/Dr Nyagumbo/Energy and protein/mn_df.csv")
##
## -- Column specification --------------------------------------------------------
## cols(
## model = col_character(),
## association = col_character(),
## difference = col_double()
## )
a4 = ggplot(mn_df,
aes(x = reorder(association,difference),y = difference, ymin=min(difference), ymax = max (difference)))+ geom_pointrange(aes(col=association))+
xlab('')+ ylab("Association type (Mean)-Maize sole (Mean)")+
geom_errorbar(aes(ymin=min(difference), ymax=max(difference),col=association),width=0.5,cex=1)+
facet_wrap(~model)+
theme_bw()+
theme(plot.title=element_text(size=10,face="bold"),
axis.text.x=element_text(size=10,face="bold",color="black"),
axis.text.y=element_text(size=10,color="black"),
axis.title=element_text(size=10,face="bold"),
strip.text.y = element_text(hjust=0,vjust = 0,face="bold"))+ coord_flip()+theme(legend.position = "")+
theme(plot.title = element_text(hjust = 0.5,color="black", size=12, face="bold"))+
geom_hline(aes(yintercept = 0),linetype="dashed",color = "red", size=0.5) +
geom_text(aes(label =difference), vjust = -0.5, size = 3)
a4
ba=a4+plot_layout(ncol = 1)
ggsave("Mean_Differences.pdf",ba, width =9, height =2.5, dpi=700)