Wednesday, July 20, 2016

Running an identity/equivalency test in the ENMTools R package

Okay, let's say we've got two enmtools.species objects: ahli and allogus.  How can we run an identity test?

Here's what we need:


  1. Our two species
  2. Our RasterStack of environmental layers
  3. The type of model we'd like ("glm", "bc", "dm", or "mx", for GLM, Bioclim, Domain, or Maxent)
  4. A formula (GLM only)
  5. The number of reps to perform

So here's how we'd run an identity test using GLM for our two species. 




id.glm = identity.test(species.1 = ahli, species.2 = allogus, env = env, type = "glm", f = presence ~ layer.1 + layer.2 + layer.3 + layer.4, nreps = 99)

Doing 99 reps takes a while, but when you're done, you get an "identity.test" object.  That contains all sorts of useful information.  A quick summary will show you some of it:

id.glm
## 
## 
##  
## 
## Identity test ahli vs. allogus
## 
## Identity test p-values:
##        D        I rank.cor 
##      0.01      0.01      0.01 
## 
## 
## Replicates:
## 
## 
## |          |         D|         I|   rank.cor|
## |:---------|---------:|---------:|----------:|
## |empirical | 0.2221752| 0.4661581| -0.4761597|
## |rep 1     | 0.8883545| 0.9899271|  0.8942366|
## |rep 2     | 0.8486324| 0.9828760|  0.9315827|
## |rep 3     | 0.8227838| 0.9742077|  0.8881490|
## |rep 4     | 0.7255044| 0.9469161|  0.5551645|
plot of chunk unnamed-chunk-14
If you want to access the empirical or replicate models, those are stored in that object as well:

names(id.glm)
[1] "description"               "reps.overlap"              "p.values"                  "empirical.species.1.model" "empirical.species.2.model"
[6] "replicate.models"          "d.plot"                    "i.plot"                    "cor.plot"     

As with building species models, identity.test works pretty much the same for Domain, Bioclim, and Maxent models with the exception that you don't need to supply a formula.

No comments:

Post a Comment