Two-species occupancy models

Sometimes we want to include the presence of another species in our occupancy model. For example, leopard occupancy might be lower at sites where tigers are present because of competitive interactions, or it might be lower at sites where a preferred prey species such as muntjac is absent.

These we can model as one-way interactions: leopards avoid tigers but tigers don’t care about leopards; leopards need muntjac but muntjac surely don’t need leopards. In these systems, the leopards are the subordinate species and tigers or muntjac are dominant. (It may seem odd to describe prey as dominant and predators as subordinate, but that’s often the ecological reality.)

If we had reliable information on the presence or absence of the dominant species, we could include that in our leopard model as a normal covariate. But usually we only have detection/nondetection data, and need to allow for sites where the dominant species is present but not detected. We thus need to model occupancy for both species simultaneously. We’ll see how to do it using JAGS on this page.

The example data set

Richmond et al (2010)1Richmond, O.M.W., Hines, J.E., & Beissinger, S.R. (2010) Two-species occupancy models: a new parameterization applied to co-occurrence of secretive rails. Ecological Applications, 20, 2036-2046 describe a study of rails at freshwater marshes in California, USA. They hypothesized that occupancy for the sparrow-sized Black Rail would be lower at sites where the larger Virginia Rail was present.

The data collected by Richmond et al are not available, but the wiqid package has a set of simulated data based on their results. You can load and look at the data with the following code:

library(wiqid)
?railSims
data(railSims)
head(railSims)
#   A1 A2 A3 B1 B2 B3 logArea reeds
# 1  0  0  0  1  1  1  -0.585  TRUE
# 2  1  1  1  0  1  0  -0.231 FALSE
# 3  1  0  0  0  1  1   1.686 FALSE
# 4  0  0  0  1  1  0   0.092  TRUE
# 5  0  0  0  1  1  1   0.155 FALSE
# 6  1  1  1  0  0  0   1.854 FALSE

# Separate the two detection histories
DHA <- as.matrix(railSims[, 1:3])
DHB <- as.matrix(railSims[, 4:6])

The first 3 columns in the data frame are detection/nondetection data for the dominant species (Virginia rail) for 3 surveys at each site, and the next 3 columns are the detection/nondetection data for the Black rail for the same 3 surveys.

logArea and reeds are site covariates which we will discuss later.

A model with no interaction

Without interaction, the model is equivalent to separate models for each species. The JAGS code below is the same as that used for salamanders, except that each line is repeated for each species:

# File name "2sps_noInt.jags"
model{
  # Likelihood
  for(i in 1:nSites) {
    # biological model
    zA[i] ~ dbern(psiA) # 1 if sps A present
    zB[i] ~ dbern(psiB) # 1 if sps B present
    # observation model
    yA[i] ~ dbin(pA * zA[i], n)
    yB[i] ~ dbin(pB * zB[i], n)
  }

  # Priors
  pA ~ dbeta(1, 1) # Uninformative priors
  pB ~ dbeta(1, 1)
  psiA ~ dbeta(1, 1)
  psiB ~ dbeta(1, 1)
}

Running this model is straightforward and produces the same output that we would get if we ran two single-species models. We will go on now to the more interesting model with interaction.

The model with occupancy interaction

Recall that we are interested in one-way interaction: occupancy of species B is affected by the presence of species A, but B has no effect on A. Thus, the code for species A is exactly as before.

For species B, we have 2 values for probability of occupancy, psiB[1] for sites where species A is absent, psiB[2] for sites where it’s present. Which to use at a particular site, i, depends on the value of zA[i], or rather zA[i] + 1.

To get decent estimates for psiB[1] and psiB[2], we need adequate sample sizes for each. That means at least 30 sites where species A is present and 30 sites where it is absent. The railSims data set has 160 sites and roughly half are occupied by the big rail, so we should have no problem.

# File name "2sps_Int.jags"
model{
  # Likelihood
  for(i in 1:nSites) {
    # biological model
    zA[i] ~ dbern(psiA)
    zB[i] ~ dbern(psiB[zA[i] + 1])
    # observation model
    yA[i] ~ dbin(pA * zA[i], n)
    yB[i] ~ dbin(pB * zB[i], n)
  }

  # Priors
  pA ~ dbeta(1, 1) # Uninformative priors
  pB ~ dbeta(1, 1)
  psiA ~ dbeta(1, 1)
  psiB[1] ~ dbeta(1, 1)
  psiB[2] ~ dbeta(1, 1)
}

Now we prepare the data; as before we can aggregate the detections across the 3 survey occasions. Then we can run the model fit:

# Aggregate detection data
yA <- rowSums(DHA)
yB <- rowSums(DHB)

nSites <- nrow(railSims)
zA <- ifelse(yA > 0, 1, NA)
zB <- ifelse(yB > 0, 1, NA)
jagsData <- list(yA = yA, yB = yB, n = 3, nSites = nSites,
    zA = zA, zB = zB)
str(jagsData)

wanted <- c("pA", "pB", "psiA", "psiB")

library(jagsUI)
out <- jags(jagsData, NULL, wanted, model="2sps_Int.jags",
    n.chains=3, n.adapt=1000, n.iter=10000, DIC=FALSE,
    parallel=TRUE)
mc2sps1 <- mcmcOutput(out)
summary(mc2sps1)
#          mean    sd median   l95   u95  Rhat MCEpc
# pA      0.796 0.028  0.797 0.742 0.852 1.000 0.797
# pB      0.717 0.029  0.717 0.657 0.770 1.000 0.846
# psiA    0.467 0.040  0.467 0.389 0.543 1.000 0.564
# psiB[1] 0.818 0.045  0.821 0.729 0.902 0.999 0.690
# psiB[2] 0.391 0.057  0.390 0.279 0.502 1.000 0.605

diagPlot(mc2sps1, max=5)
postPlot(mc2sps1)
postPlot(mc2sps1, "psiB", xlim=c(0.25, 0.95))

The diagnostic plots and posterior plots (not shown) look fine. Plotting the posterior distributions of the two psiB parameters with the same x axis clearly shows the difference in occupancy, with black rail occupancy higher where the Virginia rail is absent (psiB[1]) than where it is present (psiB[2]).

Interaction and detection

One species may behave differently when the other is present, and this may affect the probability of detection. Rails are secretive but have loud calls, so are mostly detected from calls. A change in the frequency of calling will affect detection probability.

We can model this by having 2 values for pA and pB, the first used for sites where the other species is absent and the second when it’s present. The observation model then becomes:

yA[i] ~ dbin(pA[zB[i] + 1] * zA[i], n)
yB[i] ~ dbin(pB[zA[i] + 1] * zB[i], n)

A further possibility is that such effects may be occasion-specific: in this case, that the small rails keep quiet (and are not detected) when the big rails are making a noise (and are detected). We then have 3 values for pB depending on whether:

  • the big rails are absent, zA = 0
  • the big rails are present but not calling (= not detected), zA = 1, DHA = 0
  • the big rails are present and calling (= detected), zA = 1, DHA = 1

We can’t use the aggregated data for this, as we have to deal with the survey occasions one by one. Here is the JAGS code:

# File name "2sps_detect.jags"
model{
  # Likelihood
  for(i in 1:nSites) {
    # biological model
    zA[i] ~ dbern(psiA)
    zB[i] ~ dbern(psiB[zA[i] + 1])
    # observation model
    for(j in 1:n) {
      DHA[i, j] ~ dbern(pA[zB[i] + 1] * zA[i])
      DHB[i, j] ~ dbern(pB[zA[i] + DHA[i, j] + 1] * zB[i])
    }
  }

  # Priors
  pA[1] ~ dbeta(1, 1) # spsB absent
  pA[2] ~ dbeta(1, 1) # spsB present
  pB[1] ~ dbeta(1, 1) # spsA absent
  pB[2] ~ dbeta(1, 1) # spsA present, not detected
  pB[3] ~ dbeta(1, 1) # spsA present, detected
  psiA ~ dbeta(1, 1)
  psiB[1] ~ dbeta(1, 1)
  psiB[2] ~ dbeta(1, 1)
}

And here’s the code to prepare everything and run the model:

jagsData <- list(DHA = DHA, DHB = DHB, n = 3, nSites = nSites,
    zA = zA, zB = zB)
str(jagsData)

wanted <- c("pA", "pB", "psiA", "psiB")

out <- jags(jagsData, NULL, wanted, model="2sps_detect.jags",
  n.chains=3, n.adapt=1000, n.iter=10000, DIC=FALSE,
  parallel=TRUE)
mc2sps2 <- mcmcOutput(out)
summary(mc2sps2)
#          mean    sd median   l95   u95  Rhat MCEpc
# pA[1]   0.846 0.039  0.848 0.766 0.919 1.000 1.015
# pA[2]   0.729 0.053  0.732 0.623 0.829 1.000 1.129
# pB[1]   0.809 0.029  0.810 0.750 0.864 1.001 0.834
# pB[2]   0.679 0.104  0.685 0.482 0.881 1.000 1.156
# pB[3]   0.350 0.069  0.349 0.218 0.486 0.999 1.271
# psiA    0.469 0.040  0.469 0.392 0.547 1.001 0.579
# psiB[1] 0.802 0.044  0.804 0.716 0.885 1.000 0.531
# psiB[2] 0.481 0.081  0.476 0.331 0.646 1.000 1.403

diagPlot(mc2sps2)
plot(mc2sps2)
plot(mc2sps2, "pA", xlim=c(0.6, 0.95))
plot(mc2sps2, "pB", xlim=c(0.2, 0.95))

The model runs quickly and the diagnostic plots and posterior plots (not shown) look fine.

As we suspected, the little rails are more difficult to detect when the big rails are around, and especially when they are actively calling.

More surprising is that the big rails behaviour changes when the small ones are present, making them less detectable. But the difference is small and the posterior distributions overlap, so let’s check on the difference:

diff <- apply(mc2sps2$pA, 1, diff)
postPlot(diff, compVal=0)

So indeed, we can be 95% certain that detection does decline, based on our data and uniform priors.

We’ll include these observation process interactions in the models that follow.

Model with covariates

The area of the habitat patch may affect occupancy (especially for the big rail which has a larger home range), but the difference in occupancy depends on the ratio of the areas of the sites, rather than the difference in area. Hence we work with the logarithm of the area; logArea is the covariate standardized to mean 0 and SD 1.

We incorporate area into the model with a logistic link:
logit(psiA[i]) <- a0 + aArea * area[i]

Richmond et al (2010) also hypothesized that the presence of reed-beds in the wetland would give the little rails an advantage and allow them to coexist with the big rails.

Since the reeds covariate has only 2 values, TRUE/FALSE, we could implement this by having 3 values for psiB: species A absent; A present with without reeds; A present with reeds. Picking out the right value to use for each is tricky, however.

reeds is a 0/1 vector, so reeds+1 is a 1/2 vector. Multiply by zA and we have a 0/1/2 vector, with 0 when species A is absent. Add 1 to get a 1/2/3 vector we can use as the index into a 3-valued psiB vector.

Here’s the complete JAGS model:

# File name "2sps_covars.jags"
model{
  # Likelihood
  for(i in 1:nSites) {
     # Ecological model
     logit(psiA[i]) <- a0 + aArea * area[i]
     zA[i] ~ dbern(psiA[i])
     zB[i] ~ dbern(psiB[(reeds[i]+1)*zA[i] + 1])
     # Observation model
     for(j in 1:n) {
       DHA[i, j] ~ dbern(pA[zB[i] + 1] * zA[i])
       DHB[i, j] ~ dbern(pB[zA[i] + DHA[i, j] + 1] * zB[i])
     }
  }

  # Priors
  pA[1] ~ dbeta(1, 1) # as in previous model
  pA[2] ~ dbeta(1, 1)
  pB[1] ~ dbeta(1, 1)
  pB[2] ~ dbeta(1, 1)
  pB[3] ~ dbeta(1, 1)

  a0 ~ dunif(-5, 5)
  aArea ~ dunif(-5, 5)
  psiB[1] ~ dbeta(1, 1)  # when sps A absent
  psiB[2] ~ dbeta(1, 1)  # when sps A present, no reeds
  psiB[3] ~ dbeta(1, 1)  # when sps A present, reeds
}

Let’s check that we have enough data to estimate the reeds effects: only sites with species A present count, and we want at least 30 sites with reeds and 30 without. Big rails were detected at 74 sites, 32 with reeds and 42 without, so we have no worries there.

Now we can prepare the data and run the model. In the version of wiqid on CRAN, railSims$logArea is a 1-column matrix instead of a vector; we can fix that with c(). We also change railSims$reeds from TRUE/FALSE to 1/0.

jagsData <- list(DHA = DHA, DHB = DHB, n = 3, nSites = nSites,
    zA = zA, zB = zB, area=c(railSims$logArea),
    reeds=as.numeric(railSims$reeds))
str(jagsData)

wanted <- c("pA", "pB", "a0", "aArea", "psiB")

out <- jags(jagsData, NULL, wanted, model="2sps_covars.jags",
  n.chains=3, n.adapt=1000, n.iter=10000, DIC=FALSE,
  parallel=TRUE)
mc2sps3 <- mcmcOutput(out)
summary(mc2sps3)
#           mean    sd median    l95   u95  Rhat MCEpc
# pA[1]    0.834 0.041  0.837  0.752 0.912 0.998 1.064
# pA[2]    0.738 0.053  0.741  0.630 0.835 0.998 1.184
# pB[1]    0.809 0.029  0.810  0.751 0.862 0.999 0.813
# pB[2]    0.691 0.098  0.696  0.496 0.873 1.001 1.015
# pB[3]    0.345 0.070  0.343  0.211 0.484 0.999 1.409
# a0      -0.150 0.195 -0.150 -0.522 0.243 1.000 0.708
# aArea    1.477 0.266  1.467  0.954 1.998 1.000 0.708
# psiB[1]  0.803 0.044  0.806  0.718 0.888 1.001 0.574
# psiB[2]  0.313 0.086  0.306  0.155 0.487 1.000 1.068
# psiB[3]  0.701 0.116  0.698  0.492 0.950 1.001 1.442

diagPlot(mc2sps3)
plot(mc2sps3)
plot(mc2sps3, "psiB", xlim=c(0.1, 1))
plot(mc2sps3, "aArea")

As before, the runs quickly and converges nicely; the diagnostic and posterior plots (not shown) reveal no problems.

Reed beds clearly enable the smaller Black rail to survive in wetlands where the Virginia rail occurs, giving occupancy estimates almost as high as sites without the big rail.

Wetland area is an important determinant of occupancy for the Virginia rail: a coefficient value of 1.5 is high for a standardised covariate in a logistic regression. We’ll do a plot of probability of occupancy against log(area):

range(railSims$logArea) # +/- 2.5
xx <- seq(-2.5, 2.5, length=101)
toPlot <- matrix(NA, 101, 3)
colnames(toPlot) <- c("est", "lo", "up")
for(i in 1:101){
  psi <- plogis(mc2sps3$a0 + mc2sps3$aArea * xx[i])
  toPlot[i, ] <- c(mean(psi), hdi(psi))
}

plot(xx, toPlot[, 1], type='n', ylim=range(toPlot),
    xlab="standardised log area", ylab="Occupancy")
polygon(c(xx, rev(xx)), c(toPlot[,2], rev(toPlot[,3])),
    col='skyblue', border=NA)
lines(xx, toPlot[,1], lwd=2)

You can download a ZIP file with all the code here (includes NIMBLE code).

31 thoughts on “Two-species occupancy models”

  1. Febri Widodo emailed me with a question about calculation of the Species Interaction Factor (SIF) (MacKenzie et al, 2004) from the output of wiqid::occ2sps. That’s a mess, but it’s very easy with Bayesian MCMC output. Here’s the calculation based on the mc2sps1 output:
    SIF <- mc2sps1$psiB[,2] / (mc2sps1$psiA*mc2sps1$psiB[,2] +
    (1-mc2sps1$psiA)*mc2sps1$psiB[,1])
    postPlot(SIF)

    1. Hi Mike!

      Thank you for the extremely helpful tutorials. Do you know if there is a way to calculate SIF for each sampling site, or is it only possible to calculate one value across sites?

      1. Zach,
        I’m assuming you have site covariates, so the values of psiA and/or psiB differ across sites. It should be possible to calculate for each site: you’d have to extract the values of psiA for each site, and those for psiB when species A is present and absent. Then plug the relevant values into the same equation for SIF. You’d need to loop over all sites and gather up the SIF values.

        1. Hi Mike,

          Thanks for your helpful response. I was able to add a continuous covariate (Forest Cover) to derive estimates for each site, but my results seem a bit off. For instance, the SIF value I derive from beta estimates is about 1.45, whereas it was just 0.84 for the same species pair in a model with interaction but no covariates. That kind of jump makes me think my model is suspect. Would it be alright to just glance over my JAGS model code? This is the model structure I’m using:

          model{
            # Likelihood
            for(i in 1:nSites) {
               # Ecological model
               logit(psiBA[i]) <- b0 + bFc * fc[i]
               zA[i] ~ dbern(psiA)
               zB[i] ~ dbern(ifelse(zA[i]==1, psiBA[i], psiBa))
               # Observation model
               yA[i] ~ dbin(pA * zA[i], n)
               yB[i] ~ dbin(pB * zB[i], n)
            }
          
            # Priors
            pA ~ dbeta(1, 1)
            pB ~ dbeta(1, 1)
            
            b0 ~ dunif(0, 1)
            bFc ~ dunif(0, 1)
            psiA ~ dbeta(1, 1)
            psiBa ~ dbeta(1, 1)
          }
          1. Zach,
            Re the JAGS code: your priors for b0 and bFc are unusual: they may work ok, you’d have to look at the diagnostic plots, but I suspect they are actually constraining the posterior. dunif(-5,5) is likely better (assuming your fc covariate is standardized). A neat way to specify a prior for the intercept (when all covariates are standardised) is to first specify a prior on the probability scale then convert to logit:
            psiBA0 ~ dbeta(1,1)
            b0 <- logit(psiBA0)

            Re SIF: SIF depends on psiBA, so will have a different value for each site. I don't see where your value of 1.84 comes from. If the model doesn't take too long to run, it might be easiest to include the calculation of site-wise SIF in the JAGS code. Inside the i loop, include
            SIF[i] <- psiBA[i] / (psiA*psiBA[i] + (1-psiA)*psiBa)

  2. Mike,

    Thanks again for the blog posts, they’ve been extremely helpful.
    For the species-interaction model with covariates on psiA (area) and psiB (reed), how might the structure of the model differ if reed was a continuous covariate?

    1. Joshua,
      Good question! Suppose we have a continuous variable reediness which goes from 0 to 100. We simulate some random values and standardize with
      set.seed(27)
      reediness <- runif(nSites, 0, 100)

      reedinessS <- standardize(reediness)

      In the JAGS code, in the "sites" loop, we need
      logit(psiBA[i]) <- b0 + bRd * reediness[i] # when sps A present
      zB[i] ~ dbern(ifelse(zA[i]==1, psiBA[i], psiBa))

      [Edit: the first argument to ifelse should be a condition, so corrected it to zA[i]==1.]
      The first line there calculates psiBA[i] for each site, depending on the value of reediness[i], and the second line uses ifelse to select psiBA[i] or psiBa depending on the value of zA[i]. psiBa is the same for all sites.

      Then we need priors for our new parameters:
      b0 ~ dunif(-5, 5)
      bRd ~ dunif(-5, 5)
      psiBa ~ dbeta(1, 1) # when sps A absent

      In the R code, we need to change jagsData to include reediness=reedinessS and wanted to include b0, bRd and psiBa.

      The rest of the code should then work ok. Of course, reediness is fiction, so we'd expect bRd to be near zero, but with only 160 sites, spurious correlations can easily arise.

  3. Mike,

    One more quick question. I also wondered how the 2-species model structure might change if a categorical factor was included in either the “no-interaction” or “interaction” models. In this case, “Time” is the categorical factor.

    For the no-interaction model: I assumed it would look like the following:

    #Ecological Process
      for(i in 1:nSites) {
        zA[i] ~ dbern(psiA[Time[i]]) 
        zB[i] ~ dbern(psiB[Time[i]]) 
        #Observation Process
        yA[i] ~ dbin(pA * zA[i], n)
        yB[i] ~ dbin(pB * zB[i], n)
      }
      # Priors
      pA ~ dbeta(1, 1) # Uninformative priors
      pB ~ dbeta(1, 1)
      for (i in 1:nCat) {
        psiA[i] ~ dbeta(1, 1)
        psiB[i] ~ dbeta(1, 1)
      }
    }

    And for the Interaction Model, I assumed it would look something like the following (I’m much less confident in this
    structure though!):

    for(i in 1:nSites) {
        # Ecological Model
        zA[i] ~ dbern(psiA[Time[i]])
        zB[i] ~ dbern(psiB[zA[i] + 1])
        
        # Observation model
        yA[i] ~ dbin(pA * zA[i], n)
        yB[i] ~ dbin(pB * zB[i], n)
      }
      
      # Priors
      pA ~ dbeta(1, 1) 
      pB ~ dbeta(1, 1)
      for(i in 1:nCat){
        psiA[i] ~ dbeta(1, 1)
        psiB[1] ~ dbeta(1, 1)
        psiB[2] ~ dbeta(1, 1)
      }
    }

    Any assistance on this would be much appreciated, Thanks!

    1. I edited your post to include some formatting which you can’t include and corrected typos (italicized).

      I need to run things to see if they work! Let’s assume the first 40 sites for the Rails data were surveyed in (say) January, the next 40 in March, then June and September, and we think occupancy changes through the year. So our Time covariate will be:
      Time <- rep(1:4, each=40)

      Your first model works fine (after corrections) and gives 4 values for each of psiA and psiB.

      The second model won’t work: you have psiB[1] and psiB[2] inside the nCat loop, so will get an error saying you are redefining a node. Take those outside the nCat loop will prevent the error, but psiB will not change with Time.

      psiB needs to be a matrix with different values for each Time and for each zA. Change the ecological model to
      zB[i] ~ dbern(psiB[zA[i] + 1, Time[i]])
      and in the nCat loop have
      psiB[1, i] ~ dbeta(1, 1)
      psiB[2, i] ~ dbeta(1, 1)

      Do remember in all this that our sample sizes are getting pretty small; ok if your original sample is huge.

  4. A colleague emailed me asking “how [can we] tweak the model to deal with NA data?”

    NAs in the detection data are not a problem, provided the two species have the same pattern of NAs. If the analysis used the detection matrices, DHA and DHB above, everything should just work.

    If you are using the aggregated data, yA and yB, it’s a tad more complicated. You need to include na.rm=TRUE in the call to rowSums. Also the number of occasions is no longer the same for all sites, so you need to do:
    n <- rowSums(!is.na(DHA))
    and in the JAGS code use n[i] instead of n.

    JAGS doesn't allow NAs in covariates. If these are survey-level covariates (eg, survey date) and the NAs correspond to NAs in the detection matrices, the values don't matter, so just plug in 0 or the mean of the other values or whatever.

  5. Hi Mike,

    I have a problem with the script used to plot the relation between occupancy and standardized logArea. When I run the plot it returs the following warning:

    Error in plot.window(…) : need finite ‘ylim’ values

    I tryed to solved it in several ways but I didn’t find a solution.
    Moreover, just out of curiosity, in the “model with covariates” scenario I would like to consider four cases rather than three:

    1. absence of species A, no reeds;
    2. absence of species A, reeds;
    3. presence of species A, no reeds;
    4. presence of species A, reeds.

    Would I be correct in calculating the zB[i] variable as follows?

    zB[i] ~ dbern(psiB[(reeds[i]+2)*zA[i] + 1])

    Many thanks,
    Marcello

    1. Re the error: The code has ylim=range(toPlot), so do range(toPlot) and see if there are NAs in there. If so, something went wrong with the code for the values in toPlot; do head(toPlot) to see if specific columns have issues.

      Re 4 cases: Your code won’t work, [(reeds[i]+2)*zA[i] + 1] will return 1 if zA == 0, otherwise 3 or 4.

      Better to use a 2×2 matrix for psiB, with columns for Sps A present/absent and rows for reeds present/absent (or the other way around). The priors then look like this:
      psiB[1,1] ~ dbeta(1, 1) # when sps A absent, no reeds
      psiB[2,1] ~ dbeta(1, 1) # when sps A absent, reeds
      psiB[1,2] ~ dbeta(1, 1) # when sps A present, no reeds
      psiB[2,2] ~ dbeta(1, 1) # when sps A present, reeds

      Then in the Ecological model, use
      zB[i] ~ dbern(psiB[reeds[i]+1, zA[i]+1])
      The rest of the code is the same, but you will get 4 values for psiB in the output.

  6. Hi Mike,

    Re: Re the error: I don’t know if I correctly interpreted what you meant but the toplot vector has only NA values as was created as matrix of NA values. So, I don’t know how to solve the problem.

    Re: Re for cases: many thanks!!

    Best,
    Marcello

  7. Hi Mike,

    Do you know how to create a plot akin to the standardized log area-occupancy plot, except with abundance (count data) instead of presence/absence?

    I am using a static binomial N-mixture model with two species and directional interactions with PLN distribution (expansion of Waddle et al. 2010 in AHM2 chapter 8.5.2). Script for the final model is attached.

    Thanks!

    PLN N-mixture two species

  8. Hi Mike

    I have a problem with the code you provided to count the occupancy of the two species, psiA, psiBA(below code). I want to know how to write the code when I want to know the covariates’ effect on the detection, like pA, pB, rB, rA,rBA. can you give me some sample codes?

    logit(psiBA[i]) <- b0 + bRd * reediness[i] # when sps A present
    zB[i] ~ dbern(ifelse(zA[i]==1, psiBA[i], psiBa))

    thanks

    1. Well yes, you’ve answered the question! If it’s a straightforward continuous covariate just use a logistic expression. But you would have to turn pA or pB into a matrix with a row for each site. Eg: this inside the i loop:
      logit(pA[i,1]) <- c0 + cWh * whatever[i]
      logit(pA[i,2]) <- d0 + dWh * whatever[i]

      and provide priors for the new coefficients.

      Where pA[zB[i] + 1] appears in the code, you’d need to use pA[i, zB[i] + 1] to index into the matrix.

      1. Hi Mike

        Thanks for your reply, but there still has some questions confused me. I write the code for the coefficient of detection, the model doesn’t work. The code I have written is below, could you tell me what is wrong and how to correct it?

        mod_string = ” model{
        # Likelihood
        for(i in 1:nSites) {
        # Ecological model
        logit(psiA[i]) <- a0 + aHM * HM[i] + aPeopleEr * peopleEr[i]+aGrazeEr * grazeEr[i]+aElevation * elevation[i]+aDistSettle *distSettle[i]

        logit(psiBA[i]) <- b0 + bHM * HM[i] + bPeopleEr * peopleEr[i]+bGrazeEr * grazeEr[i]+bElevation * elevation[i]+bDistSettle *distSettle[i] # when sps A present

        logit(psiBa[i]) <- c0 + cHM * HM[i] + cPeopleEr * peopleEr[i]+cGrazeEr * grazeEr[i]+cElevation * elevation[i]+cDistSettle *distSettle[i]

        logit(pA[i,1]) <- d0 + dHM * HM[i]
        logit(pA[i,2]) <- e0 + eHM * HM[i]

        zA[i] ~ dbern(psiA[i])
        zB[i] ~ dbern(ifelse(zA[i]==1, psiBA[i], psiBa[i]))

        pA[i,1] ~ dbeta(1, 1)

        # Observation model
        for(j in 1:n) {

        DHA[i, j] ~ dbern(pA[zB[i] + 1] * zA[i])
        DHB[i, j] ~ dbern(pB[zA[i] + DHA[i, j] + 1] * zB[i])
        }
        }

        # Priors
        pA[1] ~ dbeta(1, 1) # spsB absent
        pA[2] ~ dbeta(1, 1) # spsB present
        pB[1] ~ dbeta(1, 1) # spsA absent
        pB[2] ~ dbeta(1, 1) # spsA present, not detected
        pB[3] ~ dbeta(1, 1) # spsA present, detected

        a0 ~ dunif(-5, 5)
        aHM ~ dunif(-5, 5)
        aPeopleEr ~ dunif(-5, 5)
        aGrazeEr ~ dunif(-5, 5)
        aElevation ~ dunif(-5, 5)
        aDistSettle ~ dunif(-5, 5)
        b0 ~ dunif(-5, 5)
        bHM ~ dunif(-5, 5)
        bPeopleEr ~ dunif(-5, 5)
        bGrazeEr ~ dunif(-5, 5)
        bElevation ~ dunif(-5, 5)
        bDistSettle ~ dunif(-5, 5)
        c0 ~ dunif(-5, 5)
        cHM ~ dunif(-5, 5)
        cPeopleEr ~ dunif(-5, 5)
        cGrazeEr ~ dunif(-5, 5)
        cElevation ~ dunif(-5, 5)
        cDistSettle ~ dunif(-5, 5)
        d0 ~ dunif(-5, 5)
        dHM ~ dunif(-5, 5)
        e0 ~ dunif(-5, 5)
        eHM ~ dunif(-5, 5)

        }"

        Truly I am not know very well for the Bayesian models, and also not skill at writing code in R. If you don't mind, you can also send me the email, my email is: jichengp321@163.com

        So thanks again for your help.

        Best wishes

        Ji

        1. I replied by email. “My code doesn’t work” is not of much interest to anyone else!

  9. Hi Mike,

    When it comes to the co-occurrence models with several covariates, I would like to ask if you, or anybody, know(s) of some function to draw every possible variable combination (a bit like the dredge function from the MuMIn package) that one can use with co-occurrence jags models.

    Thanks,

    Luis

    1. Luis,
      A couple of general points on dredging: (1) Throwing everything in and hoping AIC or similar will sort them out is a terrible idea. You are wide open to spurious correlations in your data, and information criteria or LOO cannot be relied upon to pick those out, see here. You are far better off considering what models would be biologically reasonable, not forgetting models with quadratic terms or interactions, which are ignored when dredging. (2) Once you’ve fitted all the models in JAGS, what do you do next? They are hierarchical models so you can’t use DIC; you’d have to use WAIC or LOOAIC, and those are not easy to code. You are better off using a Bayesian model averaging approach, see Hooten & Hobbs (2015).

      Specifically with 2-species models, the basic model with occupancy interaction but no covariates or detection interactions has 5 parameters, and they are binary (occupied/not occupied). You need lots of sites to estimate those properly. If you add in “several” covariates, you’ll need data for many hundreds of sites.

      If you really want to go ahead with this, the wiqid package has a function allCombinations which will do it for you. Plug the resulting formulae into model.matrix to get model matrices. Write your JAGS code using vector multiplication so that the same code can be used no matter how many columns in the model matrix, see the JAGS code here.

  10. Hi all.

    Take the modelling steps with covariates, for instance. Is there a way to reach the individual occupancy estimates of each trapping site?

    Thanks in advance.

      1. The JAGS output gives you psiBA[i] (applicable when species A is present at site i) and psiBa[i] (when species A is absent). You may want an “overall” probability that species B is present. You can get that as a weighted mean of psiBA and psiBa, where the weights are just the probability that species A is present or absent. So
        psiB[i] <- psiBA[i] * psiA[i] + psiBa[i] * (1 - psiA[i])
        You can do that inside the JAGS code, or pull out the MCMC chains for the relevant nodes and do the arithmetic after the JAGS run.

        1. Hi Mike,
          I think perhaps I was not very clear, but no problem. I found the answer to my question and you ended up answering to another question I had in my regarding the estimate of psiB[i].
          Thank you very much

    1. Luis, Not quite sure what you want here. The psi[i] parameters give the probability of occupancy for a site with the same covariate values as site i. I suspect you want the conditional occupancy probability, where conditional = conditional on the data. So if you detected the species at site i, you know it is occupied, so psiCond = 1. If you surveyed the site but did not detect the species, it may still be occupied: if you made many visits with high detection probability, intuitively it probably isn’t occupied; if you made only 1 or 2 visits with low probability of detection, you don’t have much information, so psiCond[i] will be only a little lower than psi[i]. If you didn’t survey the site at all, then psiCond[i] = psi[i]. In a JAGS analysis, psiCond[i] is given by the mean value of z[i].

      For a two-species model this would apply to psiCondA[i] = zA[i] and psiCondB[i] = zB[i].

  11. Hi Mike,
    What is the best way to assess model fit for the two-species model? Can we use GOF test for the two-species occupancy model?
    I saw cross-validation posted for multi-species, but with a small number of sites ~60 sites (with a low number of detections), it won’t work, right?
    Regards,
    Soy

    1. Soy, Should be no problem using the usual GOF tests. You would calculate separate discrepancies for species A and B then add them up to get an overall discrepancy figure. Cross-validation is good for comparing models, but not always clear what value of, say, AUC corresponds to a “good” fit.

Leave a Reply to admin Cancel reply

Your email address will not be published.

The maximum upload file size: 1 MB. You can upload: image, document, text, archive. Drop file here