QK Figs Box 6.8

Peake and Quinn (1993) investigated the relationship between the number of species of macroinvertebrates, the total abundance of macroinvertebrates, and area of clumps of mussels on a rocky shore in southern Australia. The variables of interest are clump area (dm²), number of species, and number of individuals.

Clump of mussels, with barnacles and gastropods in the clump. Mick Keough

The data set was used in the first edition; get it here Peake, A. J. & Quinn, G. P. (1993). Temporal variation in species-area curves for invertebrates in clumps of an intertidal mussel. Ecography, 16, 269-77.

Preliminaries

First, load the required packages (pwr)

devtools::source_url("https://raw.githubusercontent.com/mjkeough/mjkeough.github.io/refs/heads/main/R/libraries.R")

library(pwr)

Load the data file and transform area, species and indiv

peakquinn <- read_csv("../data/peakquinn.csv")

Rows: 25 Columns: 3── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (3): area, indiv, species
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(peakquinn)

peakquinn$larea<-log10(peakquinn$area)
peakquinn$lspecies<-log10(peakquinn$species)
peakquinn$lindiv<-log10(peakquinn$indiv)

The relationship between species and area

scatterplot (species~area, data=peakquinn)

Problems with distributions, so use logged values

scatterplot (lspecies~larea, data=peakquinn)

pq.lm<-lm(lspecies~larea, data=peakquinn)
plot(pq.lm)

Residuals, etc., look much better now

Look at model fit

options(digits=3)   # Tidy up output
augment(pq.lm)

glance(pq.lm, digits=3)

tidy(pq.lm, conf.int = TRUE, digits=3)

And if you want the traditional ANOVA table…

anova(pq.lm)

Analysis of Variance Table

Response: lspecies
          Df Sum Sq Mean Sq F value  Pr(>F)    
larea      1  1.027    1.03     104 5.1e-10 ***
Residuals 23  0.226    0.01                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Note that all the information is in more compact form in the earlier output

Now look at number of individuals vs clump area

scatterplot (indiv~area, data=peakquinn)

pq2.lm<-lm(indiv~area,data=peakquinn)
plot(pq2.lm)

Scatterplot suggests non-linear relationship and the residual plot shows a pattern. Use log-transformed data

pq3.lm<-lm(lindiv~larea, data=peakquinn)
plot(pq3.lm)

augment(pq3.lm)

glance(pq3.lm)

tidy(pq3.lm, conf.int = TRUE, digits=3)

#And if you want the traditional ANOVA table...
anova(pq3.lm)

Analysis of Variance Table

Response: lindiv
          Df Sum Sq Mean Sq F value Pr(>F)    
larea      1   4.81    4.81     140  3e-11 ***
Residuals 23   0.79    0.03                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Nicer graphic

Fig 6.13

First load script that produces standard graphics appearance (code hidden)

p<-ggplot(peakquinn, aes(x=larea, y=lspecies))+
  geom_point(color=sc)+
  geom_smooth(method="lm", se=FALSE, color=lc)+
  theme_classic(base_size = 10)+
  theme(
    axis.text = element_text(colour = ac),
    axis.line = element_line(color = ac),
    axis.ticks = element_line(color = ac),
        )+
  xlab("Log area")+ylab("Log no. species")
#Add marginal box plots
p1<-ggMarginal(p,size=20, type="boxplot")

`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'

p1

Residual plot

p2<-ggplot(pq.lm, aes(x = pq.lm$fitted.values, y = pq.lm$residuals)) +
  geom_point(color=sc) +
  geom_smooth(se=FALSE, colour=lc)+
  theme_classic(base_size = 10)+
  theme(
    axis.text = element_text(colour = ac),
    axis.line = element_line(color = ac),
    axis.ticks = element_line(color = ac),
        )+labs(x = "Predicted log no. species", y = "Residuals", 
       )
p2

Fig 6.14 - indiv vs area

p3<-ggplot(peakquinn, aes(x=area, y=indiv))+
  geom_point(color = sc)+
  geom_smooth(se = FALSE, color = lc)+
  theme_classic(base_size = 10)+
  theme(
    axis.text = element_text(colour = ac),
    axis.line = element_line(color = ac),
    axis.ticks = element_line(color = ac),
        )+
  xlab("Area")+ylab("No. individuals")
#Add marginal box plots
p4<-ggMarginal(p3,size=20, type="boxplot")

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

p4

Residual plot

p5<-ggplot(pq2.lm, aes(x = pq2.lm$fitted.values, y = pq2.lm$residuals)) +
  geom_point(color=sc) +
  theme_classic(base_size = 10)+
  theme(
    axis.text = element_text(colour = ac),
    axis.line = element_line(color = ac),
    axis.ticks = element_line(color = ac),
        )+labs(x = "Predicted log no. species", y = "Residuals", 
       )
p5

Fig 6.15

p6<-ggplot(peakquinn, aes(x=larea, y=lindiv))+
  geom_point(color=sc)+
  geom_smooth(method="lm", color="black")+
  theme_classic(base_size = 10)+
  theme(
    axis.text = element_text(colour = ac),
    axis.line = element_line(color = ac),
    axis.ticks = element_line(color = ac),
        )+
  xlab("Log area")+ylab("Log no. indiv.")
#Add marginal box plots
p7<-ggMarginal(p6,size=20, type="boxplot")

`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'

p7

Residual plot

p8<-ggplot(pq3.lm, aes(x = pq3.lm$fitted.values, y = pq3.lm$residuals)) +
  geom_point(color=sc) +
  theme_classic(base_size = 10)+
  theme(
    axis.text = element_text(colour = ac),
    axis.line = element_line(color = ac),
    axis.ticks = element_line(color = ac),
        )+labs(x = "Predicted log no. indiv.", y = "Residuals", 
       )
p8