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
(dm2), 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

LS0tCnRpdGxlOiAiUUsgRmlncyBCb3ggNi44IgoKb3V0cHV0OiAKICBodG1sX25vdGVib29rOgogICAgdGhlbWU6IGZsYXRseQotLS0KCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQoZWNobyA9IFRSVUUpCmBgYAoKUGVha2UgYW5kIFF1aW5uICgxOTkzKSBpbnZlc3RpZ2F0ZWQgdGhlIHJlbGF0aW9uc2hpcCBiZXR3ZWVuIHRoZSBudW1iZXIgb2Ygc3BlY2llcyBvZiBtYWNyb2ludmVydGVicmF0ZXMsIHRoZSB0b3RhbCBhYnVuZGFuY2Ugb2YgbWFjcm9pbnZlcnRlYnJhdGVzLCBhbmQgYXJlYSBvZiBjbHVtcHMgb2YgbXVzc2VscyBvbiBhIHJvY2t5IHNob3JlIGluIHNvdXRoZXJuIEF1c3RyYWxpYS4gVGhlIHZhcmlhYmxlcyBvZiBpbnRlcmVzdCBhcmUgY2x1bXAgYXJlYSAoZG1eMl4pLCBudW1iZXIgb2Ygc3BlY2llcywgYW5kIG51bWJlciBvZiBpbmRpdmlkdWFscy4KCiFbXSguLi9tZWRpYS9tdXNzZWxzLmpwZyl7d2lkdGg9IjgwMCJ9CgpDbHVtcCBvZiBtdXNzZWxzLCB3aXRoIGJhcm5hY2xlcyBhbmQgZ2FzdHJvcG9kcyBpbiB0aGUgY2x1bXAuIE1pY2sgS2VvdWdoIFshW10oaW1hZ2VzL2J5LTAzLnBuZyl7d2lkdGg9IjU3In1dKCMwKQoKVGhlIGRhdGEgc2V0IHdhcyB1c2VkIGluIHRoZSBmaXJzdCBlZGl0aW9uOyBnZXQgaXQgW2hlcmVdKGRhdGEvcGVha3F1aW5uLmNzdikgUGVha2UsIEEuIEouICYgUXVpbm4sIEcuIFAuICgxOTkzKS4gVGVtcG9yYWwgdmFyaWF0aW9uIGluIHNwZWNpZXMtYXJlYSBjdXJ2ZXMgZm9yIGludmVydGVicmF0ZXMgaW4gY2x1bXBzIG9mIGFuIGludGVydGlkYWwgbXVzc2VsLiAqRWNvZ3JhcGh5KiwgMTYsIDI2OS03Ny4KCiMjIyBQcmVsaW1pbmFyaWVzCgpGaXJzdCwgbG9hZCB0aGUgcmVxdWlyZWQgcGFja2FnZXMgKHB3cikKCmBgYHtyIHJlc3VsdHM9J2hpZGUnfQpkZXZ0b29sczo6c291cmNlX3VybCgiaHR0cHM6Ly9yYXcuZ2l0aHVidXNlcmNvbnRlbnQuY29tL21qa2VvdWdoL21qa2VvdWdoLmdpdGh1Yi5pby9yZWZzL2hlYWRzL21haW4vUi9saWJyYXJpZXMuUiIpCmxpYnJhcnkocHdyKQpgYGAKCkxvYWQgdGhlIGRhdGEgZmlsZSBhbmQgdHJhbnNmb3JtIGFyZWEsIHNwZWNpZXMgYW5kIGluZGl2CgpgYGB7cn0KcGVha3F1aW5uIDwtIHJlYWRfY3N2KCIuLi9kYXRhL3BlYWtxdWlubi5jc3YiKQpoZWFkKHBlYWtxdWlubikKcGVha3F1aW5uJGxhcmVhPC1sb2cxMChwZWFrcXVpbm4kYXJlYSkKcGVha3F1aW5uJGxzcGVjaWVzPC1sb2cxMChwZWFrcXVpbm4kc3BlY2llcykKcGVha3F1aW5uJGxpbmRpdjwtbG9nMTAocGVha3F1aW5uJGluZGl2KQpgYGAKCiMjIyBUaGUgcmVsYXRpb25zaGlwIGJldHdlZW4gc3BlY2llcyBhbmQgYXJlYQoKYGBge3J9CnNjYXR0ZXJwbG90IChzcGVjaWVzfmFyZWEsIGRhdGE9cGVha3F1aW5uKQpgYGAKClByb2JsZW1zIHdpdGggZGlzdHJpYnV0aW9ucywgc28gdXNlIGxvZ2dlZCB2YWx1ZXMKCmBgYHtyfQpzY2F0dGVycGxvdCAobHNwZWNpZXN+bGFyZWEsIGRhdGE9cGVha3F1aW5uKQpwcS5sbTwtbG0obHNwZWNpZXN+bGFyZWEsIGRhdGE9cGVha3F1aW5uKQpwbG90KHBxLmxtKQpgYGAKUmVzaWR1YWxzLCBldGMuLCBsb29rIG11Y2ggYmV0dGVyIG5vdwoKTG9vayBhdCBtb2RlbCBmaXQKYGBge3J9CmF1Z21lbnQocHEubG0pCmdsYW5jZShwcS5sbSwgZGlnaXRzPTMpCnRpZHkocHEubG0sIGNvbmYuaW50ID0gVFJVRSwgZGlnaXRzPTMpCmBgYAoKQW5kIGlmIHlvdSB3YW50IHRoZSB0cmFkaXRpb25hbCBBTk9WQSB0YWJsZS4uLgoKYGBge3J9CmFub3ZhKHBxLmxtKQpgYGAKCk5vdGUgdGhhdCBhbGwgdGhlIGluZm9ybWF0aW9uIGlzIGluIG1vcmUgY29tcGFjdCBmb3JtIGluIHRoZSBlYXJsaWVyIG91dHB1dAoKCiMjIyBOb3cgbG9vayBhdCBudW1iZXIgb2YgaW5kaXZpZHVhbHMgdnMgY2x1bXAgYXJlYQoKYGBge3J9CnNjYXR0ZXJwbG90IChpbmRpdn5hcmVhLCBkYXRhPXBlYWtxdWlubikKcHEyLmxtPC1sbShpbmRpdn5hcmVhLGRhdGE9cGVha3F1aW5uKQpwbG90KHBxMi5sbSkKYGBgCgpTY2F0dGVycGxvdCBzdWdnZXN0cyBub24tbGluZWFyIHJlbGF0aW9uc2hpcCBhbmQgdGhlIHJlc2lkdWFsIHBsb3Qgc2hvd3MgYSBwYXR0ZXJuLiBVc2UgbG9nLXRyYW5zZm9ybWVkIGRhdGEKCmBgYHtyfQpwcTMubG08LWxtKGxpbmRpdn5sYXJlYSwgZGF0YT1wZWFrcXVpbm4pCnBsb3QocHEzLmxtKQphdWdtZW50KHBxMy5sbSkKZ2xhbmNlKHBxMy5sbSkKdGlkeShwcTMubG0sIGNvbmYuaW50ID0gVFJVRSwgZGlnaXRzPTMpCiNBbmQgaWYgeW91IHdhbnQgdGhlIHRyYWRpdGlvbmFsIEFOT1ZBIHRhYmxlLi4uCmFub3ZhKHBxMy5sbSkKCmBgYAoKIyMgTmljZXIgZ3JhcGhpYwoKIyMjIEZpZyA2LjEzCgpGaXJzdCBsb2FkIHNjcmlwdCB0aGF0IHByb2R1Y2VzIHN0YW5kYXJkIGdyYXBoaWNzIGFwcGVhcmFuY2UgKGNvZGUgaGlkZGVuKQoKYGBge3IgaW5jbHVkZT1GQUxTRSwgcmVzdWx0cz0naGlkZSd9CnNvdXJjZSgiLi4vUi9hcHBlYXJhbmNlLlIiKSAjVGhpcyBpcyB0aGUgY29tbW9uIGxpYnJhcnkgCmBgYAoKYGBge3J9CnA8LWdncGxvdChwZWFrcXVpbm4sIGFlcyh4PWxhcmVhLCB5PWxzcGVjaWVzKSkrCiAgZ2VvbV9wb2ludChjb2xvcj1zYykrCiAgZ2VvbV9zbW9vdGgobWV0aG9kPSJsbSIsIHNlPUZBTFNFLCBjb2xvcj1sYykrCiAgdGhlbWVfY2xhc3NpYyhiYXNlX3NpemUgPSAxMCkrCiAgdGhlbWUoCiAgICBheGlzLnRleHQgPSBlbGVtZW50X3RleHQoY29sb3VyID0gYWMpLAogICAgYXhpcy5saW5lID0gZWxlbWVudF9saW5lKGNvbG9yID0gYWMpLAogICAgYXhpcy50aWNrcyA9IGVsZW1lbnRfbGluZShjb2xvciA9IGFjKSwKICAgICAgICApKwogIHhsYWIoIkxvZyBhcmVhIikreWxhYigiTG9nIG5vLiBzcGVjaWVzIikKI0FkZCBtYXJnaW5hbCBib3ggcGxvdHMKcDE8LWdnTWFyZ2luYWwocCxzaXplPTIwLCB0eXBlPSJib3hwbG90IikKcDEKYGBgCgpSZXNpZHVhbCBwbG90CgpgYGB7cn0KcDI8LWdncGxvdChwcS5sbSwgYWVzKHggPSBwcS5sbSRmaXR0ZWQudmFsdWVzLCB5ID0gcHEubG0kcmVzaWR1YWxzKSkgKwogIGdlb21fcG9pbnQoY29sb3I9c2MpICsKICBnZW9tX3Ntb290aChzZT1GQUxTRSwgY29sb3VyPWxjKSsKICB0aGVtZV9jbGFzc2ljKGJhc2Vfc2l6ZSA9IDEwKSsKICB0aGVtZSgKICAgIGF4aXMudGV4dCA9IGVsZW1lbnRfdGV4dChjb2xvdXIgPSBhYyksCiAgICBheGlzLmxpbmUgPSBlbGVtZW50X2xpbmUoY29sb3IgPSBhYyksCiAgICBheGlzLnRpY2tzID0gZWxlbWVudF9saW5lKGNvbG9yID0gYWMpLAogICAgICAgICkrbGFicyh4ID0gIlByZWRpY3RlZCBsb2cgbm8uIHNwZWNpZXMiLCB5ID0gIlJlc2lkdWFscyIsIAogICAgICAgKQpwMgpgYGAKCiMjIyMgRmlnIDYuMTQgLSBpbmRpdiB2cyBhcmVhCgpgYGB7cn0KcDM8LWdncGxvdChwZWFrcXVpbm4sIGFlcyh4PWFyZWEsIHk9aW5kaXYpKSsKICBnZW9tX3BvaW50KGNvbG9yID0gc2MpKwogIGdlb21fc21vb3RoKHNlID0gRkFMU0UsIGNvbG9yID0gbGMpKwogIHRoZW1lX2NsYXNzaWMoYmFzZV9zaXplID0gMTApKwogIHRoZW1lKAogICAgYXhpcy50ZXh0ID0gZWxlbWVudF90ZXh0KGNvbG91ciA9IGFjKSwKICAgIGF4aXMubGluZSA9IGVsZW1lbnRfbGluZShjb2xvciA9IGFjKSwKICAgIGF4aXMudGlja3MgPSBlbGVtZW50X2xpbmUoY29sb3IgPSBhYyksCiAgICAgICAgKSsKICB4bGFiKCJBcmVhIikreWxhYigiTm8uIGluZGl2aWR1YWxzIikKI0FkZCBtYXJnaW5hbCBib3ggcGxvdHMKcDQ8LWdnTWFyZ2luYWwocDMsc2l6ZT0yMCwgdHlwZT0iYm94cGxvdCIpCnA0CmBgYAoKUmVzaWR1YWwgcGxvdAoKYGBge3J9CnA1PC1nZ3Bsb3QocHEyLmxtLCBhZXMoeCA9IHBxMi5sbSRmaXR0ZWQudmFsdWVzLCB5ID0gcHEyLmxtJHJlc2lkdWFscykpICsKICBnZW9tX3BvaW50KGNvbG9yPXNjKSArCiAgdGhlbWVfY2xhc3NpYyhiYXNlX3NpemUgPSAxMCkrCiAgdGhlbWUoCiAgICBheGlzLnRleHQgPSBlbGVtZW50X3RleHQoY29sb3VyID0gYWMpLAogICAgYXhpcy5saW5lID0gZWxlbWVudF9saW5lKGNvbG9yID0gYWMpLAogICAgYXhpcy50aWNrcyA9IGVsZW1lbnRfbGluZShjb2xvciA9IGFjKSwKICAgICAgICApK2xhYnMoeCA9ICJQcmVkaWN0ZWQgbG9nIG5vLiBzcGVjaWVzIiwgeSA9ICJSZXNpZHVhbHMiLCAKICAgICAgICkKcDUKYGBgCgojIyMjIEZpZyA2LjE1CgpgYGB7cn0KcDY8LWdncGxvdChwZWFrcXVpbm4sIGFlcyh4PWxhcmVhLCB5PWxpbmRpdikpKwogIGdlb21fcG9pbnQoY29sb3I9c2MpKwogIGdlb21fc21vb3RoKG1ldGhvZD0ibG0iLCBjb2xvcj0iYmxhY2siKSsKICB0aGVtZV9jbGFzc2ljKGJhc2Vfc2l6ZSA9IDEwKSsKICB0aGVtZSgKICAgIGF4aXMudGV4dCA9IGVsZW1lbnRfdGV4dChjb2xvdXIgPSBhYyksCiAgICBheGlzLmxpbmUgPSBlbGVtZW50X2xpbmUoY29sb3IgPSBhYyksCiAgICBheGlzLnRpY2tzID0gZWxlbWVudF9saW5lKGNvbG9yID0gYWMpLAogICAgICAgICkrCiAgeGxhYigiTG9nIGFyZWEiKSt5bGFiKCJMb2cgbm8uIGluZGl2LiIpCiNBZGQgbWFyZ2luYWwgYm94IHBsb3RzCnA3PC1nZ01hcmdpbmFsKHA2LHNpemU9MjAsIHR5cGU9ImJveHBsb3QiKQpwNwpgYGAKClJlc2lkdWFsIHBsb3QKCmBgYHtyfQpwODwtZ2dwbG90KHBxMy5sbSwgYWVzKHggPSBwcTMubG0kZml0dGVkLnZhbHVlcywgeSA9IHBxMy5sbSRyZXNpZHVhbHMpKSArCiAgZ2VvbV9wb2ludChjb2xvcj1zYykgKwogIHRoZW1lX2NsYXNzaWMoYmFzZV9zaXplID0gMTApKwogIHRoZW1lKAogICAgYXhpcy50ZXh0ID0gZWxlbWVudF90ZXh0KGNvbG91ciA9IGFjKSwKICAgIGF4aXMubGluZSA9IGVsZW1lbnRfbGluZShjb2xvciA9IGFjKSwKICAgIGF4aXMudGlja3MgPSBlbGVtZW50X2xpbmUoY29sb3IgPSBhYyksCiAgICAgICAgKStsYWJzKHggPSAiUHJlZGljdGVkIGxvZyBuby4gaW5kaXYuIiwgeSA9ICJSZXNpZHVhbHMiLCAKICAgICAgICkKcDgKYGBgCg==