QK Box 8.8

Paruelo and Lauenroth (1996) analyzed the geographic distribution and the effects of climate variables on the relative abundance of a number of plant functional types (PFTs) including shrubs, forbs, succulents (e.g. cacti), C3 grasses and C4 grasses. There were 73 sites across North America. The response variable we will focus on is the relative abundance of C3 plants and there were six potential predictors: the latitude in centesimal degrees (LAT), the longitude in centesimal degrees (LONG), the mean annual precipitation in mm (MAP), the mean annual temperature in °C (MAT), the proportion of MAP that fell in June, July and August (JJAMAP) and the proportion of MAP that fell in December, January and February (DJFMAP).

This is an example from the first edition; the data file is here

Paruelo, J. M. & Lauenroth, W. K. (1996). Relative abundance of plant functional types in grasslands and shrublands of North America. Ecological Applications, 6, 1212-24.

Preliminaries

First, load the required packages (car, lm.beta, reghelper, Rmisc and rgl)

Import paruelo data file (paruelo.csv)

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

Diagnostic checks

We will use the scatterplotMatrix function from the car package

scatterplotMatrix(~c3+lat+long+map+mat+jjamap+djfmap, data=paruelo, cex=0.33, regLine=FALSE, diagonal=list(method='boxplot'))

C3 positively skewed but not too bad; no obvious non-linearities

Now check for collinearity among all predictors - correlations and VIFs

cor(paruelo[,c('lat','long','map','mat','jjamap','djfmap')])

           lat    long    map      mat  jjamap   djfmap
lat     1.0000  0.0966 -0.247 -0.83859  0.0742 -0.06512
long    0.0966  1.0000 -0.734 -0.21311 -0.4916  0.77074
map    -0.2465 -0.7337  1.000  0.35509  0.1123 -0.40451
mat    -0.8386 -0.2131  0.355  1.00000 -0.0808  0.00148
jjamap  0.0742 -0.4916  0.112 -0.08077  1.0000 -0.79154
djfmap -0.0651  0.7707 -0.405  0.00148 -0.7915  1.00000

vif(lm(c3~lat+long+map+mat+jjamap+djfmap, data=paruelo))

   lat   long    map    mat jjamap djfmap 
  3.50   5.27   2.80   3.74   3.16   5.71 

Obvious collinearity with environmental variables tending to be correlated strongly with lat or long, so we’ll split predictors into 2 groups, one with lat and long

First model c3 against lat and long and lat*long and check for collinearity

vif(lm(c3~lat+long+lat*long, data=paruelo))

there are higher-order terms (interactions) in this model
consider setting type = 'predictor'; see ?vif

     lat     long lat:long 
   307.7     66.8    400.9 

VIFs very high indicating obvious collinearity as expected with interaction term

Run uncentred model

paruelo.lm <- lm(c3~lat+long+lat*long, data=paruelo)

Get regression diagnostics (residual plot and Cooks D)

plot(paruelo.lm)

augment(paruelo.lm)

Display results of model fitting

glance(paruelo.lm)

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

Get standardized regression coefficients, use lm.beta function from lm.beta package

lm.beta.paruelo <- lm.beta(paruelo.lm)
tidy(lm.beta.paruelo, conf.int=TRUE)

Centre predictors and refit model

paruelo$clat <- scale(paruelo$lat, center=TRUE, scale=FALSE)
paruelo$clong <- scale(paruelo$long, center=TRUE, scale=FALSE)
paruelo1.lm <- lm(c3~clat+clong+clat*clong, data=paruelo)

Recheck collinearity

vif(lm(c3~clat+clong+clat*clong, data=paruelo))

there are higher-order terms (interactions) in this model
consider setting type = 'predictor'; see ?vif

      clat      clong clat:clong 
      1.21       1.02       1.22 

Regression diagnostics (residual plot and Cooks D)

plot(paruelo1.lm)

augment(paruelo1.lm)

Diagnostics look better now

Display results of model fitting

glance(paruelo1.lm)

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

Get standardized regression coefficients, use lm.beta function from lm.beta package

lm.beta(paruelo1.lm)

Call:
lm(formula = c3 ~ clat + clong + clat * clong, data = paruelo)

Standardized Coefficients::
(Intercept)        clat       clong  clat:clong 
         NA    0.770191    0.000703    0.249471 

Standardise predictors and refit model

paruelo$slat <- scale(paruelo$lat, center=TRUE, scale=TRUE)
paruelo$slong <- scale(paruelo$long, center=TRUE, scale=TRUE)
paruelo2.lm <- lm(c3~slat+slong+slat*slong, data=paruelo)

Recheck collinearity

vif(lm(c3~slat+slong+slat*slong, data=paruelo))

there are higher-order terms (interactions) in this model
consider setting type = 'predictor'; see ?vif

      slat      slong slat:slong 
      1.21       1.02       1.22 

Regression diagnostics (residual plot and Cooks D)

plot(paruelo2.lm)

augment(paruelo2.lm)

Display results of model fitting

glance(paruelo2.lm)

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

Get standardised regression coefficients, use lm.beta function from lm.beta package

lm.beta(paruelo2.lm)

Call:
lm(formula = c3 ~ slat + slong + slat * slong, data = paruelo)

Standardized Coefficients::
(Intercept)        slat       slong  slat:slong 
         NA    0.770191    0.000703    0.249471 

3d scatterplot with smoothing surface

Note on mac OS - need xquartz installed. Graphic opens in separate window in Mac/Xquartz; haven’t checked on other operating systems.

Rotate the graph using the mouse or cursor keys

scatter3d(c3~lat+long, grid = FALSE, fit = "smooth", data=paruelo)

Examine simple slopes using reghelper package

simple_slopes(paruelo.lm)

        lat       long Test Estimate Std. Error t value df Pr(>|t|) Sig.
1 34.801242     sstest        -0.010      0.005  -2.039 69   0.0453    *
2 40.104247     sstest         0.000      0.004   0.008 69   0.9935     
3 45.407251     sstest         0.010      0.005   1.848 69   0.0689    .
4    sstest   99.96472         0.026      0.005   5.159 69 2.27e-06  ***
5    sstest 106.400137         0.038      0.005   8.198 69 8.69e-12  ***
6    sstest 112.835554         0.050      0.008   6.484 69 1.13e-08  ***