generated from rstudio/bookdown-demo
-
Notifications
You must be signed in to change notification settings - Fork 3
/
07-habitat_use.Rmd
590 lines (377 loc) · 25.5 KB
/
07-habitat_use.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
---
output: html_document
editor_options:
chunk_output_type: console
---
# Habitat use {#habitat-use}
Camera traps are well suited for the quantification of habitat use across multiple species. To assess habitat use, we typically quantify the detection rate - the number of detections divided by the time interval of interest. As detection rates are fairly simple to estimate and conceptually simple to understand, thus their use is widespread in the camera trap literature.
In its simplest form habitat use represents the number of independent events of a given species at a given camera, divided by the number of days that camera was active during that period of interest.
This 'detection rate' is thought to reflect the habitat use of a species at a given location. Extreme care should be taken if you want to equate use with abundance or density - something we discuss a little more in the [density chapter](#density).
Detection rates are typically analysed in a linear modelling framework, and come in single species and multi-species versions (see below).
*Create a new .R script*
Call it `05_example_habitat_use.R`.
*Load the required packages*
```{r ch7_1, echo=T, results='hide', message =F, warning=F, class.source="Rmain"}
# Check you have them and load them
list.of.packages <- c("kableExtra", "tidyr", "ggplot2", "gridExtra", "lme4", "dplyr", "Hmsc", "jtools", "lubridate", "corrplot", "MuMIn")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
lapply(list.of.packages, require, character.only = TRUE)
```
## Calculating capture rate
We will start by using the `total_obs` dataframe we have used in previous chapters:
```{r ch7_2, class.source="Rmain"}
# Import the total observations dataset
total_obs <- read.csv("data/processed_data/AlgarRestorationProject_30min_independent_total_observations.csv", header=T)
# Import your species list
sp_summary <- read.csv("data/processed_data/AlgarRestorationProject_species_list.csv", header=T)
```
Which, as a quick reminder, looks like this:
```{r ch7_3, echo=F}
kbl(head(total_obs))%>%
kable_paper() %>%
scroll_box(height = "200px")
```
So within each row, we have the location (`placename`), the survey effort at each given location (camera `days`), and the number of independent records of each species. This is very close to the format which most linear model analysis packages require. Easy!
Our next step to create the capture rate - our proxy for habitat use. We will divide each count by the number of days cameras were active in that location, then multiply by the a standardized number of days - often people use 100.
In R this would look like:
```{r ch7_4, class.source="Rmain"}
# Create a dataframe to store these detection rates
total_cr <- total_obs
# Divide the species abundances (which start in column four), by the amount of camera effort
total_cr[ ,sp_summary$sp ] <- (total_cr[ , sp_summary$sp]/total_cr$days)*100
```
We can then examine the relationship between raw counts (on the x-axis) with our detection rate (on the y-axis), using *Odocoileus virginianus* as an example. In the plot below each black dot represents a `placename` where camera trapping has occurred.
```{r ch7_5, class.source="Rinfo"}
plot(total_cr$Odocoileus.virginianus ~ total_obs$Odocoileus.virginianus,
las=1, pch=19,
ylab="Capture rate per 100 days",
xlab="Number of independent records")
```
As you can see they are not a perfect match as the capture rate accounts for the variation in effort between different sites.
### Examples from the literature
[Palmer, Meredith S., et al. "Evaluating relative abundance indices for terrestrial herbivores from large‐scale camera trap surveys." African journal of ecology 56.4 (2018): 791-803.](https://onlinelibrary.wiley.com/doi/abs/10.1111/aje.12566)
## Single-species models
The most common way to analyse habitat-use data is through linear models. Linear models typically relate a continuous response variable - in our case capture rate - to a set of one or more discrete or continuous predictor variables. In this simple example we will explore the relationship between the capture rate of a species with the categorical 'feature_type' variable and the continuous `line_of_sight_m` variables.
There are a variety if different frameworks to fit and compare different linear models to address a host of different hypotheses, but if you are just starting out you should investigate two widely used packages:
- `lme4` -> frequentest and information theoretic approaches
- `brms` -> Bayesian approaches
There is no right or wrong about which package and which approach you use to test your hypotheses. Some packages have functionalities that others don't, which may force your hand. Just make sure you understand the implications of your choices when it comes to reporting your results!
### Simple linear model
We will start by analyzing a frequentest linear model with a single observation for each camera location.
In this worked example we will analyse how habitat use varies using a linear model `lm()`. The model takes the form:
Response term (y) ~ fixed effect 1 (x1) + fixed effect 2 (x2), data frame (data=)
It is beyond the scope of this course to test the model assumptions or interrogate the findings, there are better resources to allow you to do that (e.g. we highly recommend reading [Gałecki, Andrzej, and Tomasz Burzykowski. "Linear mixed-effects model." Linear mixed-effects models using R. Springer, New York, NY, 2013. 245-273](https://link.springer.com/book/10.1007/978-1-4614-3900-4)).
In this example we will explore if the habitat use of *Odocoileus virginianus* varies based on the `feature_type' the line of sight where it is found.
**Preparing our data**
Recall that the information about each location is recorded in the file:
```{r ch7_6, class.source="Rmain"}
locs <- read.csv("data/processed_data/AlgarRestorationProject_camera_locations_and_covariates.csv", header=T)
# Convert to categorical factors
locs <- locs %>%
mutate_if(is.character,as.factor)
# You should also standardize your covariates - it helps models coverage an facillitates comparison of effects sizes
library(MuMIn)
z_locs <- stdize(locs)
```
Take a look at it to see what it has done!
```{r ch7_6a, echo=F}
kbl(head(z_locs))%>%
kable_paper() %>%
scroll_box(width = "750px", height = "200px")
```
Each location (`placename`) has a single row in the dataframe.
We will now add our covariates to our capture rate dataframe:
```{r ch7_7, class.source="Rmain", message=F, warning=F}
mod_dat <- left_join(total_cr, z_locs) # from the dplyr package
```
### Catagorical predictor
So we start by exploring the influence of 'feature_type' on our response term.
What are the categories we have in our feature type variable?
```{r ch7_8, class.source="Rinfo"}
table(z_locs$feature_type)
```
`feature_type` is a a categorical variable which reflects strata where the camera trap was deployed:
- HumanUse = a camera on a seismic line used and maintained in an "open" state by humans
- Offline = a camera in contiguous forest >200m from a seismic line
- NatRegen = a seismic line which is naturally regenerating
Lets do a quick raw data plot to see what results we might expect:
```{r ch7_9, class.source="Rinfo"}
boxplot(mod_dat$Odocoileus.virginianus~mod_dat$feature_type,
las=1,
xlab="feature_type",
ylab="Habitat use")
```
It looks like white-tailed deer habitat use may be higher in naturally regenerating areas, but there is a lot of overlap between sites.
Next we will fit a simple linear model using the `lm()' function in base R.
```{r ch7_10, class.source="Rmain"}
# model results <- lm( Y data ~ x Data, data= dataframe source)
lm_cat <- lm(Odocoileus.virginianus ~ feature_type, data = mod_dat)
```
Lets looks at the model summary:
```{r ch7_11, class.source="Rmain"}
summary(lm_cat)
```
Categorical covariates are show as contrasts from the reference level (in this case HumanUse), and the p-value relate to testing whether the other categories are significantly different from the reference level. Other things to note are that our R-squared value (how much variation the model explains) is fairly low - but that is common in camera trap models.
We can take a quick look at the predictions using the `jtools` package. More examples of its use are can be found in the `[Visualizing regression model predictions vignette](https://jtools.jacob-long.com/articles/effect_plot.html) associated with the package.
```{r ch7_12, class.source="Rmain"}
effect_plot(lm_cat, # The model object
pred = feature_type, # The variable you want to predict
interval = TRUE, # Whether you want confidence intervals (default = 0.95)
partial.residuals = T, # Show the residual variation -after accounting for fixed effects
y.label = "Habitat use") # Change the y axis label
```
### Continuous predictor
Let's also explore a continuous predictor `line_of_sight_m':
```{r ch7_13, class.source="Rinfo"}
plot(mod_dat$Odocoileus.virginianus~mod_dat$z.line_of_sight_m,
las=1,
xlab="line_of_sight_m",
ylab="Habitat use")
```
Next we will fit a simple linear model using the `lm()' function in base R.
```{r ch7_14, class.source="Rmain"}
# model results <- lm( Y data ~ x Data, data= dataframe source)
lm_con <- lm(Odocoileus.virginianus ~ z.line_of_sight_m, data = mod_dat)
```
Lets looks at the model summary:
```{r ch7_15, class.source="Rmain"}
summary(lm_con)
```
Here the effect represents the gradient of the relationship between `line_of_sight_m` and the habitat use of white-tailed deer. The effect is negative, and the p-value is below the arbitrary 0.05 threshold, which suggests it my be an important predictor of white-tailed deer habitat use.
It will make more sense if we plot it - again using `jtools`
```{r ch7_16, class.source="Rmain"}
effect_plot(lm_con, # The model object
pred = z.line_of_sight_m, # The variable you want to predict
interval = TRUE, # Whether you want confidence intervals (default = 0.95)
partial.residuals = T, # Show the residual variation -after accounting for fixed effects
y.label = "Habitat use") # Change the y axis label
```
### Model comparisons
There are times when we may want to compare which model is "the best", or which model is the most parsimonious. One way to do this is through the use of Information Theory - we can compare which model explains the most amount of variation after applying a penalty for how complex it is (more complex models will always explain more variation, even if just by chance).
One useful package for this is the `MuMIn` package and the function `model.sel()` for model selection:
```{r ch7_17, class.source="Rmain"}
library(MuMIn)
# Lets also create a "null model" something without any predictors in at all, to compare these models to:
lm_null <- lm(Odocoileus.virginianus ~ 1, data = mod_dat)
```
Now compare the three alternative models:
```{r ch7_18, class.source="Rmain"}
model.sel(lm_null, lm_cat, lm_con)
```
Whilst both models improve on the null model, there is stronger support for `line_of_sight_m` than for our feature types in influencing white-tailed deer habitat use. Cool!
### Problems with these models
But can you see any problems with this type of model?
We probably should be concerned about the fact that:
- There are negative predictions for both sets of confidence intervals - but you can't get a negative capture rate!
- We do not account for seasonality - we saw species detection rates change with time of year in the [data exploration section](#exploration)
And more besides!
### Mixed-effects models
Let's build a more robust habitat-use model which addresses some of the issues highlighted here. To do this we will take advantage of a type of analysis called 'mixed effects modelling'. Mixed effects models allow us to perform robust analysis of populations which have been repeatedly sampled through time. As such, we can break our data set down into months without violating the assumptions of the models.
If you are new to mixed effects models you must try this fantastic interactive aid to help you understand how they work: [Michael Freeman's 'An Introduction to Hierarchical Modeling'](http://mfviz.com/hierarchical-models/)
And for a deep-dive into the inner workings of mixed effects models and their assumptions, see the following paper:
[Harrison, Xavier A., et al. "A brief introduction to mixed effects modelling and multi-model inference in ecology." PeerJ 6 (2018): e4794.](https://peerj.com/articles/4794/)
First we must install the packages we require: 'lme4' and `tidyr':
```{r ch7_19, class.source="Rmain"}
library(lme4); library(tidyr)
```
The lme4 package requires a dataframe format (as above), with the response term and the predictor variables all included in the same location.
Second, lets create our monthly analysis dataframe:
```{r ch7_20, class.source="Rmain"}
# Import the total observations dataset
monthly_obs <- read.csv("data/processed_data/AlgarRestorationProject_30min_independent_monthly_observations.csv", header=T)
```
Let's join this dataframe with the location data, as before:
```{r ch7_21, class.source="Rmain", message=F, warning=F}
mod_dat <- left_join(monthly_obs, z_locs)
```
And let's do another raw data check:
```{r ch7_22, class.source="Rinfo"}
boxplot(mod_dat$Odocoileus.virginianus~mod_dat$feature_type,
las=1,
xlab="feature_type",
ylab="Habitat use")
```
The patterns are similar to the `total_obs` dataset, however there is more noise!
Now that we have monthly data, we might also want to control for some element of seasonality in our models. We can extract the month from our date column.
```{r ch7_23, class.source="Rmain"}
mod_dat$date <- ym(mod_dat$date)
mod_dat$month<- month(mod_dat$date, label=T)
```
Next we will fit a mixed effects model to this data set using `lme4`. You may have noticed that we haven't calculated a separate capture rate dataframe as we did in the simple example! That is because we can create a relative abundance index within the model itself by providing an `offset()` term . An offset term serves to scale the response term based on the amount of survey effort, and preserves the original units of the observations (counts).
The model takes the form:
Response term ~ fixed effect + offset() + (1|random intercept), data frame, distribution
We include `placename` as the random intercept, as camera locations are repeatedly sampled at monthly intervals and thus our data (rows in the dataframe) are not independent. We use the `poisson` family, as our response term is a count.
```{r ch7_24,eval=F, class.source="Rmain"}
glmm_cat <- glmer(Odocoileus.virginianus ~
feature_type + month + offset(log(days)) + (1|placename) , data=mod_dat, family="poisson")
```
Oh that warning doesn't look friendly. Convergence errors often arise as when we fit a model which is too complicated! It basically says the the model hasn't completely "solved" the parameters within it. In it's current form the model is estimating and effect for every single month of the year, that's a lot of thing to estimate with such a small dataset.
Let's simplify things to `summer` and `winter`!
```{r ch7_25, class.source="Rmain"}
# Lets create a new column and give it the value summer
mod_dat$season <- "summer"
mod_dat$season[month(mod_dat$date) %in% c(10,11,12,1,2,3)] <- "winter"
# make it a factor factors
mod_dat <- mod_dat %>%
mutate_if(is.character,as.factor)
```
And re-run our model - also this time with a negative binomial distribution to account for excess zeros:
```{r ch7_26, class.source="Rmain"}
glmm_cat <- glmer.nb(Odocoileus.virginianus ~
feature_type + season + offset(log(days)) + (1|placename) , data=mod_dat)
```
We can view a summary of the model fit using:
```{r ch7_27, class.source="Rmain"}
summary(glmm_cat)
```
We can plot the predictions from these models using the `jtools` package.
First lets look at the effects of `feature_type`:
```{r ch7_28, class.source="Rmain"}
effect_plot(glmm_cat, pred = feature_type, interval = TRUE, y.label = "Habitat use",
, data=mod_dat)
```
As with our simple linear model, the mixed effects model also suggests a difference between the different `feature_type` strata for white-tailed deer.
Lets also look at the effect of our new `season` variable:
```{r ch7_29, class.source="Rmain", message=F, warning=F}
effect_plot(glmm_cat, pred = season, interval = TRUE, y.label = "Habitat use",
, data=mod_dat)
```
Which suggests habitat use is slightly lower habitat use in winter then in summer.
### On your own
Using this code as a scaffold, explore some more of the patterns we explored in the [data exploration chapter](#exploration).
Remember we have the following species:
```{r ch7_30, echo=F}
sp_summary$sp
```
And the following covariates:
```{r ch7_31, echo=F}
colnames(z_locs)[colnames(z_locs) %in% c("feature_type","z.line_of_sight_m", "z.water_depth_m", "z.elevation", "z.road_dist_m", "z.water_dist_m", "z.mean_ndvi")]
```
**A note of caution**
As stated at the start of this guide, we are not focusing on whether the models we apply are appropriate or finding "the best" models for this datasheet, so do not spend too much time trying to interpret this information!
### Advanced mixed-model predictions
Tools such as `jtools` are great for generating simple predictions from mixed models, however the more complex the models get, the more you may want to specify your own prediction dataframes.
If you want more applied examples of generating predictions from mixed effects models, check out [Ben Bolkers workbook](https://bbolker.github.io/mixedmodels-misc/ecostats_chap.html). There is also some great discussion about model selection and r-squared values.
### Examples in the literature
[Tattersall, E. R., Burgar, J. M., Fisher, J. T., & Burton, A. C. (2020). Mammal seismic line use varies with restoration: Applying habitat restoration to species at risk conservation in a working landscape. Biological Conservation, 241, 108295.](https://www.sciencedirect.com/science/article/abs/pii/S0006320719307013)
## Multispecies models
In the above examples, we analyse each individual species separately. This is great if you only care about one species, however we often want a more holistic understanding of wildlife communities! Recent advances in computer power and analytic approaches mean it is becoming increasingly popular to model multiple species within the same framework! This opens up a variety of things not previously possible.
**A note of caution** In experimenting with single species models you may have realized it can sometimes be hard to build a sensible and robust model. Now do this for >10 species in the **same model**, and the potential to get silly results increases. Tread carefully!
As with single species linear models, there are many choices available for modeling multiple species in the same framework. Two notable options are:
- [GJAM](https://cran.r-project.org/web/packages/gjam/vignettes/gjamVignette.html)
- [HMSc](https://www2.helsinki.fi/en/researchgroups/statistical-ecology/hmsc)
In this example we will use the `Hmsc` package.
```{r ch7_32, class.source="Rmain"}
library(Hmsc)
```
**Preparing our data**
The format of data required for joint species distribution models is very similar to the data required for single species models. However, rather than storing the response term and fixed effects within the the same data frame (as with `mod_dat` above), we need a separate `Y` matrix of site_time x species, and a separate `Xdata` dataframe containing the fixed and random effects.
```{r ch7_33, class.source="Rmain"}
# Pull the count data into its own matrix
Y <- as.matrix(monthly_obs[,sp_summary$sp])
# Give the row names a useful label, in this case the site_date values
# (just in case you want to check things)
row.names(Y) <- paste(monthly_obs$placename, monthly_obs$date, sep="_")
```
Which looks like this:
```{r ch7_34, echo=F}
kbl(head(Y))%>%
kable_paper() %>%
scroll_box(height = "200px")
```
We then create the XData in a similar way to before, but this time dropping the species information:
```{r ch7_35, class.source="Rmain", message=F, warning=F}
Xdat <- left_join(monthly_obs[c("placename", "date", "days")], z_locs)
# All XData must be numeric or factors, so lets check what we have
```
Whose output looks like this:
```{r ch7_36, echo=F}
kbl(head(Xdat))%>%
kable_paper() %>%
scroll_box(height = "200px")
```
With Bayesian approaches we need to set up our sampling conditions
```{r ch7_37, class.source="Rmain"}
nChains = 2 # How many total repeats to run
thin = 5 # How often to thin the samples
samples = 100 # How many samples to take
transient = 1000 # How long should the "warm up" be
verbose = T # Give reports on model progress
```
Setup our random effect:
```{r ch7_38, class.source="Rmain"}
# Add a station-level random effect (for the co-variances)
studyDesign = data.frame(station = as.factor(Xdat$placename))
rL = HmscRandomLevel(units = studyDesign$station)
```
Specify our model""
```{r ch7_39, message =F, warning=F, eval=T, class.source="Rmain"}
# Model specification
mod <- Hmsc(Y = Y,
XData = Xdat[,c("z.line_of_sight_m", "z.water_depth_m", "days")],
XFormula = ~z.line_of_sight_m + z.water_depth_m + log(days),
studyDesign = studyDesign,
ranLevels = list(station = rL),
distr="poisson")
```
And fit the model:
```{r ch7_39a, eval=T,results="hide", class.source="Rmain"}
out <- sampleMcmc(mod, thin = thin, samples = samples, transient = transient,
nChains = nChains, verbose = verbose)
```
We can plot a basic summary of the modeled effects using the following code.
```{r ch7_40, class.source="Rmain"}
postBeta = getPostEstimate(out, parName = "Beta")
par(mar=c(8,12,1,1))
plotBeta(out, post = postBeta, param = "Support", supportLevel = 0.95)
```
We the colors denote the size and magnitude of the effect of proportion of lowland habitat. *NOTE* treat these results with caution as the number of model runs is very low (to increase speed) and the model assumptions have not been interrogated.
```{r ch7_41, class.source="Rmain"}
OmegaCor = computeAssociations(out)
supportLevel = 0.0
toPlot = ((OmegaCor[[1]]$support>supportLevel)
+ (OmegaCor[[1]]$support<(1-supportLevel))>0)*OmegaCor[[1]]$mean
corrplot(toPlot,
method = "color",
type="upper",
order = "FPC",
col = colorRampPalette(c("blue","white","red"))(200),
title = paste("random effect level:",
mod$rLNames[1]), mar=c(0,0,1,0))
```
### Potential dangers
The analysis has worked and we have some really stylish output! But - take screenshots of the output and run it again. Compare your screen shots.
Bayesian solvers don't work the same way as frequentist approaches. With frequentist approaches you get the same result every time, with bayesian approaches a solver explores the parameter space to "find" the right solution. If you do not give time for the solver to coverage on the right solution, you will get a result that is not in the slightest bit reliable!
For a nice overview on assessing Bayesian model convergence see Michael Clark's [bayseian model diagnostics page](https://m-clark.github.io/bayesian-basics/diagnostics.html).
Let's have a look at our traceplots - these are plots which show the Bayesian solvers efforts to converge on the answer for each parameter with each iteration of the model (red and black done the different runs). If they have converged on a solution they should be steady and stable, the coloured lines on the left should overlap and the density plot on the right should be uni-modal.
First for the fixed effects in the model:
```{r ch7_42a, eval=F}
mpost = convertToCodaObject(out)
plot(mpost$Beta)
```
What do you think?
These sampling chains will have to be much longer for these models to converge!
### Further reading
The best place for examples of HMSC analyses right now are package vignettes:
[Getting started with HMSC-R: univariate models](https://cran.r-project.org/web/packages/Hmsc/vignettes/vignette_1_univariate.pdf)
[Getting started with HMSC-R: low-dimensional multivariate models](https://cran.r-project.org/web/packages/Hmsc/vignettes/vignette_2_multivariate_low.pdf)
[Getting started with HMSC-R: high-dimensional multivariate models](https://cran.r-project.org/web/packages/Hmsc/vignettes/vignette_3_multivariate_high.pdf)
[Getting started with HMSC-R: spatial models](https://cran.r-project.org/web/packages/Hmsc/vignettes/vignette_4_spatial.pdf)
### Examples in the literature
[Carvalho Jr, Elildo AR, et al. "Effects of illegal logging on Amazonian medium and large-sized terrestrial vertebrates." Forest Ecology and Management 466 (2020): 118105.](https://www.sciencedirect.com/science/article/pii/S0378112720300803)
[Beirne, Christopher, et al. "Multispecies modelling reveals potential for habitat restoration to re‐establish boreal vertebrate community dynamics." Journal of Applied Ecology 58.12 (2021): 2821-2832.](https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/1365-2664.14020)
```{r ch7_42, echo=F, eval=F}
# Maybe add in the future:
#Bootstrapped confidence intervals from the Bolker link
# Bootstrap some CI's
#set.seed(101)
#m_bb <- bootMer(m1,
# FUN=function(x)
# predict(x,re.form=NA,newdata=newDat,
# type="response"),
# nsim=400)
#
#m_CI<- t(apply(g_bb$t,2,quantile,c(0.025,0.975),na.rm=TRUE))
#newDat <- cbind(newDat, m_CI)
```