#loading the data set 
library(readr)
intercrop <- read_csv("~/Bledza/intercrop.csv")
## Parsed with column specification:
## cols(
##   Season = col_character(),
##   Site = col_character(),
##   Plot = col_double(),
##   Treatment = col_double(),
##   Subtreatment = col_character(),
##   Hybrid = col_character(),
##   Rep = col_double(),
##   Crop = col_character(),
##   Biomass = col_double(),
##   Total_Biomass = col_double(),
##   Grain_yield = col_double(),
##   Core_yield = col_double(),
##   Plant_population = col_double()
## )
View(intercrop)
attach(intercrop)
names(intercrop)
##  [1] "Season"           "Site"             "Plot"            
##  [4] "Treatment"        "Subtreatment"     "Hybrid"          
##  [7] "Rep"              "Crop"             "Biomass"         
## [10] "Total_Biomass"    "Grain_yield"      "Core_yield"      
## [13] "Plant_population"
Season=as.factor(Season)
Site=as.factor(Site)
Treatment=as.factor(Treatment)
Subtreatment=as.factor(Subtreatment)
Hybrid=as.factor(Hybrid)
Rep=as.factor(Rep)
library(ggplot2)
library(maps)
library(ggalt)
library(extrafontdb)
library(MASS)
library(pscl)
## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis
library(psych)
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(gridExtra)
library(lme4)
## Loading required package: Matrix
library(lmerTest)
## Warning: package 'lmerTest' was built under R version 3.5.3
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(agricolae)
library(psych)

Descriptive Statistics

summary(intercrop)
##     Season              Site                Plot         Treatment     
##  Length:624         Length:624         Min.   : 1.00   Min.   : 1.000  
##  Class :character   Class :character   1st Qu.:13.75   1st Qu.: 5.000  
##  Mode  :character   Mode  :character   Median :29.00   Median : 9.500  
##                                        Mean   :31.64   Mean   : 9.542  
##                                        3rd Qu.:49.00   3rd Qu.:14.000  
##                                        Max.   :72.00   Max.   :18.000  
##                                                                        
##  Subtreatment          Hybrid               Rep           Crop          
##  Length:624         Length:624         Min.   :1.00   Length:624        
##  Class :character   Class :character   1st Qu.:2.00   Class :character  
##  Mode  :character   Mode  :character   Median :2.00   Mode  :character  
##                                        Mean   :2.50                     
##                                        3rd Qu.:3.25                     
##                                        Max.   :4.00                     
##                                                                         
##     Biomass         Total_Biomass       Grain_yield        Core_yield    
##  Min.   :   68.38   Min.   :   68.38   Min.   :    0.0   Min.   :   0.0  
##  1st Qu.: 1224.49   1st Qu.: 2476.77   1st Qu.:  852.1   1st Qu.: 321.4  
##  Median : 2478.44   Median : 7239.49   Median : 3154.1   Median : 807.2  
##  Mean   : 3697.97   Mean   : 8640.26   Mean   : 3782.1   Mean   :1043.8  
##  3rd Qu.: 4396.91   3rd Qu.:11687.42   3rd Qu.: 5957.4   3rd Qu.:1414.3  
##  Max.   :23541.48   Max.   :38653.41   Max.   :13863.5   Max.   :5527.4  
##                     NA's   :34         NA's   :34                        
##  Plant_population 
##  Min.   :  935.8  
##  1st Qu.:32865.8  
##  Median :40382.8  
##  Mean   :38662.9  
##  3rd Qu.:50702.1  
##  Max.   :92903.2  
##  NA's   :72
describeBy(Grain_yield,Hybrid)
## 
##  Descriptive statistics by group 
## group: CZH131001
##    vars  n    mean      sd  median trimmed     mad min      max    range
## X1    1 32 5027.49 3412.29 5961.41 4976.26 3249.06   0 11890.89 11890.89
##     skew kurtosis     se
## X1 -0.18    -1.05 603.21
## -------------------------------------------------------- 
## group: CZH131011
##    vars  n   mean     sd  median trimmed     mad min      max    range
## X1    1 98 4238.9 3564.3 3220.94 3911.34 3690.33   0 13010.84 13010.84
##    skew kurtosis     se
## X1 0.66    -0.71 360.05
## -------------------------------------------------------- 
## group: CZH131014
##    vars  n    mean      sd  median trimmed     mad min     max   range
## X1    1 70 3390.45 3059.14 2440.89 3071.82 2678.97   0 9927.53 9927.53
##    skew kurtosis     se
## X1 0.78    -0.75 365.64
## -------------------------------------------------------- 
## group: PAN413
##    vars  n    mean     sd  median trimmed     mad min     max   range skew
## X1    1 96 2973.66 2221.7 3240.79 2897.25 3257.02   0 7644.11 7644.11 0.14
##    kurtosis     se
## X1    -1.43 226.75
## -------------------------------------------------------- 
## group: SC608
##    vars   n    mean     sd  median trimmed     mad min      max    range
## X1    1 103 3874.86 3210.2 2713.51 3670.72 3598.47   0 10389.13 10389.13
##    skew kurtosis     se
## X1 0.45    -1.24 316.31
## -------------------------------------------------------- 
## group: SC635
##    vars  n   mean      sd  median trimmed     mad min  max range skew
## X1    1 44 2993.1 2434.21 2818.83 2905.45 3447.71   0 7064  7064 0.15
##    kurtosis     se
## X1    -1.62 366.97
## -------------------------------------------------------- 
## group: SC637
##    vars  n    mean      sd  median trimmed     mad min      max    range
## X1    1 53 3861.95 3406.76 3413.01 3569.87 3932.54   0 10741.52 10741.52
##    skew kurtosis     se
## X1 0.65    -0.92 467.95
## -------------------------------------------------------- 
## group: SC727
##    vars  n    mean      sd  median trimmed     mad min      max    range
## X1    1 94 4221.72 3632.34 3388.15 3774.52 3805.24   0 13863.54 13863.54
##    skew kurtosis     se
## X1 0.83    -0.25 374.65
describeBy(Grain_yield,Treatment)
## 
##  Descriptive statistics by group 
## group: 1
##    vars  n    mean      sd  median trimmed     mad min      max    range
## X1    1 32 4536.57 3493.51 3715.11 4335.05 3870.85   0 11890.89 11890.89
##    skew kurtosis     se
## X1 0.43    -1.13 617.57
## -------------------------------------------------------- 
## group: 2
##    vars  n   mean      sd  median trimmed     mad min     max   range skew
## X1    1 35 3483.5 3269.16 2679.21  3241.2 3781.61   0 9927.53 9927.53 0.52
##    kurtosis     se
## X1    -1.22 552.59
## -------------------------------------------------------- 
## group: 3
##    vars  n    mean      sd  median trimmed    mad min     max   range skew
## X1    1 35 3746.23 2989.95 3835.83 3618.58 4239.7   0 8858.83 8858.83 0.26
##    kurtosis     se
## X1     -1.5 505.39
## -------------------------------------------------------- 
## group: 4
##    vars  n    mean      sd median trimmed     mad    min      max   range
## X1    1 32 4523.69 3197.38 3798.1 4355.97 3747.17 402.92 10142.55 9739.63
##    skew kurtosis     se
## X1 0.39    -1.37 565.22
## -------------------------------------------------------- 
## group: 5
##    vars  n   mean      sd  median trimmed     mad min     max   range skew
## X1    1 36 3356.7 3106.29 2314.41 3168.19 3213.77   0 8669.36 8669.36 0.51
##    kurtosis     se
## X1    -1.31 517.72
## -------------------------------------------------------- 
## group: 6
##    vars  n    mean      sd  median trimmed     mad   min      max    range
## X1    1 35 3814.62 3313.29 2892.75  3608.9 3928.34 73.19 10389.13 10315.93
##    skew kurtosis     se
## X1 0.41    -1.36 560.05
## -------------------------------------------------------- 
## group: 7
##    vars  n    mean      sd  median trimmed     mad    min     max   range
## X1    1 31 3739.38 2994.36 3674.05 3423.69 3954.98 263.56 9915.07 9651.52
##    skew kurtosis    se
## X1 0.58    -0.81 537.8
## -------------------------------------------------------- 
## group: 8
##    vars  n    mean      sd  median trimmed     mad min      max    range
## X1    1 34 3185.19 3166.39 2993.12 2763.23 3727.14   0 10741.52 10741.52
##    skew kurtosis     se
## X1 0.84    -0.23 543.03
## -------------------------------------------------------- 
## group: 9
##    vars  n    mean      sd  median trimmed     mad    min     max   range
## X1    1 32 3505.09 2959.73 3198.71 3210.75 3732.21 156.71 9758.76 9602.05
##    skew kurtosis     se
## X1 0.54    -0.93 523.21
## -------------------------------------------------------- 
## group: 10
##    vars  n    mean      sd  median trimmed     mad min      max    range
## X1    1 30 4984.11 3516.29 4120.87  4661.6 3747.11 840 13010.84 12170.84
##    skew kurtosis     se
## X1 0.63    -0.81 641.98
## -------------------------------------------------------- 
## group: 11
##    vars  n    mean      sd  median trimmed     mad min      max    range
## X1    1 35 3752.06 3528.57 2542.68 3436.66 3496.18   0 11421.95 11421.95
##    skew kurtosis     se
## X1 0.67     -0.9 596.44
## -------------------------------------------------------- 
## group: 12
##    vars  n    mean      sd median trimmed    mad    min     max    range
## X1    1 33 4077.78 3642.67 3285.3 3668.85 4037.7 161.03 11924.4 11763.37
##    skew kurtosis     se
## X1 0.69    -0.72 634.11
## -------------------------------------------------------- 
## group: 13
##    vars  n    mean      sd  median trimmed     mad    min      max
## X1    1 28 4602.07 3664.51 4575.12 4300.36 4269.64 391.48 12567.53
##       range skew kurtosis     se
## X1 12176.05 0.69     -0.6 692.53
## -------------------------------------------------------- 
## group: 14
##    vars  n    mean      sd  median trimmed    mad min      max    range
## X1    1 33 4106.87 3557.92 3450.68 3732.75 3777.4   0 11929.12 11929.12
##    skew kurtosis     se
## X1 0.61     -0.7 619.36
## -------------------------------------------------------- 
## group: 15
##    vars  n    mean      sd  median trimmed     mad    min      max
## X1    1 33 4013.84 3764.42 2844.44 3490.01 3220.47 361.03 13863.54
##       range skew kurtosis    se
## X1 13502.52 1.07     0.13 655.3
## -------------------------------------------------------- 
## group: 16
##    vars  n    mean      sd  median trimmed     mad   min     max   range
## X1    1 29 3327.53 2030.37 3946.69 3329.64 1959.09 437.5 6180.56 5743.06
##     skew kurtosis     se
## X1 -0.26    -1.53 377.03
## -------------------------------------------------------- 
## group: 17
##    vars  n    mean      sd  median trimmed     mad min     max   range
## X1    1 33 2514.97 2099.28 1680.07 2415.59 2460.34   0 6129.98 6129.98
##    skew kurtosis     se
## X1 0.25    -1.56 365.44
## -------------------------------------------------------- 
## group: 18
##    vars  n    mean      sd median trimmed     mad    min     max   range
## X1    1 34 3117.03 2469.15   2823 2994.63 3415.75 194.96 7644.11 7449.15
##    skew kurtosis     se
## X1 0.23     -1.5 423.46

YIELD #########################################

fix_model<-lmer(Grain_yield~Treatment+Hybrid +Site+ Plot+Rep+(1|Season), data=intercrop)
summary(fix_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## Grain_yield ~ Treatment + Hybrid + Site + Plot + Rep + (1 | Season)
##    Data: intercrop
## 
## REML criterion at convergence: 10490.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.42334 -0.52067 -0.00042  0.59904  2.79777 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Season   (Intercept) 4444846  2108    
##  Residual             4002774  2001    
## Number of obs: 590, groups:  Season, 4
## 
## Fixed effects:
##                 Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)     2228.868   1177.049    4.596   1.894 0.121906    
## Treatment       -346.201    103.032  574.140  -3.360 0.000831 ***
## HybridCZH131011 2817.709   1020.538  574.267   2.761 0.005947 ** 
## HybridCZH131014 -742.918    468.119  574.569  -1.587 0.113057    
## HybridPAN413    3531.469   1606.535  574.229   2.198 0.028334 *  
## HybridSC608      541.096    526.798  574.401   1.027 0.304786    
## HybridSC635      831.157    778.312  574.332   1.068 0.286015    
## HybridSC637     1031.330    784.887  574.421   1.314 0.189375    
## HybridSC727     3663.461   1310.794  574.245   2.795 0.005366 ** 
## SiteD.T.C       1186.217    283.962  576.830   4.177 3.41e-05 ***
## SiteUZ          5208.217    261.003  575.580  19.955  < 2e-16 ***
## Plot              -5.556      6.089  576.345  -0.912 0.361922    
## Rep               -7.619    101.665  575.290  -0.075 0.940286    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
y=anova(fix_model)
y
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq   Mean Sq NumDF  DenDF  F value    Pr(>F)    
## Treatment   45193027  45193027     1 574.14  11.2904 0.0008308 ***
## Hybrid     143065981  20437997     7 574.25   5.1060 1.207e-05 ***
## Site      1792978420 896489210     2 575.88 223.9670 < 2.2e-16 ***
## Plot         3332455   3332455     1 576.34   0.8325 0.3619216    
## Rep            22481     22481     1 575.29   0.0056 0.9402861    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova_model<-aov(Grain_yield~Treatment+Hybrid +Site+Subtreatment+ Plot+Rep+Season, data=intercrop)
anova_model
## Call:
##    aov(formula = Grain_yield ~ Treatment + Hybrid + Site + Subtreatment + 
##     Plot + Rep + Season, data = intercrop)
## 
## Terms:
##                  Treatment     Hybrid       Site Subtreatment       Plot
## Sum of Squares    17302030  203758179 1638035679     76085964    1353542
## Deg. of Freedom          1          7          2            2          1
##                        Rep     Season  Residuals
## Sum of Squares     3600928 1849006131 2245594939
## Deg. of Freedom          1          3        572
## 
## Residual standard error: 1981.38
## Estimated effects may be unbalanced
## 34 observations deleted due to missingness
###################################################################
h<-HSD.test(anova_model,"Hybrid",group =TRUE)
h
## $statistics
##   MSerror  Df     Mean       CV
##   3925865 572 3782.079 52.38864
## 
## $parameters
##    test name.t ntr StudentizedRange alpha
##   Tukey Hybrid   8          4.30229  0.05
## 
## $means
##           Grain_yield      std   r Min       Max       Q25      Q50
## CZH131001    5027.489 3412.291  32   0 11890.885 2238.9084 5961.407
## CZH131011    4238.899 3564.302  98   0 13010.836 1057.5282 3220.943
## CZH131014    3390.447 3059.138  70   0  9927.533  958.9497 2440.890
## PAN413       2973.660 2221.696  96   0  7644.114  706.3166 3240.795
## SC608        3874.863 3210.201 103   0 10389.126 1200.1617 2713.514
## SC635        2993.103 2434.205  44   0  7064.005  493.5683 2818.832
## SC637        3861.952 3406.760  53   0 10741.517  793.3846 3413.007
## SC727        4221.716 3632.344  94   0 13863.544  912.2385 3388.147
##                Q75
## CZH131001 6913.336
## CZH131011 6932.362
## CZH131014 5182.967
## PAN413    4931.045
## SC608     7069.760
## SC635     5121.421
## SC637     6228.761
## SC727     6222.880
## 
## $comparison
## NULL
## 
## $groups
##           Grain_yield groups
## CZH131001    5027.489      a
## CZH131011    4238.899     ab
## SC727        4221.716     ab
## SC608        3874.863    abc
## SC637        3861.952   abcd
## CZH131014    3390.447    bcd
## SC635        2993.103     cd
## PAN413       2973.660      d
## 
## attr(,"class")
## [1] "group"
bar.group(h$groups,ylim=c(0,7000),density=4,border="blue")

bar.err(h$means,variation="range",ylim=c(0,7000),bar=FALSE,col=0)

######################################################################
t<-HSD.test(anova_model,"Treatment",group =TRUE)
t
## $statistics
##   MSerror  Df     Mean       CV
##   3925865 572 3782.079 52.38864
## 
## $parameters
##    test    name.t ntr StudentizedRange alpha
##   Tukey Treatment  18         4.956669  0.05
## 
## $means
##    Grain_yield      std  r       Min       Max       Q25      Q50      Q75
## 1     4536.569 3493.513 32   0.00000 11890.885 1537.5305 3715.108 6822.262
## 2     3483.502 3269.161 35   0.00000  9927.533  223.3012 2679.215 6298.822
## 3     3746.233 2989.949 35   0.00000  8858.825  919.5455 3835.830 6449.111
## 4     4523.685 3197.380 32 402.91814 10142.550 2048.3806 3798.097 7271.503
## 5     3356.699 3106.292 36   0.00000  8669.364  360.1456 2314.409 5912.156
## 6     3814.623 3313.292 35  73.19351 10389.126  502.3568 2892.748 7271.685
## 7     3739.381 2994.365 31 263.55630  9915.074  916.6386 3674.049 5357.061
## 8     3185.189 3166.390 34   0.00000 10741.517  395.7053 2993.115 4937.052
## 9     3505.085 2959.728 32 156.70893  9758.758  610.9249 3198.709 5576.527
## 10    4984.112 3516.288 30 840.00000 13010.836 2099.3252 4120.871 7861.589
## 11    3752.057 3528.570 35   0.00000 11421.953  731.8476 2542.679 6736.111
## 12    4077.780 3642.666 33 161.03226 11924.399  782.8027 3285.297 6963.066
## 13    4602.067 3664.513 28 391.47660 12567.529 1130.0318 4575.120 6454.020
## 14    4106.872 3557.924 33   0.00000 11929.123  924.7947 3450.679 6280.473
## 15    4013.839 3764.419 33 361.02857 13863.544  774.3590 2844.444 5957.864
## 16    3327.529 2030.368 29 437.50000  6180.564  772.4138 3946.695 5023.208
## 17    2514.973 2099.278 33   0.00000  6129.979  661.1111 1680.068 4650.652
## 18    3117.027 2469.146 34 194.96153  7644.114  544.4444 2823.004 5321.183
## 
## $comparison
## NULL
## 
## $groups
##    Grain_yield groups
## 10    4984.112      a
## 13    4602.067     ab
## 1     4536.569     ab
## 4     4523.685     ab
## 14    4106.872    abc
## 12    4077.780    abc
## 15    4013.839    abc
## 6     3814.623    abc
## 11    3752.057    abc
## 3     3746.233    abc
## 7     3739.381    abc
## 9     3505.085    abc
## 2     3483.502    abc
## 5     3356.699    abc
## 16    3327.529    abc
## 8     3185.189     bc
## 18    3117.027     bc
## 17    2514.973      c
## 
## attr(,"class")
## [1] "group"
bar.group(t$groups,ylim=c(0,10000),density=4,border="blue")

bar.err(t$means,variation="range",ylim=c(0,7000),bar=FALSE,col=0)

####################################################################
s<-HSD.test(anova_model,"Subtreatment",group =TRUE)
s
## $statistics
##   MSerror  Df     Mean       CV
##   3925865 572 3782.079 52.38864
## 
## $parameters
##    test       name.t ntr StudentizedRange alpha
##   Tukey Subtreatment   3         3.323158  0.05
## 
## $means
##             Grain_yield      std   r Min      Max       Q25      Q50
## CBC2           3686.925 3189.556 203   0 13863.54  687.2060 2892.748
## No_cowpea      4347.105 3182.570 181   0 13010.84 1553.4884 3968.114
## Traditional    3379.393 3172.700 206   0 11929.12  499.3448 2440.708
##                  Q75
## CBC2        6045.185
## No_cowpea   6228.761
## Traditional 5593.567
## 
## $comparison
## NULL
## 
## $groups
##             Grain_yield groups
## No_cowpea      4347.105      a
## CBC2           3686.925      b
## Traditional    3379.393      b
## 
## attr(,"class")
## [1] "group"
bar.group(s$groups,ylim=c(0,10000),density=4,border="blue")

bar.err(s$means,variation="range",ylim=c(0,7000),bar=FALSE,col=0)

TOTAL BIOMASS ##################################

#################################################################
##MIXED MODEL 

fix_model<-lmer(Total_Biomass~Treatment+Hybrid +Site+ Plot+Rep+(1|Season), data=intercrop)
summary(fix_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Total_Biomass ~ Treatment + Hybrid + Site + Plot + Rep + (1 |  
##     Season)
##    Data: intercrop
## 
## REML criterion at convergence: 11320.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1684 -0.6063 -0.0312  0.5269  3.8390 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Season   (Intercept) 28339256 5323    
##  Residual             16827352 4102    
## Number of obs: 590, groups:  Season, 4
## 
## Fixed effects:
##                  Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)      5229.853   2870.172     4.007   1.822  0.14239    
## Treatment        -576.597    211.256   574.086  -2.729  0.00654 ** 
## HybridCZH131011  4892.918   2092.534   574.172   2.338  0.01972 *  
## HybridCZH131014   862.802    959.889   574.377   0.899  0.36911    
## HybridPAN413     6057.355   3294.057   574.146   1.839  0.06645 .  
## HybridSC608      1368.966   1080.181   574.262   1.267  0.20554    
## HybridSC635      1033.951   1595.884   574.216   0.648  0.51732    
## HybridSC637      2795.151   1609.388   574.275   1.737  0.08296 .  
## HybridSC727      7389.615   2687.673   574.157   2.749  0.00616 ** 
## SiteD.T.C        -505.865    582.611   576.267  -0.868  0.38561    
## SiteUZ          11632.296    535.294   575.108  21.731  < 2e-16 ***
## Plot              -13.412     12.490   575.745  -1.074  0.28338    
## Rep               -66.509    208.493   574.893  -0.319  0.74984    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
y=anova(fix_model)
y
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq    Mean Sq NumDF  DenDF  F value    Pr(>F)    
## Treatment 1.2536e+08  125355562     1 574.09   7.4495   0.00654 ** 
## Hybrid    6.5225e+08   93177943     7 574.16   5.5373 3.499e-06 ***
## Site      1.1227e+10 5613330990     2 575.98 333.5837 < 2.2e-16 ***
## Plot      1.9402e+07   19401600     1 575.75   1.1530   0.28338    
## Rep       1.7124e+06    1712383     1 574.89   0.1018   0.74984    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#############################################################
anova_model<-aov(Total_Biomass~Treatment+Hybrid +Site+Subtreatment+ Plot+Rep+Season, data=intercrop)
anova_model
## Call:
##    aov(formula = Total_Biomass ~ Treatment + Hybrid + Site + Subtreatment + 
##     Plot + Rep + Season, data = intercrop)
## 
## Terms:
##                   Treatment      Hybrid        Site Subtreatment
## Sum of Squares    126588368   832594826 10246880581    218923565
## Deg. of Freedom           1           7           2            2
##                        Plot         Rep      Season   Residuals
## Sum of Squares      3597296   100216753 10502920007  9537511971
## Deg. of Freedom           1           1           3         572
## 
## Residual standard error: 4083.378
## Estimated effects may be unbalanced
## 34 observations deleted due to missingness
###################################################################
h<-HSD.test(anova_model,"Hybrid",group =TRUE)
h
## $statistics
##    MSerror  Df     Mean       CV
##   16673972 572 8640.255 47.25992
## 
## $parameters
##    test name.t ntr StudentizedRange alpha
##   Tukey Hybrid   8          4.30229  0.05
## 
## $means
##           Total_Biomass      std   r        Min      Max      Q25
## CZH131001     10475.461 6447.703  32 1182.38905 21428.19 4258.580
## CZH131011      8918.295 7615.827  98  502.72870 34013.33 2806.597
## CZH131014      9172.861 9044.117  70   68.37607 33530.69 2941.758
## PAN413         6833.955 4776.305  96  322.66936 19698.65 1954.167
## SC608          8544.548 6515.643 103  336.16673 24324.60 3561.726
## SC635          6397.786 4024.105  44 1148.87299 13787.94 2354.716
## SC637          9141.818 8330.636  53  349.66850 28661.06 2232.616
## SC727         10045.486 8882.917  94  967.16566 38653.41 2805.616
##                 Q50       Q75
## CZH131001 11792.198 14399.539
## CZH131011  6556.786 11762.199
## CZH131014  5966.031 10093.937
## PAN413     6832.440  9644.611
## SC608      6734.929 12110.683
## SC635      6380.349  9568.150
## SC637      7226.479 10816.730
## SC727      8247.398 12288.946
## 
## $comparison
## NULL
## 
## $groups
##           Total_Biomass groups
## CZH131001     10475.461      a
## SC727         10045.486      a
## CZH131014      9172.861      a
## SC637          9141.818      a
## CZH131011      8918.295      a
## SC608          8544.548     ab
## PAN413         6833.955      b
## SC635          6397.786      b
## 
## attr(,"class")
## [1] "group"
bar.group(h$groups,ylim=c(0,18000),density=4,border="blue")

bar.err(h$means,variation="range",ylim=c(0,7000),bar=FALSE,col=0)

######################################################################
t<-HSD.test(anova_model,"Treatment",group =TRUE)
t
## $statistics
##    MSerror  Df     Mean       CV
##   16673972 572 8640.255 47.25992
## 
## $parameters
##    test    name.t ntr StudentizedRange alpha
##   Tukey Treatment  18         4.956669  0.05
## 
## $means
##    Total_Biomass      std  r        Min      Max      Q25      Q50
## 1      11012.718 9192.313 32   68.37607 33530.69 4622.294 8721.596
## 2       8846.766 8249.611 35  366.08236 28940.95 1973.974 7558.907
## 3       9007.748 7554.938 35  546.63671 27375.42 2984.181 7230.542
## 4       9606.445 6542.406 32 1149.43586 23347.10 4766.318 8243.327
## 5       7744.688 6402.296 36  783.26675 23030.08 1499.991 5974.593
## 6       8396.385 6661.696 35  336.16673 24324.60 2175.097 6734.929
## 7       8155.209 6383.593 31  671.06494 23532.33 3016.679 7744.554
## 8       7496.737 7392.720 34  349.66850 28661.06 2185.886 6484.913
## 9       8072.450 6842.896 32  560.06199 27870.61 2350.737 6758.590
## 10     10035.673 7725.797 30 1731.25615 34013.33 5102.513 8187.803
## 11      8223.048 7545.163 35  502.72870 28186.89 1961.497 5161.087
## 12      8639.880 7710.779 33  714.20570 29726.60 2469.543 6342.859
## 13     10822.618 9572.453 28 1363.12467 38653.41 5148.446 8986.590
## 14     10064.063 8640.932 33  967.16566 32171.69 2956.369 8700.897
## 15      9367.524 8737.179 33 1185.30957 36251.96 2618.315 7463.042
## 16      7204.574 4353.726 29 1198.83154 15217.92 2008.200 7790.016
## 17      6077.439 4737.872 33  322.66936 17846.22 1869.064 5514.365
## 18      7252.104 5191.406 34  586.67528 19698.65 2067.044 7146.323
##          Q75
## 1  14398.992
## 2  12751.346
## 3  12335.158
## 4  12379.837
## 5  10605.672
## 6  12932.904
## 7  10551.089
## 8   9079.449
## 9  10286.167
## 10 13570.861
## 11 11263.103
## 12 11714.362
## 13 12266.001
## 14 12798.699
## 15 12056.520
## 16  9489.514
## 17  9344.572
## 18 10880.401
## 
## $comparison
## NULL
## 
## $groups
##    Total_Biomass groups
## 1      11012.718      a
## 13     10822.618     ab
## 14     10064.063     ab
## 10     10035.673     ab
## 4       9606.445    abc
## 15      9367.524    abc
## 3       9007.748    abc
## 2       8846.766    abc
## 12      8639.880    abc
## 6       8396.385    abc
## 11      8223.048    abc
## 7       8155.209    abc
## 9       8072.450    abc
## 5       7744.688    abc
## 8       7496.737    abc
## 18      7252.104     bc
## 16      7204.574     bc
## 17      6077.439      c
## 
## attr(,"class")
## [1] "group"
bar.group(t$groups,ylim=c(0,18000),density=4,border="blue")

bar.err(t$means,variation="range",ylim=c(0,7000),bar=FALSE,col=0)

####################################################################
s<-HSD.test(anova_model,"Subtreatment",group =TRUE)
s
## $statistics
##    MSerror  Df     Mean       CV
##   16673972 572 8640.255 47.25992
## 
## $parameters
##    test       name.t ntr StudentizedRange alpha
##   Tukey Subtreatment   3         3.323158  0.05
## 
## $means
##             Total_Biomass      std   r       Min      Max      Q25
## CBC2             8465.228 7110.723 203 336.16673 36251.96 2415.868
## No_cowpea        9595.596 7529.264 181  68.37607 38653.41 4270.865
## Traditional      7973.333 7287.954 206 322.66936 32171.69 2034.889
##                  Q50      Q75
## CBC2        7119.483 11756.62
## No_cowpea   8479.359 11845.38
## Traditional 5972.182 11062.46
## 
## $comparison
## NULL
## 
## $groups
##             Total_Biomass groups
## No_cowpea        9595.596      a
## CBC2             8465.228      b
## Traditional      7973.333      b
## 
## attr(,"class")
## [1] "group"
bar.group(s$groups,ylim=c(0,18000),density=4,border="blue")

bar.err(s$means,variation="range",ylim=c(0,7000),bar=FALSE,col=0)

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

Graphical analysis

##### Grain yield and plant population

ggplot(intercrop, aes(x = Plant_population, y = Grain_yield, color = Site)) +  
  geom_point(size=0.5, aes(shape=Site)) + 
  geom_smooth(method=lm, position = "jitter", aes(fill=Site), level = 0.95)+ylab("Grain yield [kg/ha]") + 
  xlab("Plant Population")+facet_wrap(.~Hybrid)
## Warning: Removed 72 rows containing non-finite values (stat_smooth).
## Warning: Removed 72 rows containing missing values (geom_point).

###################################################################
ggplot(intercrop, aes(x = Core_yield, y = Grain_yield, color = Site)) +  
  geom_point(size=0.5, aes(shape=Site)) + 
  geom_smooth(method=lm, position = "jitter", aes(fill=Site), level = 0.95)+ylab("Grain yield [kg/ha]") + 
  xlab("Core yield [kg/ha]")+facet_wrap(.~Hybrid)
## Warning: Removed 34 rows containing non-finite values (stat_smooth).
## Warning: Removed 34 rows containing missing values (geom_point).

#### Boxplots Hybrid and Grain yield 
theme_set(theme_gray(base_size =9))
ggplot(intercrop, aes(x = Hybrid, y = Grain_yield)) + 
    geom_boxplot(size=1.2,varwidth = TRUE) + 
    geom_point(data = intercrop, aes(y = mean(Grain_yield))) +
    geom_line(data = intercrop, aes(y =mean(Grain_yield)))+ylab("Maize Grain yield [kg/ha]") + xlab("Hybrids")
## Warning: Removed 34 rows containing non-finite values (stat_boxplot).
## Warning: Removed 624 rows containing missing values (geom_point).
## Warning: Removed 624 rows containing missing values (geom_path).

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

theme_set(theme_gray(base_size =7))
ggplot(intercrop, aes(x = Hybrid, y = Grain_yield)) + 
    geom_boxplot(size=0.7,varwidth = TRUE) + 
    geom_point(data = intercrop, aes(y = mean(Grain_yield))) +
    geom_line(data = intercrop, aes(y =mean(Grain_yield)))+ylab("Maize Grain yield [kg/ha]") + xlab("Hybrids")+facet_wrap(.~Season)
## Warning: Removed 34 rows containing non-finite values (stat_boxplot).
## Warning: Removed 624 rows containing missing values (geom_point).
## Warning: Removed 624 rows containing missing values (geom_path).