Here's what we need:
- Our two species
- Our RasterStack of environmental layers
- The type of model we'd like ("glm", "bc", "dm", or "mx", for GLM, Bioclim, Domain, or Maxent)
- A formula (GLM only)
- 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|
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.
Hi,
ReplyDeleteHow to interpret the result obtained from this?
Thank you
Hi,
ReplyDeleteHow to interpret the results obtained from this?
Thank you