QK Box 6.2

Christensen et al. (1996) studied the relationships between coarse woody debris (CWD) and shoreline vegetation and lake development in a sample of 16 lakes in North America. The main variables of interest are the density of cabins (no. km⁻¹), density of riparian trees (trees km⁻¹), the basal area of riparian trees (m² km⁻¹), density of coarse woody debris (no. km⁻¹), basal area of coarse woody debris (m² km⁻¹).

Coarse woody debris provides important habitat structure. Not quite a lake in Wisconsin, but World Heritage rainforest in northern Australia! Photo M. Keough

Coarse woody debris provides important habitat structure. Not quite a lake in Wisconsin, but World Heritage rainforest in northern Australia! Photo M. Keough

Christensen, D. L., Herwig, B. R., Schindler, D. E. & Carpenter, S. R. (1996). Impacts of lakeshore residential development on coarse woody debris in north temperate lakes. Ecological Applications, 6, 1143-49.

The original paper: https://doi.org/10.2307/2269598

This data set was used for the first edition, and is available here christ.csv

Preliminaries

First, load the required packages (car, lm.beta).

Import Christensen data (christ.csv)

christ <- read.csv("../data/christ.csv")
head(christ,10)

NA

CWD versus riparian density first

scatterplot(cwdbasal~ripdens, data=christ)

christ.lm <- lm(cwdbasal ~ ripdens, data=christ)

Check diagnostics

plot(christ.lm)

augment(christ.lm)

Examine model output

glance(christ.lm)

tidy(christ.lm, conf.int = TRUE)

anova(christ.lm)

Analysis of Variance Table

Response: cwdbasal
          Df Sum Sq Mean Sq F value    Pr(>F)    
ripdens    1  32054   32054  24.303 0.0002216 ***
Residuals 14  18466    1319                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get standardized coefficients

lm.beta(christ.lm)

Call:
lm(formula = cwdbasal ~ ripdens, data = christ)

Standardized Coefficients::
(Intercept)     ripdens 
         NA   0.7965489 

Predict new CWD values for riparian densities of 1000 and 2000 trees per km

predict(christ.lm, data.frame(ripdens=c(1500)), interval="prediction", se=T)

$fit
       fit      lwr     upr
1 96.17503 14.88709 177.463

$se.fit
[1] 10.83789

$df
[1] 14

$residual.scale
[1] 36.31761

Now CWD versus cabin density

scatterplot(cwdbasal~cabin, data=christ)

christ1.lm <- lm(cwdbasal ~ cabin, data=christ)

Check diagnostics

plot(christ1.lm)

augment(christ1.lm)

anova(christ1.lm)

Analysis of Variance Table

Response: cwdbasal
          Df Sum Sq Mean Sq F value   Pr(>F)   
cabin      1  25440 25440.1  14.201 0.002076 **
Residuals 14  25080  1791.4                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Examine model output

Transform cabin density to log10 (cabin density + 1) and refit model

christ$lcabin <- log10(christ$cabin+1)
christ2.lm <- lm(cwdbasal ~ lcabin, data=christ)

Get diagnostics

plot(christ2.lm)

augment(christ2.lm)

Examine model output

glance(christ2.lm)

tidy(christ2.lm, conf.int = TRUE)

anova(christ2.lm)

Analysis of Variance Table

Response: cwdbasal
          Df Sum Sq Mean Sq F value    Pr(>F)    
lcabin     1  32840   32840  26.004 0.0001619 ***
Residuals 14  17680    1263                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get standardized coefficients

lm.beta(christ2.lm)

Call:
lm(formula = cwdbasal ~ lcabin, data = christ)

Standardized Coefficients::
(Intercept)      lcabin 
         NA  -0.8062482 

Fit log cwd vs log cabin model to get residual plots, etc. for figure used in book

christ$lcwd<-log10(christ$cwdbasal)
christ3.lm<-lm(christ$lcwd~christ$lcabin)
glance(christ3.lm)

tidy(christ3.lm)

anova(christ3.lm)

Analysis of Variance Table

Response: christ$lcwd
              Df Sum Sq Mean Sq F value    Pr(>F)    
christ$lcabin  1 5.7336  5.7336  26.141 0.0001579 ***
Residuals     14 3.0707  0.2193                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1