data=read.table("food.txt",head=TRUE)
str(data)
## 'data.frame': 144 obs. of 5 variables:
## $ nema : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
## $ fungus: int 1 1 1 1 1 1 2 2 2 2 ...
## $ exp : int 1 1 1 1 1 1 1 1 1 1 ...
## $ rep : int 1 2 3 4 5 6 1 2 3 4 ...
## $ final : num 133 264 302 14 368 ...
data$fungus = as.factor(data$fungus)
data$exp = as.factor(data$exp)
data$rep = as.factor(data$rep)
attach(data)
data
## nema fungus exp rep final
## 1 A 1 1 1 133.0000000
## 2 A 1 1 2 264.0000000
## 3 A 1 1 3 302.0000000
## 4 A 1 1 4 14.0000000
## 5 A 1 1 5 368.0000000
## 6 A 1 1 6 193.0000000
## 7 A 2 1 1 0.9259259
## 8 A 2 1 2 5.5555556
## 9 A 2 1 3 30.6666667
## 10 A 2 1 4 0.0000000
## 11 A 2 1 5 36.1111111
## 12 A 2 1 6 14.1666667
## 13 A 3 1 1 831.0000000
## 14 A 3 1 2 1111.0000000
## 15 A 3 1 3 616.0000000
## 16 A 3 1 4 3837.0000000
## 17 A 3 1 5 2240.0000000
## 18 A 3 1 6 412.6666667
## 19 A 4 1 1 1227.0000000
## 20 A 4 1 2 1822.2222220
## 21 A 4 1 3 2442.5925930
## 22 A 4 1 4 956.0000000
## 23 A 4 1 5 914.6666667
## 24 A 4 1 6 1078.6666670
## 25 A 5 1 1 5.0000000
## 26 A 5 1 2 46.0000000
## 27 A 5 1 3 3.0000000
## 28 A 5 1 4 0.9259259
## 29 A 5 1 5 8.0000000
## 30 A 5 1 6 32.6666667
## 31 A 6 1 1 820.5000000
## 32 A 6 1 2 1570.5000000
## 33 A 6 1 3 691.5000000
## 34 A 6 1 4 923.0000000
## 35 A 6 1 5 712.0000000
## 36 A 6 1 6 897.0000000
## 37 B 1 1 1 879.0849673
## 38 B 1 1 2 1452.6666670
## 39 B 1 1 3 2506.9444440
## 40 B 1 1 4 1850.0000000
## 41 B 1 1 5 311.7647059
## 42 B 1 1 6 1514.6666670
## 43 B 2 1 1 1070.6666670
## 44 B 2 1 2 1245.3333330
## 45 B 2 1 3 1388.6666670
## 46 B 2 1 4 636.0000000
## 47 B 2 1 5 575.3333333
## 48 B 2 1 6 1140.5228760
## 49 B 3 1 1 336.6666667
## 50 B 3 1 2 474.6666667
## 51 B 3 1 3 1308.4848480
## 52 B 3 1 4 1177.6666670
## 53 B 3 1 5 1735.6666670
## 54 B 3 1 6 1684.6666670
## 55 B 4 1 1 451.6666667
## 56 B 4 1 2 401.3333333
## 57 B 4 1 3 2022.6666670
## 58 B 4 1 4 1079.3333330
## 59 B 4 1 5 1075.3333330
## 60 B 4 1 6 1535.3333330
## 61 B 5 1 1 15.3333333
## 62 B 5 1 2 14.0000000
## 63 B 5 1 3 20.0000000
## 64 B 5 1 4 6.6666667
## 65 B 5 1 5 35.3333333
## 66 B 5 1 6 52.0000000
## 67 B 6 1 1 400.3333333
## 68 B 6 1 2 NA
## 69 B 6 1 3 1151.6666670
## 70 B 6 1 4 1474.3333330
## 71 B 6 1 5 1986.6666670
## 72 B 6 1 6 1094.3333330
## 73 A 1 2 1 60.3448276
## 74 A 1 2 2 10.8333333
## 75 A 1 2 3 59.1666667
## 76 A 1 2 4 158.3333333
## 77 A 1 2 5 180.4687500
## 78 A 1 2 6 80.0000000
## 79 A 2 2 1 0.0000000
## 80 A 2 2 2 1.7241379
## 81 A 2 2 3 14.6551724
## 82 A 2 2 4 0.0000000
## 83 A 2 2 5 6.0344828
## 84 A 2 2 6 0.0000000
## 85 A 3 2 1 806.0000000
## 86 A 3 2 2 587.0000000
## 87 A 3 2 3 326.5000000
## 88 A 3 2 4 430.5000000
## 89 A 3 2 5 897.0000000
## 90 A 3 2 6 1802.0000000
## 91 A 4 2 1 1085.0000000
## 92 A 4 2 2 2511.1111110
## 93 A 4 2 3 349.0740741
## 94 A 4 2 4 903.0000000
## 95 A 4 2 5 308.6206897
## 96 A 4 2 6 2402.6666670
## 97 A 5 2 1 9.0000000
## 98 A 5 2 2 9.0000000
## 99 A 5 2 3 26.0000000
## 100 A 5 2 4 28.0000000
## 101 A 5 2 5 120.0000000
## 102 A 5 2 6 NA
## 103 A 6 2 1 854.0000000
## 104 A 6 2 2 154.5000000
## 105 A 6 2 3 882.5000000
## 106 A 6 2 4 600.0000000
## 107 A 6 2 5 914.1666667
## 108 A 6 2 6 1720.8333330
## 109 B 1 2 1 743.7908497
## 110 B 1 2 2 166.0000000
## 111 B 1 2 3 1606.0000000
## 112 B 1 2 4 1244.0000000
## 113 B 1 2 5 592.6666667
## 114 B 1 2 6 2080.0000000
## 115 B 2 2 1 57.3333333
## 116 B 2 2 2 283.0065359
## 117 B 2 2 3 3450.0000000
## 118 B 2 2 4 2534.0000000
## 119 B 2 2 5 1628.6666670
## 120 B 2 2 6 1078.6666670
## 121 B 3 2 1 218.0000000
## 122 B 3 2 2 575.6666667
## 123 B 3 2 3 798.6666667
## 124 B 3 2 4 400.3333333
## 125 B 3 2 5 1320.6666670
## 126 B 3 2 6 1419.3333330
## 127 B 4 2 1 732.2404372
## 128 B 4 2 2 576.0000000
## 129 B 4 2 3 1672.0000000
## 130 B 4 2 4 1546.6666670
## 131 B 4 2 5 675.3333333
## 132 B 4 2 6 859.3333333
## 133 B 5 2 1 30.0000000
## 134 B 5 2 2 64.0000000
## 135 B 5 2 3 48.6666667
## 136 B 5 2 4 48.0000000
## 137 B 5 2 5 NA
## 138 B 5 2 6 NA
## 139 B 6 2 1 880.0000000
## 140 B 6 2 2 460.6666667
## 141 B 6 2 3 1217.6666670
## 142 B 6 2 4 2590.3333330
## 143 B 6 2 5 708.6666667
## 144 B 6 2 6 NA
boxplot(data$final~data$fungus)
boxplot(data$final~data$nema)
data$logfinal=log(data$final+1)
interaction_exp<- aov(logfinal~exp+nema+fungus+nema*fungus, data=data)
summary(interaction_exp)
## Df Sum Sq Mean Sq F value Pr(>F)
## exp 1 0.39 0.39 0.49 0.485
## nema 1 68.07 68.07 86.36 5.74e-16 ***
## fungus 5 298.33 59.67 75.70 < 2e-16 ***
## nema:fungus 5 132.95 26.59 33.74 < 2e-16 ***
## Residuals 126 99.31 0.79
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 5 observations deleted due to missingness
plot(interaction_exp)
data=data[-c(101,28,115),]
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.
interaction_exp<- aov(logfinal~exp+nema+fungus+nema*fungus, data=data)
summary(interaction_exp)
## Df Sum Sq Mean Sq F value Pr(>F)
## exp 1 0.70 0.70 1.038 0.31
## nema 1 63.70 63.70 95.044 <2e-16 ***
## fungus 5 285.05 57.01 85.067 <2e-16 ***
## nema:fungus 5 139.43 27.89 41.610 <2e-16 ***
## Residuals 123 82.43 0.67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 5 observations deleted due to missingness
plot(interaction_exp)
data=data[-c(11,74,4),]
library(lsmeans)
## Warning: package 'lsmeans' was built under R version 3.4.4
## Loading required package: emmeans
## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.
lsmeans(interaction_exp, ~fungus)
## NOTE: Results may be misleading due to involvement in interactions
## fungus lsmean SE df lower.CL upper.CL
## 1 5.76 0.167 123 5.43 6.09
## 2 4.26 0.171 123 3.92 4.59
## 3 6.72 0.167 123 6.39 7.05
## 4 6.92 0.167 123 6.59 7.25
## 5 2.99 0.188 123 2.62 3.37
## 6 6.80 0.175 123 6.46 7.15
##
## Results are averaged over the levels of: exp, nema
## Confidence level used: 0.95
lsmeans(interaction_exp, pairwise~fungus, adjust='Tukey')
## NOTE: Results may be misleading due to involvement in interactions
## $lsmeans
## fungus lsmean SE df lower.CL upper.CL
## 1 5.76 0.167 123 5.43 6.09
## 2 4.26 0.171 123 3.92 4.59
## 3 6.72 0.167 123 6.39 7.05
## 4 6.92 0.167 123 6.59 7.25
## 5 2.99 0.188 123 2.62 3.37
## 6 6.80 0.175 123 6.46 7.15
##
## Results are averaged over the levels of: exp, nema
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## 1 - 2 1.5021 0.239 123 6.284 <.0001
## 1 - 3 -0.9667 0.236 123 -4.090 0.0011
## 1 - 4 -1.1642 0.236 123 -4.926 <.0001
## 1 - 5 2.7642 0.252 123 10.977 <.0001
## 1 - 6 -1.0475 0.242 123 -4.326 0.0004
## 2 - 3 -2.4687 0.239 123 -10.329 <.0001
## 2 - 4 -2.6663 0.239 123 -11.155 <.0001
## 2 - 5 1.2621 0.254 123 4.965 <.0001
## 2 - 6 -2.5496 0.245 123 -10.415 <.0001
## 3 - 4 -0.1976 0.236 123 -0.836 0.9602
## 3 - 5 3.7309 0.252 123 14.815 <.0001
## 3 - 6 -0.0808 0.242 123 -0.334 0.9994
## 4 - 5 3.9285 0.252 123 15.600 <.0001
## 4 - 6 0.1168 0.242 123 0.482 0.9967
## 5 - 6 -3.8117 0.257 123 -14.814 <.0001
##
## Results are averaged over the levels of: exp, nema
## P value adjustment: tukey method for comparing a family of 6 estimates
library(agricolae)
## Warning: package 'agricolae' was built under R version 3.4.4
(HSD.test(interaction_exp,"fungus"))
## $statistics
## MSerror Df Mean CV
## 0.67019 123 5.646431 14.49856
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey fungus 6 4.094364 0.05
##
## $means
## logfinal std r Min Max Q25 Q50 Q75
## 1 5.757406 1.5042929 24 2.470920 7.827219 5.027709 5.729592 7.165629
## 2 4.140739 3.0233541 23 0.000000 8.146419 1.441232 3.613916 7.012263
## 3 6.724075 0.7055343 24 5.389072 8.252707 6.140355 6.708578 7.204648
## 4 6.921650 0.6080119 24 5.735348 7.828879 6.577277 6.982862 7.363971
## 5 3.029550 0.8303066 19 1.386294 4.174387 2.302585 3.295837 3.721442
## 6 6.792773 0.5955118 22 5.046646 7.859928 6.565967 6.792031 7.091593
##
## $comparison
## NULL
##
## $groups
## logfinal groups
## 4 6.921650 a
## 6 6.792773 a
## 3 6.724075 a
## 1 5.757406 b
## 2 4.140739 c
## 5 3.029550 d
##
## attr(,"class")
## [1] "group"
library(agricolae)
(HSD.test(interaction_exp,"nema"))
## $statistics
## MSerror Df Mean CV
## 0.67019 123 5.646431 14.49856
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey nema 2 2.799349 0.05
##
## $means
## logfinal std r Min Max Q25 Q50 Q75
## A 4.971061 2.344447 69 0.000000 8.252707 3.295837 5.791488 6.806829
## B 6.341961 1.423132 67 2.036882 8.146419 5.996037 6.781058 7.289242
##
## $comparison
## NULL
##
## $groups
## logfinal groups
## B 6.341961 a
## A 4.971061 b
##
## attr(,"class")
## [1] "group"
library(agricolae)
(HSD.test(interaction_exp,"exp"))
## $statistics
## MSerror Df Mean CV
## 0.67019 123 5.646431 14.49856
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey exp 2 2.799349 0.05
##
## $means
## logfinal std r Min Max Q25 Q50 Q75
## 1 5.715874 2.046272 70 0 8.252707 4.202179 6.640307 7.124256
## 2 5.572780 2.081780 66 0 8.146419 4.229403 6.381522 6.988794
##
## $comparison
## NULL
##
## $groups
## logfinal groups
## 1 5.715874 a
## 2 5.572780 a
##
## attr(,"class")
## [1] "group"
library(summariser)
## Warning: package 'summariser' was built under R version 3.4.4
## Loading required package: dplyr
## Warning: package 'dplyr' was built under R version 3.4.4
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Loading required package: ggplot2
## Loading required package: lazyeval
## Warning: package 'lazyeval' was built under R version 3.4.4
## Loading required package: plotrix
## Warning: package 'plotrix' was built under R version 3.4.3
nema1<- summarise(group_by(data,fungus),
mu = mean(final),
sd = sd(final),
n = length(final),
se = sd/sqrt(n),
ciu = mu + (qt(0.025, df = n-1)* se),
cil = mu - (qt(0.025, df = n-1)* se))
## Warning: package 'bindrcpp' was built under R version 3.4.4
nema1
## # A tibble: 6 x 7
## fungus mu sd n se ciu cil
## <fct> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 1 759. 759. 22 162. 422. 1096.
## 2 2 687. 939. 22 200. 270. 1103.
## 3 3 1056. 808. 24 165. 715. 1397.
## 4 4 1193. 667. 24 136. 911. 1474.
## 5 5 NA NaN 22 NaN NA NA
## 6 6 NA NaN 24 NaN NA NA
box1 <- ggplot(data=data , aes(fungus, nema, fill = final))+
geom_jitter(width = 0.05) +
geom_boxplot()+
labs(x ="Fungi", y = "Nematoide", title = "Feeding preference of Foliar nematode") +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
box1
ggsave("box_fresh_nema.png", width = 7, height=5, dpi = 300)