packages

library(tidyverse)
library(ggpubr)
library(sandwich)
library(rcompanion)

Deviance

y<-c(3,6,4,7)
mean(y)
## [1] 5
plot(y , pch=20,col=2, cex= 3)
points(c(3:4), y[c(3,4)] , 
       pch=20,col= 4, cex= 3) 
abline(h= mean(y), col=3, lty="dashed")
abline( h= mean(c(3,6)), col=2)  # block 1
abline( h= mean(c(4,7)), col=4)   # block 2

Observed Vrs. Expected (null nodel)

TOTAL DEVIANCE (Crawley 2005)

null model [Ha: μ1 = μ2]
2*(3*log(3/5)+6*log(6/5)+
     4*log(4/5)+7*log(7/5))
## [1] 2.048368
#plot
plot(y , pch=20,col=2, cex= 3)
points(c(3:4), y[c(3,4)] , 
       pch=20,col= 4, cex= 3) 
abline(h= mean(y), col=3, lty="dashed")
segments(x0= c(1:4), y0= y, x= c(1:4), y=rep(mean(y),4))

Residual Deviance

2*(3*log(3/4.5)+6*log(6/4.5)+
     4*log(4/5.5)+7*log(7/5.5))
## [1] 1.848033
#plot
plot(y , pch=20,col=2, cex= 3)
points(c(3:4), y[c(3,4)] , 
       pch=20,col= 4, cex= 3) 
abline(h= mean(y), col=3, lty="dashed")
abline( h= mean(c(3,6)), col=2)  # block 1
abline( h= mean(c(4,7)), col=4)   # block 2
segments(x0= c(1:4), y0= y, x= c(1:4), y=c(4.5,4.5,5.5,5.5))

GLM
##   y location
## 1 3        A
## 2 6        A
## 3 4        B
## 4 7        B
glm <- glm(y ~ location, family= poisson)
coef(glm)
## (Intercept)   locationB 
##   1.5040774   0.2006707
exp(coef(glm)[1])   # Expected value for residual Deviance for location "A"
## (Intercept) 
##         4.5
exp(coef(glm)[1] + coef(glm)[2])   # Expected value for resid.Dev for location "B"
## (Intercept) 
##         5.5
glm$deviance        # Residual Deviance
## [1] 1.848033
glm$null.deviance   # Null Deviance
## [1] 2.048368
library(ggplot2)
ggplot(data=dev, aes(x= location, y= y , col= location)) +
  geom_point(size= 4) +
  geom_hline(yintercept = mean(dev$y), col="dark green", linetype="dashed") +
  geom_curve(x= dev$location[1] ,y=4.481689, xend= dev$location[3], yend= 5.47778, 
             curvature = 0.15,
             size=1.3,
             col= "dark gray" , linetype="dashed") +
  geom_errorbar(aes(ymin=c(3,3,4,4), ymax= c(6,6,7,7)), width= 0.2)

EXAMPLE

“slugsurvey.txt”

read.table("Crawley/slugsurvey.txt", header=TRUE, stringsAsFactors=TRUE) -> slug
names(slug)
## [1] "slugs" "field"
head(slug)
##   slugs   field
## 1     3 Nursery
## 2     0 Nursery
## 3     0 Nursery
## 4     0 Nursery
## 5     0 Nursery
## 6     1 Nursery
tail(slug)
##    slugs   field
## 75     3 Rookery
## 76     6 Rookery
## 77     2 Rookery
## 78     4 Rookery
## 79     1 Rookery
## 80     1 Rookery
ggdensity(slug$slugs)

glm “slugsurvey.txt”

glm.slug <- glm(slug$slugs~ slug$field, family = quasipoisson)
summary(glm.slug)
## 
## Call:
## glm(formula = slug$slugs ~ slug$field, family = quasipoisson)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1331  -1.5969  -0.9519   0.4580   4.8727  
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         0.2429     0.2494   0.974   0.3331  
## slug$fieldRookery   0.5790     0.3116   1.858   0.0669 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 3.17311)
## 
##     Null deviance: 224.86  on 79  degrees of freedom
## Residual deviance: 213.44  on 78  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 6

ggplot SLUGS

ggplot(data= slug, aes(x= field, y= slugs, col= field)) + 
  geom_boxplot() +
  geom_point(position = position_dodge2(width = 0.3)) +
  geom_smooth(method = "glm", method.args = list(family = 'poisson')) +
  geom_curve(x= slug$field[1] ,y= exp(coef(glm.slug)[1]), xend=slug$field[41], yend= exp(coef(glm.slug)[1]+coef(glm.slug)[2]),
             curvature = 0.07, col="dark gray", size= 1.2)

AOV approach (General linear model)

library(car)
shapiro.test(slug$slugs)
## 
##  Shapiro-Wilk normality test
## 
## data:  slug$slugs
## W = 0.77782, p-value = 1.187e-09
leveneTest(slugs~ field, data= slug)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.6427 0.4252
##       78
aov.slug <- aov(slugs~ field, data= slug)
summary(aov.slug)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## field        1     20  20.000     3.9 0.0518 .
## Residuals   78    400   5.128                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1