#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).