-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathreadme.Rmd
102 lines (90 loc) · 3.42 KB
/
readme.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
---
output:
md_document:
variant: markdown_github
---
###BBGDM
[![Travis-CI Build Status](https://travis-ci.org/skiptoniam/bbgdm.svg?branch=master)](https://travis-ci.org/skiptoniam/bbgdm)
[![codecov.io](https://codecov.io/github/skiptoniam/bbgdm/coverage.svg?branch=master)](https://codecov.io/github/skiptoniam/bbgdm?branch=master)
[![License: GPL v2](https://img.shields.io/badge/License-GPL%20v2-blue.svg)](https://www.gnu.org/licenses/old-licenses/gpl-2.0.en.html)
BBGDM is a R package for running Generalized Dissimilarity Models with Bayesian Bootstrap for parameter estimation. To install package run the following command in your R terminal
```{r,eval=FALSE}
install.packages(c('devtools'))
devtools::install_github('skiptoniam/bbgdm')
```
##### Load the required libraries, we need vegan for the dune dataset.
```{r,error=FALSE,message=FALSE,warning=FALSE}
library(bbgdm)
library(vegan)
```
##### Run bbgdm on the famous dune meadow data
The dune meadow vegetation data, dune, has cover class values of 30 species on 20 sites.
Make the abundance data presence/absence.
```{r}
data(dune)
data(dune.env)
dune_pa <- ifelse(dune>0,1,0)
```
##### Fit a bbgdm
Now we have a species by sites matrix of vegetation data and the associated environmental data for these sites.
```{r,message=FALSE,warning=FALSE,error=FALSE,results='hide'}
form <- ~1+A1
fm1 <- bbgdm(form,dune_pa, dune.env,family="binomial",link='logit',
dism_metric="number_non_shared",spline_type = 'ispline',
nboot=100, geo=FALSE,optim.meth='nlmnib',control = bbgdm.control(cores = 8))
```
##### Print model summary
Here we print out the basic details of the model.
```{r}
print(fm1)
```
##### Plot diagnostics
Using the `diagnostics` function we can extract Random Qunatile Residuals for plotting.
```{r,fig.width = 6, fig.height = 6,fig.align='center'}
resids <- diagnostics(fm1)
par(mfrow=c(2,2))
plot(resids)
```
##### Plot response curves
We can use `as.response` to look at the spline responses in our BBGDM.
The black line represents the median fit, while the grey shaded area is the uncertainty in this fit.
```{r,fig.width = 4, fig.height = 4,fig.align='center'}
response <- as.response(fm1)
par(mfrow=c(1,1))
plot(response)
```
##### Run 'Wald-like' test on parameters
```{r,results='asis',message=FALSE,warning=FALSE,error=FALSE}
library(xtable)
wt <- bbgdm.wald.test(fm1)
tab <- xtable(wt)
print(tab, type = "html")
```
##### Predict BBGDM model
```{r,fig.width = 8, fig.height = 3,fig.align='center'}
#generate some random spatial autocorrelated data.
library(raster)
set.seed(123)
xy <- expand.grid(x=seq(145, 150, 0.1), y=seq(-40, -35, 0.1))
d <- as.matrix(dist(xy))
w <- exp(-1/nrow(xy) * d)
ww <- chol(w)
xy$z <- t(ww) %*% rnorm(nrow(xy), 0, 0.1)
xy$z <- scales::rescale(xy$z,range(dune.env$A1))
coordinates(xy) <- ~x+y
r <- rasterize(xy, raster(points2grid(xy)), 'z')
#give it the same name as variable in bbgdm model.
names(r)<- 'A1'
r2 <- raster(r)
res(r2) <- 0.05
r2 <- resample(r, r2)
#use this layer to predict turnover.
pred.dune.sim.dat <- predict(fm1,r2,uncertainty = TRUE)
#plot the data and the turnover predictions, and error.
colram <- colorRampPalette(c("darkblue","yellow","red"))
colram.se <- colorRampPalette(c('antiquewhite','pink','red'))
par(mfrow=c(1,3),mar=c(3,2,2,6))
plot(r2,main='Simulated A1')
plot(pred.dune.sim.dat[[1]],col=colram(100),main='BBGDM turnover')
plot(pred.dune.sim.dat[[2]],col=colram.se(100),main='BBGDM CV of turnover')
```