Monday, July 18, 2016

Using an enmtools.species object to build an ENM

Now that we've got our enmtools.species object and have assigned it presence and background data and a species name, we can use it to build models very simply!  For this functionality, ENMTools is basically just acting as a wrapper for dismo, using those functions to actually build models.  At present ENMTools only has interfaces for GLM, Maxent, Bioclim, and Domain, but that will change with time.

So let's use our enmtools.species object "ahli" to build a quick Bioclim model.  I've got a RasterStack object made up of four environmental layers.  It's named "env", and the layers are just "layer.1", etc.

Let's build a model!

ahli.bc = enmtools.bc(species = ahli, env = env)

For Bioclim, Domain, and Maxent, it's that easy!  ENMTools extracts the presence and background data from env using the data stored in the species object, builds a model, and returns some lovely formatted output.

For GLM we need to supply a formula as well, but other than that it's idential.


ahli.glm = enmtools.glm(f = pres ~ layer.1 + layer.2 + layer.3 + layer.4, species = ahli, env = env)

Let's look at the output:

ahli.glm
## 
## 
## Formula:  presence ~ layer.1 + layer.2 + layer.3 + layer.4
## 
## 
## 
## Data table (top ten lines): 
## 
## | layer.1| layer.2| layer.3| layer.4| presence|
## |-------:|-------:|-------:|-------:|--------:|
## |    2765|    1235|    1174|     252|        1|
## |    2289|    1732|     957|     231|        1|
## |    2158|    1870|     983|     253|        1|
## |    2207|    1877|     967|     259|        1|
## |    2244|    1828|     945|     249|        1|
## |    2250|    1766|     919|     235|        1|
## |    2201|    1822|     978|     277|        1|
## |    2214|    1786|     986|     284|        1|
## |    2287|    1722|     992|     266|        1|
## |    2984|     965|    1311|     237|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.67171  -0.20485  -0.14150  -0.09528   3.08762  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept) 46.178922  24.777923   1.864   0.0624 .
## layer.1     -0.013347   0.006276  -2.127   0.0334 *
## layer.2     -0.011985   0.006612  -1.813   0.0699 .
## layer.3      0.003485   0.006586   0.529   0.5967  
## layer.4     -0.009092   0.021248  -0.428   0.6687  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 164.58  on 1015  degrees of freedom
## Residual deviance: 150.15  on 1011  degrees of freedom
## AIC: 160.15
## 
## Number of Fisher Scoring iterations: 8
## 
## 
## 
## Suitability:  
## class       : RasterLayer 
## dimensions  : 418, 1535, 641630  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -86.90809, -74.11642, 19.80837, 23.2917  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 6.419793e-08, 0.999983  (min, max)
plot of chunk unnamed-chunk-8
You get pretty, formatted output and a nice plot as well.  Next up: identity/equivalency tests!



2 comments:

  1. It's really demanded to create this model? Because I already have my models from maxent and I thought I could use them to run my background.

    ReplyDelete
  2. Yes, using the R version you need to create enmtools.species objects. That is not the case with the Perl version, but it doesn't do many of the things that the R version does.

    ReplyDelete