Deviance Concept

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: μ₁ = μ₂]

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

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.5, xend= dev$location[3], yend= 5.5, 
             curvature = 0.18,
             linewidth=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)

tail(slug)

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)

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

---