Loading the data set

library(readr)
wheat<- read_csv("~/Dr Mupangwa/mechanization paper/Complete Mechinization Data Analysis/mech_wheat.csv")

## Parsed with column specification:
## cols(
##   year = col_double(),
##   region = col_character(),
##   district = col_character(),
##   agroecology = col_character(),
##   soil_type = col_character(),
##   treatment = col_character(),
##   grain_yield = col_double(),
##   stover = col_double(),
##   hi = col_double(),
##   gross_margins = col_double()
## )

attach(wheat)
names(wheat)

##  [1] "year"          "region"        "district"      "agroecology"  
##  [5] "soil_type"     "treatment"     "grain_yield"   "stover"       
##  [9] "hi"            "gross_margins"

year=as.factor(year)
treatment=as.factor(treatment)
agroecology=as.factor(agroecology)
soil_type=as.factor(soil_type)

Required Libraries

library(agricolae)
library(ggplot2)
library(Hmisc)

## Loading required package: lattice

## Loading required package: survival

## Loading required package: Formula

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:base':
## 
##     format.pval, units

library(pastecs)
library(psych)

## 
## Attaching package: 'psych'

## The following object is masked from 'package:Hmisc':
## 
##     describe

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(doBy)
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

## 
## Attaching package: 'xts'

## The following objects are masked from 'package:pastecs':
## 
##     first, last

## 
## Attaching package: 'PerformanceAnalytics'

## The following objects are masked from 'package:agricolae':
## 
##     kurtosis, skewness

## The following object is masked from 'package:graphics':
## 
##     legend

library(psych)
library(gridExtra)
library(reshape2)
library(ggpubr)

## Loading required package: magrittr

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:pastecs':
## 
##     extract

library(ggpmisc)
library(patchwork) ## new package for displaying plots

Descriptive Statistics

describe.by(grain_yield,year)

## Warning: describe.by is deprecated. Please use the describeBy function

## 
##  Descriptive statistics by group 
## group: 2018
##    vars  n    mean      sd median trimmed     mad  min  max range skew kurtosis
## X1    1 22 3455.77 1277.43   3351    3268 1025.22 1939 6699  4760 1.16     0.77
##        se
## X1 272.35
## ------------------------------------------------------------ 
## group: 2019
##    vars  n mean      sd median trimmed     mad min  max range skew kurtosis
## X1    1 30 3223 1635.06   3096 3087.42 1927.38 752 7397  6645 0.65    -0.31
##        se
## X1 298.52

describe.by(grain_yield,treatment)

## Warning: describe.by is deprecated. Please use the describeBy function

## 
##  Descriptive statistics by group 
## group: CP
##    vars  n    mean      sd median trimmed     mad min  max range skew kurtosis
## X1    1 26 3040.92 1379.88 2883.5 2959.27 1314.32 752 6523  5771 0.66    -0.09
##        se
## X1 270.62
## ------------------------------------------------------------ 
## group: No-till
##    vars  n    mean      sd median trimmed    mad  min  max range skew kurtosis
## X1    1 26 3602.04 1559.37 3351.5 3482.32 1504.1 1313 7397  6084 0.74    -0.11
##        se
## X1 305.82

describe.by(grain_yield,agroecology)

## Warning: describe.by is deprecated. Please use the describeBy function

## 
##  Descriptive statistics by group 
## group: Semi-arid
##    vars n    mean     sd median trimmed    mad  min  max range skew kurtosis
## X1    1 8 2834.38 697.69   2817 2834.38 851.01 1939 3850  1911 0.08    -1.71
##        se
## X1 246.67
## ------------------------------------------------------------ 
## group: Sub-humid
##    vars  n    mean      sd median trimmed     mad min  max range skew kurtosis
## X1    1 44 3410.05 1576.24   3196 3281.17 1638.27 752 7397  6645 0.62    -0.23
##        se
## X1 237.63

describe.by(grain_yield,soil_type)

## Warning: describe.by is deprecated. Please use the describeBy function

## 
##  Descriptive statistics by group 
## group: Clay loams
##    vars  n    mean      sd median trimmed     mad min  max range skew kurtosis
## X1    1 32 2937.56 1404.97 2688.5 2793.19 1382.52 752 6699  5947 0.93     0.57
##        se
## X1 248.37
## ------------------------------------------------------------ 
## group: Loamy sands
##    vars n    mean     sd median trimmed    mad  min  max range skew kurtosis
## X1    1 8 2834.38 697.69   2817 2834.38 851.01 1939 3850  1911 0.08    -1.71
##        se
## X1 246.67
## ------------------------------------------------------------ 
## group: Vertisol
##    vars  n mean      sd median trimmed     mad  min  max range skew kurtosis
## X1    1 12 4670 1329.71   4396  4547.2 1584.16 3171 7397  4226 0.58    -0.84
##        se
## X1 383.85

#Grain Yield model

md<-lm(grain_yield~treatment+year+agroecology+soil_type,data = wheat)
summary(md)

## 
## Call:
## lm(formula = grain_yield ~ treatment + year + agroecology + soil_type, 
##     data = wheat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1889.4  -773.2  -159.2   494.2  2992.7 
## 
## Coefficients: (1 not defined because of singularities)
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          3135546.9   819070.0   3.828 0.000381 ***
## treatmentNo-till         561.1      315.9   1.776 0.082169 .  
## year                   -1552.5      405.9  -3.825 0.000384 ***
## agroecologySub-humid     976.5      504.8   1.934 0.059102 .  
## soil_typeLoamy sands        NA         NA      NA       NA    
## soil_typeVertisol       2411.7      424.5   5.681 8.14e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1139 on 47 degrees of freedom
## Multiple R-squared:  0.4579, Adjusted R-squared:  0.4118 
## F-statistic: 9.926 on 4 and 47 DF,  p-value: 6.617e-06

anova_md<-aov(grain_yield~treatment+year+agroecology+soil_type+treatment*agroecology+treatment*soil_type,data = wheat)
summary(anova_md)

##                       Df   Sum Sq  Mean Sq F value   Pr(>F)    
## treatment              1  4093056  4093056   3.047   0.0877 .  
## year                   1   687709   687709   0.512   0.4780    
## agroecology            1  4854270  4854270   3.614   0.0637 .  
## soil_type              1 41876180 41876180  31.175 1.29e-06 ***
## treatment:agroecology  1   339274   339274   0.253   0.6177    
## treatment:soil_type    1   188908   188908   0.141   0.7094    
## Residuals             45 60446289  1343251                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#####################
h<-HSD.test(anova_md,"treatment",group =TRUE)
h

## $statistics
##   MSerror Df     Mean       CV      MSD
##   1343251 45 3321.481 34.89368 647.4238
## 
## $parameters
##    test    name.t ntr StudentizedRange alpha
##   Tukey treatment   2         2.848372  0.05
## 
## $means
##         grain_yield      std  r  Min  Max  Q25    Q50     Q75
## CP         3040.923 1379.885 26  752 6523 1991 2883.5 3708.25
## No-till    3602.038 1559.366 26 1313 7397 2611 3351.5 4376.00
## 
## $comparison
## NULL
## 
## $groups
##         grain_yield groups
## No-till    3602.038      a
## CP         3040.923      a
## 
## attr(,"class")
## [1] "group"

a<-HSD.test(anova_md,"soil_type",group =TRUE)
a

## $statistics
##   MSerror Df     Mean       CV
##   1343251 45 3321.481 34.89368
## 
## $parameters
##    test    name.t ntr StudentizedRange alpha
##   Tukey soil_type   3         3.427507  0.05
## 
## $means
##             grain_yield      std  r  Min  Max    Q25    Q50     Q75
## Clay loams     2937.562 1404.969 32  752 6699 1832.5 2688.5 3656.25
## Loamy sands    2834.375  697.686  8 1939 3850 2331.5 2817.0 3304.00
## Vertisol       4670.000 1329.712 12 3171 7397 3660.5 4396.0 5461.75
## 
## $comparison
## NULL
## 
## $groups
##             grain_yield groups
## Vertisol       4670.000      a
## Clay loams     2937.562      b
## Loamy sands    2834.375      b
## 
## attr(,"class")
## [1] "group"

Gross Margins Anova Table

gm<-aov(gross_margins~treatment+year+agroecology+soil_type+treatment*agroecology+treatment*soil_type,data = wheat)
summary(gm)

##                       Df   Sum Sq Mean Sq F value   Pr(>F)    
## treatment              1  1672829 1672829   5.826   0.0199 *  
## year                   1    49582   49582   0.173   0.6797    
## agroecology            1    12277   12277   0.043   0.8371    
## soil_type              1  9650280 9650280  33.606 6.24e-07 ***
## treatment:agroecology  1    79516   79516   0.277   0.6013    
## treatment:soil_type    1    37814   37814   0.132   0.7184    
## Residuals             45 12921996  287155                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Violin plots

### mean and sample size function 
stat_box_data <- function(grain_yield, upper_limit = max(grain_yield)*1.5) {
  return( 
    data.frame(
      y = 8000,
      label = paste(
                    'Mean =', round(mean(grain_yield), 0), '\n')
    )
  )
}

Treatment by grain yield

theme_set(theme_grey(base_size = 10)) ###This sets the font sizes of anything writt
ab<-ggplot(wheat, aes(x = treatment, y = grain_yield))+ 
geom_violin(size=0.9,shape=8)+
  geom_boxplot(size=0.5,width = 0.4,shape=6,fill=c("grey","grey")) +
geom_smooth(method=lm)+ 
ylab("Wheat Grain Yield [kg/ha]") + 
xlab("")+
stat_summary(fun.y=mean, geom="point", shape=10, size=3, color="blue", fill="red")+
stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=4
  )+ 
  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=10))+ 
  theme( axis.line = element_line(colour = "black", 
                      size = 0.8, linetype = "solid"))

## Warning: Ignoring unknown parameters: shape

ab

Grain yield by season

theme_set(theme_grey(base_size = 10)) ###This sets the font sizes of anything writt
c<-ggplot(wheat, aes(x = treatment, y = grain_yield))+ 
geom_violin(size=0.9,shape=8)+
  geom_boxplot(size=0.5,width = 0.2,shape=6,fill=c("grey","grey","grey","grey")) +
geom_smooth(method=lm)+ 
ylab("Wheat Grain Yield [kg/ha]") + 
xlab("")+
stat_summary(fun.y=mean, geom="point", shape=10, size=3, color="blue", fill="red")+
stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=4
  )+ 
  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=10))+ 
  theme( axis.line = element_line(colour = "black", 
                      size = 0.8, linetype = "solid"))+facet_wrap(.~year)

## Warning: Ignoring unknown parameters: shape

Combining the plots ab and c

(c)/(ab)

Wheat Stover plots

### mean and sample size function 
st <- function(stover, upper_limit = max(stover)*1.5) {
  return( 
    data.frame(
      y = 13000,
      label = paste(
                    'Mean =', round(mean(stover), 0), '\n')
    )
  )
}

theme_set(theme_grey(base_size = 10)) ###This sets the font sizes of anything writt
ab<-ggplot(wheat, aes(x = treatment, y = stover))+ 
geom_violin(size=0.9,shape=8)+
  geom_boxplot(size=0.5,width = 0.3,shape=6,fill=c("grey","grey")) +
geom_smooth(method=lm)+ 
ylab("Wheat Stover") + 
xlab("")+
stat_summary(fun.y=mean, geom="point", shape=10, size=3, color="blue", fill="red")+
stat_summary(
    fun.data = st, 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=4
  )+ 
  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=10))+ 
  theme( axis.line = element_line(colour = "black", 
                      size = 0.8, linetype = "solid"))+
  scale_y_continuous(breaks = c(0, 2000,4000,6000,8000,10000,12000,14000))

## Warning: Ignoring unknown parameters: shape

ab

theme_set(theme_grey(base_size = 10)) ###This sets the font sizes of anything writt
ggplot(wheat, aes(x = treatment, y = stover))+ 
geom_violin(size=0.9,shape=8)+
  geom_boxplot(size=0.5,width = 0.3,shape=6,fill=c("grey","grey","grey","grey")) +
geom_smooth(method=lm)+ 
ylab("Wheat Stover") + 
xlab("")+
stat_summary(fun.y=mean, geom="point", shape=10, size=3, color="blue", fill="red")+
stat_summary(
    fun.data = st, 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=4
  )+ 
  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=10))+ 
  theme( axis.line = element_line(colour = "black", 
                      size = 0.8, linetype = "solid"))+
  scale_y_continuous(breaks = c(0, 2000,4000,6000,8000,10000,12000,14000))+
  facet_wrap(.~year)

## Warning: Ignoring unknown parameters: shape

Treatment , grain yield and agroecology

# by agroecology
theme_set(theme_grey(base_size = 10)) ###This sets the font sizes of anything writt
ggplot(wheat, aes(x = treatment, y = grain_yield))+ 
geom_violin(size=0.9,shape=8)+
  geom_boxplot(size=0.5,width = 0.4,shape=6,fill=c("grey","grey","grey","grey")) +
geom_smooth(method=lm)+ 
ylab("Wheat Grain Yield [kg/ha]") + 
xlab("")+
stat_summary(fun.y=mean, geom="point", shape=10, size=3, color="blue", fill="red")+
stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=4
  )+ 
  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=10))+ 
  theme( axis.line = element_line(colour = "black", 
                      size = 0.8, linetype = "solid"))+facet_wrap(.~agroecology)+
  scale_y_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000))

## Warning: Ignoring unknown parameters: shape

Treatment grain yield and soil type

#soil type
theme_set(theme_grey(base_size = 10)) ###This sets the font sizes of anything writt
ggplot(wheat, aes(x = treatment, y = grain_yield))+ 
geom_violin(size=0.9,shape=8)+
  geom_boxplot(size=0.5,width = 0.3,shape=6,fill=c("grey","grey","grey","grey","grey","grey")) +
geom_smooth(method=lm)+ 
ylab("Wheat Grain Yield [kg/ha]") + 
xlab("")+
stat_summary(fun.y=mean, geom="point", shape=10, size=3, color="blue", fill="red")+
stat_summary(
    fun.data = stat_box_data, 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=3
  )+ 
  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=10))+ 
  theme( axis.line = element_line(colour = "black", 
                      size = 0.8, linetype = "solid"))+facet_wrap(.~soil_type)+
  scale_y_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000))

## Warning: Ignoring unknown parameters: shape

Scatter plots with regression output

## general maize yield and plant population graph
my.formula <- grain_yield~ stover

ggplot(wheat, aes(x =stover , y = grain_yield)) +  
  geom_point(size=3,shape=7) + 
  geom_smooth(method=lm, position = "jitter", level = 0.95)+
  stat_poly_eq(formula = my.formula, 
                aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
                parse = TRUE)+
   ylab("Wheat Grain Yield [kg/ha]") + 
  xlab("Stover Yield [kg/ha]")+
  theme(legend.position = c(0.2, 0.85))+
  theme_classic()+
  stat_cor(method = "pearson", label.x = 1500, label.y = 6500)+
  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"))+
  scale_y_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000))

## by treatment
my.formula <- grain_yield~ stover

ggplot(wheat, aes(x =stover , y = grain_yield)) +  
  geom_point(size=3,shape=7) + 
  geom_smooth(method=lm, position = "jitter", level = 0.95)+
  stat_regline_equation(label.x = 1500, label.y = 7000)+
   ylab("Wheat Grain Yield [kg/ha]") + 
  xlab("Stover Yield [kg/ha]")+
  theme(legend.position = c(0.2, 0.85))+
  theme_bw()+
  stat_cor(method = "pearson", label.x = 1500, label.y = 6500)+
  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"))+
  scale_y_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000))+
  facet_wrap(.~treatment)

## by agroecology
my.formula <- grain_yield~ stover

ggplot(wheat, aes(x =stover , y = grain_yield)) +  
  geom_point(size=3,shape=7) + 
  geom_smooth(method=lm, position = "jitter", level = 0.95)+
  stat_regline_equation(label.x = 1500, label.y = 7000)+
   ylab("Wheat Grain Yield [kg/ha]") + 
  xlab("Stover Yield [kg/ha]")+
  theme(legend.position = c(0.2, 0.85))+
  theme_bw()+
  stat_cor(method = "pearson", label.x = 1500, label.y = 6500)+
  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"))+
  scale_y_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000))+
  facet_wrap(.~agroecology)

## by soil type
my.formula <- grain_yield~ stover

ggplot(wheat, aes(x =stover , y = grain_yield)) +  
  geom_point(size=3,shape=7) + 
  geom_smooth(method=lm, position = "jitter", level = 0.95)+
  stat_regline_equation(label.x = 1500, label.y = 7000)+
   ylab("Wheat Grain Yield [kg/ha]") + 
  xlab("Stover Yield [kg/ha]")+
  theme(legend.position = c(0.2, 0.85))+
  theme_bw()+
  stat_cor(method = "pearson", label.x = 1500, label.y = 6500)+
  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"))+
  scale_y_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000))+
  facet_wrap(.~soil_type)

## by season
my.formula <- grain_yield~ stover

ggplot(wheat, aes(x =stover , y = grain_yield)) +  
  geom_point(size=3,shape=7) + 
  geom_smooth(method=lm, position = "jitter", level = 0.95)+
  stat_regline_equation(label.x = 1500, label.y = 7000)+
   ylab("Wheat Grain Yield [kg/ha]") + 
  xlab("Stover Yield [kg/ha]")+
  theme(legend.position = c(0.2, 0.85))+
  theme_bw()+
  stat_cor(method = "pearson", label.x = 1500, label.y = 6500)+
  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"))+
  scale_y_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000))

Grain yield and Harvest Index

my.formula <- grain_yield~ hi

ggplot(wheat, aes(x =hi , y = grain_yield)) +  
  geom_point(size=3,shape=7) + 
  geom_smooth(method=lm, position = "jitter", level = 0.95)+
  stat_regline_equation(label.x = 0.2, label.y = 7000)+
   ylab("Wheat Grain Yield [kg/ha]") + 
  xlab("Harvest Index")+
  theme(legend.position = c(0.2, 0.85))+
  theme_bw()+
  stat_cor(method = "pearson", label.x = 0.2, label.y = 6500)+
  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"))+
  scale_y_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000))

my.formula <- grain_yield~ hi

ggplot(wheat, aes(x =hi , y = grain_yield)) +  
  geom_point(size=3,shape=7) + 
  geom_smooth(method=lm, position = "jitter", level = 0.95)+
  stat_regline_equation(label.x = 0.2, label.y = 7000)+
   ylab("Wheat Grain Yield [kg/ha]") + 
  xlab("Harvest Index")+
  theme(legend.position = c(0.2, 0.85))+
  theme_bw()+
  stat_cor(method = "pearson", label.x = 0.2, label.y = 6500)+
  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"))+
  scale_y_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000))+
  facet_wrap(.~treatment)

Gross Margins

### mean and sample size function 
stat_box<- function(gross_margins, upper_limit = max(gross_margins)*1.5) {
  return( 
    data.frame(
      y = 2800,
      label = paste(
                    'Mean =', round(mean(gross_margins), 0), '\n')
    )
  )
}

theme_set(theme_grey(base_size = 10)) ###This sets the font sizes of anything writt
ggplot(wheat, aes(x = treatment, y = gross_margins))+ 
geom_violin(size=0.9,shape=8)+
  geom_boxplot(size=0.5,width = 0.4,shape=6,fill=c("grey","grey")) +
geom_smooth(method=lm)+ 
ylab("Gross Margins") + 
xlab("")+ 
stat_summary(fun.y=mean, geom="point", shape=10, size=3, color="blue", fill="red")+
stat_summary(
    fun.data = stat_box, 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=4
  )+
  theme_classic(base_size = 12)+
  scale_y_continuous(breaks = c(0,400,800,1200,1600,2000,2400,2800))+
  theme(axis.text.x = element_text(face="bold", color="black", 
                           size=12),
          axis.text.y = element_text(face="bold", color="black", 
                           size=12))+ 
  theme( axis.line = element_line(colour = "black", 
                      size = 0.8, linetype = "solid"))

## Warning: Ignoring unknown parameters: shape

theme_set(theme_grey(base_size = 10)) ###This sets the font sizes of anything writt
ggplot(wheat, aes(x = treatment, y = gross_margins))+ 
geom_violin(size=0.9,shape=8)+
  geom_boxplot(size=0.5,width = 0.4,shape=6,fill=c("grey","grey","grey","grey")) +
geom_smooth(method=lm)+ 
ylab("Gross Margins") + 
xlab("")+ 
stat_summary(fun.y=mean, geom="point", shape=10, size=3, color="blue", fill="red")+
stat_summary(
    fun.data = stat_box, 
    geom = "text", 
    hjust = 0.5,
    vjust = 0.9,
    size=4
  )+
  theme_classic(base_size = 12)+
  scale_y_continuous(breaks = c(0,400,800,1200,1600,2000,2400,2800))+
  theme(axis.text.x = element_text(face="bold", color="black", 
                           size=12),
          axis.text.y = element_text(face="bold", color="black", 
                           size=12))+ 
  theme( axis.line = element_line(colour = "black", 
                      size = 0.8, linetype = "solid"))+
  facet_wrap(.~year)

## Warning: Ignoring unknown parameters: shape

Wheat Mechinization Data Analysis

Chipindu Lovemore [chipindul@rpyanalytics.com]

3/6/2020

Loading the data set

Required Libraries

Descriptive Statistics

Gross Margins Anova Table

Violin plots

Treatment by grain yield

Grain yield by season

Combining the plots ab and c

Wheat Stover plots

Treatment , grain yield and agroecology

Treatment grain yield and soil type

Scatter plots with regression output

Grain yield and Harvest Index

Gross Margins