AugustT
commented
6 years ago Hi @EichenbergBEF
Some of these are captured in the output, I'll show you where they are:
library(sparta)
set.seed(125)
# Create data
n <- 15000 #size of dataset
nyr <- 20 # number of years in data
nSamples <- 100 # set number of dates
nSites <- 50 # set number of sites
# Create somes dates
first <- as.Date(strptime("2010/01/01", "%Y/%m/%d"))
last <- as.Date(strptime(paste(2010+(nyr-1),"/12/31", sep=''), "%Y/%m/%d"))
dt <- last-first
rDates <- first + (runif(nSamples)*dt)
# taxa are set as random letters
taxa <- sample(letters, size = n, TRUE)
# sites are visited randomly
site <- sample(paste('A', 1:nSites, sep=''), size = n, TRUE)
# the date of visit is selected at random from those created earlier
time_period <- sample(rDates, size = n, TRUE)
results <- occDetModel(taxa = taxa,
site = site,
time_period = time_period,
species_list = 'a',
write_results = FALSE,
n_iterations = 100,
burnin = 10,
thinning = 2,
seed = 125,
modeltype = c("indran", "intercept", "jul_date"))If we want to see how the model was coded we can look at the model in the output. Note that the results object is a list of the same length of the species, in this case only the one (a species called a)
str(results$a, 1)List of 13
$ model :List of 8
..- attr(*, "class")= chr "jags"
$ BUGSoutput :List of 24
..- attr(*, "class")= chr "bugs"
$ parameters.to.save : chr [1:13] "psi.fs" "tau2" "tau.lp" "alpha.p" ...
$ model.file : chr "C:\\Users\\tomaug\\AppData\\Local\\Temp\\RtmpQ3B0NN\\file2a98345279c.txt"
$ n.iter : num 100
$ DIC : logi TRUE
$ SPP_NAME : chr "a"
$ min_year : num 2010
$ max_year : num 2029
$ nsites : int 50
$ nvisits : int 4700
$ species_sites : int 50
$ species_observations: num 511
- attr(*, "class")= chr "occDet"The model element contains model information
results$a$modelJAGS model:
model{
# JAGS code for SPARTA model plus independent random effect prior
# on the state model year effect + intercept
# State model
for (i in 1:nsite){
for (t in 1:nyear){
z[i,t] ~ dbern(muZ[i,t])
logit(muZ[i,t])<- eta.psi0 + a[t] + eta[i]
}}
# State model priors
eta.psi0 <- log(psi0/(1-psi0))
psi0 ~ dbeta(psi0.a, psi0.b)
for(t in 1:nyear){
a[t] ~ dnorm(0, tau.a)
}
tau.a <- 1/(sd.a * sd.a)
sd.a ~ dunif(0, 5)
for (i in 1:nsite) {
eta[i] ~ dnorm(0, tau2)
}
tau2 <- 1/(sigma2 * sigma2)
sigma2 ~ dunif(0, 5)
# Observation model priors
eta.p0 <- log(p0/(1-p0))
p0 ~ dbeta(p0.a, p0.b)
for (t in 1:nyear) {
alpha.p[t] ~ dnorm(0, tau.lp)
}
tau.lp <- 1 / (sd.lp * sd.lp)
sd.lp ~ dunif(0, 5)
# Derived parameters
for (t in 1:nyear) {
psi.fs[t] <- sum(z[1:nsite,t])/nsite
} beta1 ~ dnorm(0, 0.0001)
beta2 ~ dnorm(0, 0.0001)
LL.p ~ dunif(dtype2p_min, dtype2p_max)
### Observation Model
for(j in 1:nvisit) {
y[j] ~ dbern(Py[j])
Py[j]<- z[Site[j],Year[j]]*p[j]
logit(p[j]) <- alpha.p[Year[j]] + eta.p0 + beta1*JulDate[j] + beta2*pow(JulDate[j], 2) + LL.p*logL[j]
} }
Fully observed variables:
JulDate Site Year dtype2p_max dtype2p_min logL nsite nvisit nyear p0.a p0.b psi0.a psi0.b y I'm not sure how familiar with JAGS code you are so this might all mean nothing, but I show it in case it does. The observation model deals with the detection elements.
for (t in 1:nyear) {
alpha.p[t] ~ dnorm(0, tau.lp)
}This deals with the detectability estimates of the species. To get these values from our results we can look them up in the results table:
rownames(results$a$BUGSoutput$summary) [1] "LL.p" "a[1]" "a[2]" "a[3]"
[5] "a[4]" "a[5]" "a[6]" "a[7]"
[9] "a[8]" "a[9]" "a[10]" "a[11]"
[13] "a[12]" "a[13]" "a[14]" "a[15]"
[17] "a[16]" "a[17]" "a[18]" "a[19]"
[21] "a[20]" "alpha.p[1]" "alpha.p[2]" "alpha.p[3]"
[25] "alpha.p[4]" "alpha.p[5]" "alpha.p[6]" "alpha.p[7]"
[29] "alpha.p[8]" "alpha.p[9]" "alpha.p[10]" "alpha.p[11]"
[33] "alpha.p[12]" "alpha.p[13]" "alpha.p[14]" "alpha.p[15]"
[37] "alpha.p[16]" "alpha.p[17]" "alpha.p[18]" "alpha.p[19]"
[41] "alpha.p[20]" "beta1" "beta2" "deviance"
[45] "eta.p0" "eta.psi0" "psi.fs[1]" "psi.fs[2]"
[49] "psi.fs[3]" "psi.fs[4]" "psi.fs[5]" "psi.fs[6]"
[53] "psi.fs[7]" "psi.fs[8]" "psi.fs[9]" "psi.fs[10]"
[57] "psi.fs[11]" "psi.fs[12]" "psi.fs[13]" "psi.fs[14]"
[61] "psi.fs[15]" "psi.fs[16]" "psi.fs[17]" "psi.fs[18]"
[65] "psi.fs[19]" "psi.fs[20]" "tau.a" "tau.lp"
[69] "tau2" You can see here all the alpha.p[X] rows which will give you the estimates for each year. For Julian day the model file shows there is a first order and second order coefficient for Juilan day called beta1 and beta2, you can see these also named in the summary table.
I think this is something we should elaborate on in the vignette. Maybe a section on 'How to extract coefficients'. We cover the occupancy estimates but it is reasonable to expect people to want to investigate the others and at the moment you would need to be able to find and read the JAGS model code to find them. One for us to work on.
@GPowney @CharlieOuthwaite - Please add anything I might have missed.
EichenbergBEF
Hi, in the occDetModel, the detection model is itted according to some options, depending on the model specification (e.g. List Length, Julian Date, Categorical list length...). Is there a way to have a look at the coefficients in the detection model? I have not found this. Thank you very much for your help.