library(readr)
GMCCoffstation <- read_csv("GMCCoffstation.csv")
## Parsed with column specification:
## cols(
##   season = col_double(),
##   ward = col_double(),
##   farmer = col_character(),
##   crop = col_character(),
##   treat = col_character(),
##   rep = col_double(),
##   biomass = col_double(),
##   `biomass/1000` = col_double(),
##   grain = col_double(),
##   `grain/1000` = col_double()
## )
View(GMCCoffstation)
attach(GMCCoffstation)
library(lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(ggplot2)
library(psych)
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(agricolae)
require(MASS)
## Loading required package: MASS
library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## 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
library(lubridate) # for working with dates
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
library(scales)   # to access breaks/formatting functions
## 
## Attaching package: 'scales'
## The following objects are masked from 'package:psych':
## 
##     alpha, rescale
## The following object is masked from 'package:readr':
## 
##     col_factor
library(gridExtra) # for arranging plots
library(ggthemes)
library(ggrepel)
library(tidyr)
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:Matrix':
## 
##     expand
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:gridExtra':
## 
##     combine
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## The following objects are masked from 'package:Hmisc':
## 
##     src, summarize
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(dplyr)          # for data manipulation
library(tidyr)          # for data tidying
library(ggplot2)        # for generating the visualizations
library(treemap)
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(ggalluvial)
season=as.numeric(season)
ward=as.factor(ward)
farmer=as.factor(farmer)
treat=as.factor(treat)
rep=as.factor(rep)
######################ANOVA###############
mod<-aov(grain/1000~rep+farmer+ward+season:rep+treat+season+farmer:ward+farmer:treat:ward+season:ward,data=GMCCoffstation)
anova(mod)
## Analysis of Variance Table
## 
## Response: grain/1000
##                    Df  Sum Sq Mean Sq F value    Pr(>F)    
## rep                 1   4.287  4.2872  4.5834  0.033816 *  
## farmer             20 162.584  8.1292  8.6908 < 2.2e-16 ***
## treat               5  26.772  5.3545  5.7244 6.976e-05 ***
## season              1  27.098 27.0982 28.9704 2.612e-07 ***
## rep:season          1   5.067  5.0671  5.4172  0.021208 *  
## ward:season         1   9.128  9.1284  9.7591  0.002123 ** 
## farmer:ward:treat 100  81.343  0.8134  0.8696  0.774441    
## Residuals         158 147.789  0.9354                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod)
##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## rep                 1   4.29   4.287   4.583  0.03382 *  
## farmer             20 162.58   8.129   8.691  < 2e-16 ***
## treat               5  26.77   5.354   5.724 6.98e-05 ***
## season              1  27.10  27.098  28.970 2.61e-07 ***
## rep:season          1   5.07   5.067   5.417  0.02121 *  
## ward:season         1   9.13   9.128   9.759  0.00212 ** 
## farmer:ward:treat 100  81.34   0.813   0.870  0.77444    
## Residuals         158 147.79   0.935                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
y<-HSD.test(mod,"treat",group = TRUE)
y
## $statistics
##     MSerror  Df     Mean       CV       MSD
##   0.9353752 158 2.549194 37.93937 0.5695539
## 
## $parameters
##    test name.t ntr StudentizedRange alpha
##   Tukey  treat   6         4.080022  0.05
## 
## $means
##                grain/1000      std  r    Min    Max      Q25     Q50
## COWP_INT         2.505654 1.221901 48 0.3836 6.0609 1.767800 2.54715
## JCKBN_INT        2.418983 1.008070 48 0.2537 4.9269 1.800750 2.36910
## LABLAB_INT       2.418181 1.292090 48 0.4825 6.5697 1.510825 2.39455
## PGPEA_COWP_INT   2.172323 1.112348 48 0.5073 5.1913 1.366975 1.98030
## PGPEA_INT        2.618533 1.249121 48 0.6387 6.0833 1.640575 2.35075
## SOLE_MZ          3.161487 1.525724 48 0.5355 6.6848 2.090525 2.71220
##                     Q75
## COWP_INT       3.162225
## JCKBN_INT      2.917725
## LABLAB_INT     3.011300
## PGPEA_COWP_INT 2.917450
## PGPEA_INT      3.269675
## SOLE_MZ        4.337925
## 
## $comparison
## NULL
## 
## $groups
##                grain/1000 groups
## SOLE_MZ          3.161487      a
## PGPEA_INT        2.618533     ab
## COWP_INT         2.505654      b
## JCKBN_INT        2.418983      b
## LABLAB_INT       2.418181      b
## PGPEA_COWP_INT   2.172323      b
## 
## attr(,"class")
## [1] "group"
x<-HSD.test(mod,"season",group = TRUE)
x
## $statistics
##     MSerror  Df     Mean       CV      MSD
##   0.9353752 158 2.549194 37.93937 0.330285
## 
## $parameters
##    test name.t ntr StudentizedRange alpha
##   Tukey season   3         3.346043  0.05
## 
## $means
##      grain/1000      std  r    Min    Max      Q25     Q50      Q75
## 2017   2.843977 1.366597 96 0.5073 6.6848 1.778000 2.69950 3.754375
## 2018   2.738111 1.132469 96 0.9399 6.1279 1.937725 2.52270 3.135900
## 2019   2.065493 1.172418 96 0.2537 5.1913 1.082975 1.98715 2.822875
## 
## $comparison
## NULL
## 
## $groups
##      grain/1000 groups
## 2017   2.843977      a
## 2018   2.738111      a
## 2019   2.065493      b
## 
## attr(,"class")
## [1] "group"
v<-HSD.test(mod,"farmer",group = TRUE)
v
## $statistics
##     MSerror  Df     Mean       CV
##   0.9353752 158 2.549194 37.93937
## 
## $parameters
##    test name.t ntr StudentizedRange alpha
##   Tukey farmer  21         5.135691  0.05
## 
## $means
##              grain/1000       std  r    Min    Max      Q25     Q50
## E. Maresa      2.478561 1.3101922 18 0.7193 4.8348 1.350725 2.46155
## E. Matanda     2.656083 0.5783633 12 1.9428 3.8277 2.065325 2.70585
## E. Pangwa      1.964733 0.6196082 12 0.7998 2.8411 1.651825 2.03390
## E.Matanda      1.806967 0.3775486  6 1.4516 2.3965 1.540900 1.66680
## E.Pangwa       1.079467 0.6002911  6 0.5384 2.2410 0.782175 0.93520
## F. Masango     2.564500 0.4976239 12 1.8057 3.4873 2.283400 2.52735
## F.Masango      3.014800 1.1298755  6 1.7877 5.0675 2.398425 2.86675
## J.Ackchomwe    2.232250 1.0063666  6 0.7087 3.1801 1.579375 2.66785
## J.Mukandi      1.817194 0.6088046 18 0.6928 3.1896 1.472750 1.82325
## Jakuchomwe     3.018175 0.6930173 12 2.0529 4.0870 2.447650 3.13120
## L. Kambarami   2.994161 1.9116775 18 0.4278 6.5697 1.534750 2.83300
## L. Masango     2.057233 0.8364770 12 0.9399 3.9311 1.637875 1.90390
## L.Chikaka      3.431733 1.5742558 18 0.7895 6.0833 2.346600 3.30985
## L.Masango      2.594133 1.0972136  6 0.5355 3.8251 2.626550 2.81755
## M. Masango     1.864950 0.9205786 18 0.5073 4.3559 1.244450 1.59845
## M.Chiweshe     3.637050 0.7211141 18 2.0781 4.9142 3.205700 3.69995
## S. Chikauro    4.383767 1.1424939 18 1.7859 6.6848 3.948925 4.36210
## S. Phiri       2.003256 1.0100625 18 0.2537 3.8888 1.190375 2.30385
## S.Manhembe     2.161783 0.8010658 18 0.7849 3.7792 1.684750 1.83680
## S.Vambe        2.492867 1.5323083 18 0.3836 5.3343 1.031150 2.45080
## T. Matanda     1.772089 0.3872743 18 1.0686 2.5193 1.547150 1.80515
##                   Q75
## E. Maresa    3.237825
## E. Matanda   3.026175
## E. Pangwa    2.446125
## E.Matanda    2.037575
## E.Pangwa     1.055675
## F. Masango   2.775075
## F.Masango    3.160025
## J.Ackchomwe  2.881300
## J.Mukandi    2.190225
## Jakuchomwe   3.481000
## L. Kambarami 4.379225
## L. Masango   2.413125
## L.Chikaka    4.440850
## L.Masango    2.959050
## M. Masango   2.297050
## M.Chiweshe   4.161150
## S. Chikauro  4.985475
## S. Phiri     2.602075
## S.Manhembe   2.398525
## S.Vambe      3.030775
## T. Matanda   1.950200
## 
## $comparison
## NULL
## 
## $groups
##              grain/1000 groups
## S. Chikauro    4.383767      a
## M.Chiweshe     3.637050     ab
## L.Chikaka      3.431733     ab
## Jakuchomwe     3.018175     bc
## F.Masango      3.014800    bcd
## L. Kambarami   2.994161    bcd
## E. Matanda     2.656083    bcd
## L.Masango      2.594133    bcd
## F. Masango     2.564500    bcd
## S.Vambe        2.492867    bcd
## E. Maresa      2.478561    bcd
## J.Ackchomwe    2.232250    bcd
## S.Manhembe     2.161783     cd
## L. Masango     2.057233     cd
## S. Phiri       2.003256     cd
## E. Pangwa      1.964733     cd
## M. Masango     1.864950     cd
## J.Mukandi      1.817194     cd
## E.Matanda      1.806967     cd
## T. Matanda     1.772089     cd
## E.Pangwa       1.079467      d
## 
## attr(,"class")
## [1] "group"
#################density plots#################
ggplot(data = GMCCoffstation, mapping = aes(x = grain, fill = ward)) + geom_density(alpha = 0.5)

#################density plots#################
ggplot(data = GMCCoffstation, mapping = aes(x = grain, fill = farmer)) + geom_density(alpha = 0.5)

###############boxplots##############
theme_set(theme_gray(base_size =12))
ggplot(GMCCoffstation, aes(x = treat, y =grain/1000, color = treat)) + 
    geom_boxplot(size=1.2,varwidth = TRUE) + 
    geom_point(data = GMCCoffstation, aes(y = mean(grain/1000))) +
    geom_line(data = GMCCoffstation, aes(y =mean(grain/1000)))+ylab("Grain yield [t/ha] ") + xlab("treatment")+ 
theme( axis.text.x = element_text( angle = 90,  hjust = 1))+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))

###############facets#################
theme_set(theme_gray(base_size =12))
ggplot(GMCCoffstation, aes(x = treat, y =grain/1000, color = treat)) + 
    geom_boxplot(size=1.2,varwidth = TRUE) + 
    geom_point(data = GMCCoffstation, aes(y = mean(grain/1000))) +
    geom_line(data = GMCCoffstation, aes(y =mean(grain/1000)))+ylab("Grain yield [t/ha] ") + xlab("treatment") + facet_wrap(.~season)+ 
theme( axis.text.x = element_text( angle = 90,  hjust = 1))+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7))

################################################
theme_set(theme_gray(base_size =12))
ggplot(GMCCoffstation, aes(x = treat, y =grain/1000, color = treat)) + 
    geom_boxplot(size=1.2,varwidth = TRUE) + 
    geom_point(data = GMCCoffstation, aes(y = mean(grain/1000))) +
    geom_line(data = GMCCoffstation, aes(y =mean(grain/1000)))+ylab("Grain yield [t/ha] ") + xlab("treatment") + facet_wrap(.~farmer)+ 
theme( axis.text.x = element_text( angle = 90,  hjust = 1))+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,2.5,4.5,6.5,7.5))
## 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?

######################################################
theme_set(theme_gray(base_size =12))
ggplot(GMCCoffstation, aes(x = treat, y =grain/1000, color = treat)) + 
    geom_boxplot(size=1.2,varwidth = TRUE) + 
    geom_point(data = GMCCoffstation, aes(y = mean(grain/1000))) +
    geom_line(data = GMCCoffstation, aes(y =mean(grain/1000)))+ylab("Grain yield [t/ha] ") + xlab("treatment") + facet_wrap(.~season+ward)+ 
theme( axis.text.x = element_text( angle = 90,  hjust = 1))+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,1.5,2.5,3.5,4.5,5.5,6.5,7.5))

############biomass/1000###########
theme_set(theme_gray(base_size =15))
m <- ggplot(data=GMCCoffstation,aes(x=treat, y=biomass))
m + geom_violin(size=1.3,shape=8) + geom_boxplot(width=.2, outlier.size=0,fill=c("red","yellow","grey", "blue", "gold", "green"))+ylab("biomass [t/ha]") + xlab("treatment")+ theme(axis.title = element_text(size=20, face="bold"), axis.text.x = element_text(size=12, face="bold", angle = 90, hjust = 1), axis.text.y = element_text(size=16, face="bold"))+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))
## Warning: Ignoring unknown parameters: shape

theme_set(theme_gray(base_size =15))
m <- ggplot(data=GMCCoffstation,aes(x=treat, y=grain))
m + geom_violin(size=1.3,shape=8) + geom_boxplot(width=.2, outlier.size=0,fill=c("red","yellow","grey", "blue", "gold", "green"))+ylab("grain [t/ha]") +stat_summary(fun.y=median, geom="point", size=2, color="red")+ xlab("treatment")+ geom_jitter(shape=16, position=position_jitter(0.2))+ theme(axis.title = element_text(size=20, face="bold"), axis.text.x = element_text(size=12, face="bold", angle = 90, hjust = 1), axis.text.y = element_text(size=16, face="bold"))+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))
## Warning: Ignoring unknown parameters: shape

#####################facet###################
theme_set(theme_gray(base_size =12))
ggplot(GMCCoffstation, aes(x = treat, y =biomass/1000, color = treat)) + 
    geom_boxplot(size=1.2,varwidth = TRUE) + 
    geom_point(data = GMCCoffstation, aes(y = mean(biomass/1000))) +
    geom_line(data = GMCCoffstation, aes(y =mean(biomass/1000)))+ylab("Biomass [t/ha] ") + xlab("treatment") + facet_wrap(.~season)+ 
theme( axis.text.x = element_text( angle = 90,  hjust = 1))+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7))

ggplot(GMCCoffstation, aes(treat, grain/1000)) + geom_boxplot(fill = "springgreen3") + theme_bw() +
  xlab("") + ylab("Yield (t/ha)") + 
  theme(axis.title = element_text(size=24, face="bold"), 
        axis.text.x = element_text(size=12, face="bold", angle = 90, hjust = 1),
        axis.text.y = element_text(size=20, face="bold"))+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))

ggplot(GMCCoffstation, aes(x = treat, y = grain, fill = treat)) + 
  geom_violin() + scale_y_log10() + geom_boxplot(width = 0.2) + theme(axis.title = element_text(size=20, face="bold"), axis.text.x = element_text(size=12, face="bold", angle = 90, hjust = 1), axis.text.y = element_text(size=16, face="bold")) + theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

#############basic ggplot###################
theme_set(theme_gray(base_size =15))
p<-ggplot(GMCCoffstation, aes(x = treat, y = grain, fill = treat)) + geom_boxplot(aes(fill = treat) + scale_y_log10())+ theme( axis.text.x = element_text( angle = 90,  hjust = 1))+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))
p

#################using facet grid##################
###########split in vertical direction##################
theme_set(theme_gray(base_size =15))
p + facet_grid(rows = vars(season)) + theme(legend.position = "none")+ scale_y_continuous(breaks = c(500,1500,2500,3500,4500,5500,6500))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

theme_set(theme_gray(base_size =15))
p + facet_grid(rows = vars(ward)) + theme(legend.position = "none")+ scale_y_continuous(breaks = c(500,1500,2500,3500,4500,5500,6500))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

###############split in horizontal direction#############
p + facet_grid(cols = vars(season), scales= "free", space= "free")+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

p + facet_grid(cols = vars(ward), scales= "free", space= "free")+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

p + facet_grid(cols = vars(season), scales= "free", space= "free")+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

#############basic biomass ggplot###################
theme_set(theme_gray(base_size =15))
p<-ggplot(GMCCoffstation, aes(x = treat, y = biomass, fill = treat)) + geom_boxplot(aes(fill = treat) + scale_y_log10())+ theme( axis.text.x = element_text( angle = 90,  hjust = 1))+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))
p

#################using facet grid##################
###########split in vertical direction##################
theme_set(theme_gray(base_size =15))
p + facet_grid(rows = vars(season)) + theme(legend.position = "none")+ scale_y_continuous(breaks = c(500,1500,2500,3500,4500,5500,6500))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

theme_set(theme_gray(base_size =15))
p + facet_grid(rows = vars(ward)) + theme(legend.position = "none")+ scale_y_continuous(breaks = c(500,1500,2500,3500,4500,5500,6500))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

###############split in horizontal direction#############
p + facet_grid(cols = vars(season), scales= "free", space= "free")+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

p + facet_grid(cols = vars(ward), scales= "free", space= "free")+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

p + facet_grid(vars(season), vars(ward))+ scale_y_continuous(breaks = c(500,1500,2500,3500,4500,5500,6500))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

p + facet_grid(cols = vars(season), scales= "free", space= "free")+ theme(legend.position = "none")+ scale_y_continuous(breaks = c(0,500,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

ggplot() + geom_col(data = GMCCoffstation, aes(x = treat, y = grain))+ theme_classic(base_size = 20)+ theme(axis.title = element_text(size=20, face="bold"), axis.text.x = element_text(size=12, face="bold", angle = 75, hjust = 1), axis.text.y = element_text(size=16, face="bold"))

ggplot(GMCCoffstation, aes(x=treat, y=grain, fill=treat))+
  geom_bar(stat='identity', fill="forest green")+
  ylab("grain")+ theme_classic(base_size = 20)+ theme(axis.title = element_text(size=20, face="bold"), axis.text.x = element_text(size=12, face="bold", angle = 75, hjust = 1), axis.text.y = element_text(size=16, face="bold"))+
  facet_wrap(~season)

#########################################3
ggplot(GMCCoffstation, aes(x=treat, y=grain, fill=treat))+
  geom_bar(stat='identity', fill="forest green")+
  ylab("grain")+ theme_classic(base_size = 20)+ theme(axis.title = element_text(size=20, face="bold"), axis.text.x = element_text(size=12, face="bold", angle = 75, hjust = 1), axis.text.y = element_text(size=16, face="bold"))+
  facet_wrap(~season,  ncol=1)

###############################3
ggplot(GMCCoffstation, aes(x=treat, y=grain, fill=treat))+
  geom_bar(stat='identity', fill="forest green")+
  ylab("grain")+ theme_classic(base_size = 20)+ theme(axis.title = element_text(size=20, face="bold"), axis.text.x = element_text(size=12, face="bold", angle = 75, hjust = 1), axis.text.y = element_text(size=16, face="bold"))+
  facet_wrap(~season,  ncol=1, strip.position = "left")

###################GEOM PLOTS #########################
ggplot(data = GMCCoffstation, aes(x=season,y=grain/1000,color=treat))+
  geom_point()+ geom_smooth()   
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401

#Simple alluvial diagram

ggplot(GMCCoffstation, aes(kg = grain/1000, axis1 = treat, axis2 = season)) +  
  geom_stratum() + geom_alluvium(aes(fill = treat), color = "black")

#####################
ggplot(GMCCoffstation, aes(kg = grain/1000, axis1 = treat, axis2 = season)) +   
  geom_alluvium(aes(fill = treat), color = "black") + geom_stratum() 

#Add text

ggplot(GMCCoffstation, aes(kg = grain/1000, axis1 = treat, axis2 = season)) +  
  geom_alluvium(aes(fill = treat), color = "black") + geom_stratum() +
  geom_text(stat = "stratum", label.strata = TRUE) +
  theme(legend.position = "true",
        axis.text = element_text(size=14)) 

###################GEOM PLOTS #########################
ggplot(data = GMCCoffstation, aes(x=season,y=grain/1000,color=treat))+
  geom_point()+ geom_smooth()   
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 6.7241e-017
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 2017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 6.7241e-017
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 4.0401