ANts

library(readxl)
aov <- read_excel("wrangled ant.xlsx", sheet = "for ANOVA")

str(aov)

## tibble [56 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Block    : num [1:56] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Nest     : num [1:56] 1 2 3 4 5 6 7 1 2 3 ...
##  $ Treatment: chr [1:56] "A" "A" "A" "A" ...
##  $ Count    : num [1:56] 79 17 1 7 238 23 4 68 112 0 ...
##  $ Location : chr [1:56] "Blacktown" "Blacktown" "Blacktown" "Blacktown" ...
##  $ Counry   : chr [1:56] "Australia" "Australia" "Australia" "Australia" ...

block <- as.factor(aov$Block)
nest <- as.factor(aov$Nest)
treatment <- as.factor(aov$Treatment)

count <- aov$Count
country <- as.factor(aov$Counry)







boxplot(count~treatment)

A = Honey B = Tuna

boxplot(count~country)

Assumptions

Homogenity of Variances

library(GAD)

## Loading required package: matrixStats

## Loading required package: R.methodsS3

## R.methodsS3 v1.8.1 (2020-08-26 16:20:06 UTC) successfully loaded. See ?R.methodsS3 for help.

bartlett.test(count~treatment)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  count by treatment
## Bartlett's K-squared = 0.082576, df = 1, p-value = 0.7738

bartlett.test(count~country)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  count by country
## Bartlett's K-squared = 3.0857, df = 1, p-value = 0.07899

Homogenous variances between each treatment group, but not for each location. ### Normality

shapiro.test(count)

## 
##  Shapiro-Wilk normality test
## 
## data:  count
## W = 0.70122, p-value = 2.23e-09

qqnorm(count
      )

hist(log(count+1))

antslog= log(count+1)

antsqq = lm(antslog ~ country+treatment, data = aov)
par(mfrow=c(2,2))
plot(antsqq)

library(ggplot2)

ANOVA MODELS

Selected ANOVA with blocking. We are using blocking because of two reasons:

Our observations are divided up between observers. Each observer should be a block to make sure that it is still independent observations despite our groupings.
Blocking is used to explain the variation between each grouping (observer) by assessing variation within blocks.

For our experimental design, it is important that each treatment is applied at each block.

# two way with blocking and interaction

blocking <- aov(antslog ~ treatment*country + block)

anova(blocking)

## Analysis of Variance Table
## 
## Response: antslog
##                   Df Sum Sq Mean Sq F value    Pr(>F)    
## treatment          1  0.109  0.1087  0.0671 0.7967087    
## country            1 26.190 26.1898 16.1630 0.0002001 ***
## block              3 10.713  3.5711  2.2039 0.0994391 .  
## treatment:country  1  7.140  7.1397  4.4063 0.0409814 *  
## Residuals         49 79.397  1.6203                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interaction Plot

if(!require(car)){install.packages("car")}

## Loading required package: car

## Loading required package: carData

if(!require(psych)){install.packages("psych")}

## Loading required package: psych

## 
## Attaching package: 'psych'

## The following object is masked from 'package:car':
## 
##     logit

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

if(!require(multcompView)){install.packages("multcompView")}

## Loading required package: multcompView

if(!require(lsmeans)){install.packages("lsmeans")}

## Loading required package: lsmeans

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

if(!require(FSA)){install.packages("FSA")}

## Loading required package: FSA

## ## FSA v0.8.30. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book.

## 
## Attaching package: 'FSA'

## The following object is masked from 'package:psych':
## 
##     headtail

## The following object is masked from 'package:car':
## 
##     bootCase

if(!require(ggplot2)){install.packages("ggplot2")}
if(!require(phia)){install.packages("phia")}

## Loading required package: phia

interaction.plot(x.factor = country,
                 trace.factor = treatment,
                 response = antslog,
                 fun = mean,
                 type = "b",
                 fixed = TRUE,
                 leg.bty = "o")


If you don't want to look at the interaction between treatment and country:

```r
noint <- aov(antslog~treatment+country+block)
anova(noint)

## Analysis of Variance Table
## 
## Response: antslog
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## treatment  1  0.109  0.1087  0.0628 0.8031359    
## country    1 26.190 26.1898 15.1321 0.0002967 ***
## block      3 10.713  3.5711  2.0634 0.1169261    
## Residuals 50 86.537  1.7307                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Post-Hoc

Perform Tukey’s HSD post-hoc test for pairwise comparisons.

tukey <- TukeyHSD(blocking)
tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = antslog ~ treatment * country + block)
## 
## $treatment
##            diff        lwr       upr     p adj
## B-A -0.08811808 -0.7717849 0.5955487 0.7967087
## 
## $country
##                      diff       lwr      upr     p adj
## Taiwan-Australia 1.427228 0.7138232 2.140634 0.0002001
## 
## $block
##           diff        lwr       upr     p adj
## 2-1 -0.4594003 -1.8775447 0.9587441 0.8888409
## 3-1 -0.7295290 -2.1476734 0.6886154 0.5947516
## 4-1 -1.2493699 -2.7419238 0.2431840 0.1408208
## 5-1 -0.1561581 -1.7538416 1.4415254 0.9986681
## 3-2 -0.2701287 -1.7418067 1.2015493 0.9849513
## 4-2 -0.7899696 -2.3334785 0.7535393 0.5994510
## 5-2  0.3032422 -1.3421438 1.9486283 0.9847213
## 4-3 -0.5198409 -2.0633498 1.0236680 0.8741499
## 5-3  0.5733709 -1.0720151 2.2187570 0.8599472
## 5-4  1.0932119 -0.6167235 2.8031472 0.3793378
## 
## $`treatment:country`
##                               diff         lwr       upr     p adj
## B:Australia-A:Australia -0.6203986 -1.74882477 0.5080275 0.4677723
## A:Taiwan-A:Australia     0.6820356 -0.65313616 2.0172074 0.5309673
## B:Taiwan-A:Australia     1.5520226  0.21685077 2.8871944 0.0167500
## A:Taiwan-B:Australia     1.3024343 -0.03273752 2.6376061 0.0583115
## B:Taiwan-B:Australia     2.1724212  0.83724942 3.5075930 0.0004201
## B:Taiwan-A:Taiwan        0.8699869 -0.64395558 2.3839294 0.4288065

We see that treatment B (tuna) is significantly different for Taiwan and Australia (p=0.000420).

(p = 0.0156).

tuk.aov <- aov(antslog ~ treatment:country)
tukey.plot.test <- TukeyHSD(tuk.aov)
plot(tukey.plot.test, las = 2, cex.axis =0.35)

If it does not go below zero, it is significant. We can see that B Taiwan and A Australia are very different. B Taiwan and B Australia are also different. Could this be attributed to different nutritional composition of the tuna selected for the experiment, or due to different species of ants preferring different things?

library(ggplot2)
  
twoplot <- ggplot(blocking, aes (x=country, y= antslog, group=treatment)) +
  geom_point( cex= 1.5, pch = 1.0, position = position_jitter(w =0.1, h = 0))
twoplot