occupancy estimates for single species multiseason modelling (colext)

[Vishnu r menon]

Hi all,

I'm a running a multiseason (two seasons- prefire and postfire) model for foxes. I have a covariate called treatment which are sites that are either fox 'baited' or 'unbaited'. Now I am running a simple model with just treatment affecting initial occupancy:

sitecovariate <-data.frame(treatment,region)

umf <- unmarkedMultFrame(y= fox_det, yearlySiteCovs = list( baiting=baiting),

                         siteCovs = sitecovariate,

                          obsCovs = list(time.det= time, site = site), numPrimary = 2)

mod <- colext(~treatment, ~1,~1,~1, umf) 

The model works fine (summary below), but I'm unsure as to how to get occupancy estimates across baited/unbaited for both seasons. So my output should look like:

Summary of the model:

It says that I can access occupancy probability using @projected.mean or @smoothed.mean, but I do not know which value refers to what in the output. m...@projected.mean gives me the output:

Does anyone know how to interpret this output @projected output and to figure what the respective values refer to? Apologies if the question is a bit too silly/confusing, I'm new to unmarked.

Thanks,

Vishnu

[Richard Chandler]

Hi Vishnu,

Use the 'projected' values if you want predictions of the *expected* value of occupancy. These are based on the maximum likelihood estimates and the formula:

PrOccupied(i,t) = PrOccupied(i,t-1)*[1-PrExtinction(i,t-1)] + [1-PrOccupied(i,t-1)]*PrColonized(i,t-1)

Use the 'smoothed' values if you want predictions of the *realized* value of the latent occupancy state, conditional on the entire sequence of observations. These smoothed values are computed using the forward-backward algorithm:

https://en.wikipedia.org/wiki/Forward%E2%80%93backward_algorithm

In other words, use the smoothed values if you want estimates of occupancy at the specific sites you sampled during the specific time period when you sampled. Use the projected values if you want standard predictions that might be useful for assessing what would happen if you repeated the study under similar conditions.

A great early reference on hidden Markov models and smoothing, in the context of occupancy models, is:

Fiske, I.J., Royle, J.A. and Gross, K., 2014. Inference for finite-sample trajectories in dynamic multi-state site-occupancy models using hidden Markov model smoothing. Environmental and ecological statistics, 21(2), pp.313-328.

Richard

---

From: unma...@googlegroups.com <unma...@googlegroups.com> on behalf of Vishnu r menon <vishnur...@gmail.com>
Sent: Wednesday, February 24, 2021 12:00 AM
To: unmarked <unma...@googlegroups.com>
Subject: [unmarked] occupancy estimates for single species multiseason modelling (colext)

 

[EXTERNAL SENDER - PROCEED CAUTIOUSLY]

--
You received this message because you are subscribed to the Google Groups "unmarked" group.
To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/851d37b1-11f3-4fa6-8e5c-492e608a78d7n%40googlegroups.com.

[Vishnu r menon]

Hi Richard,

Thanks for your reply and explanation.

I know that I want to use projected, but my issue was that I couldn't clearly interpret the output of projected values from my model (named as mod here):

mod <-  colext(~treatment, ~1,~1,~1, umf) 

projected(mod) gives the output

I understand that 1 and 2 refers to season one and season two, but does 'unoccupied' and 'occupied' refer to the occupancies of sites that are occupied/not occupied by that species? If so, how do I get occupancy estimates relating to the variable affecting initial occupancy, that being 'treatment' in this case. treatment is a dataframe with whether a site is 'baited' or 'unbaited'. So I would like to get the occupancy estimates in baited and unbaited sites in season 1 and then in season 2. I hope my question makes sense, thanks a lot!

Regards,

Vishnu

[Richard Chandler]

Hi Vishnu,

In that case, you can use `predict` to give you estimates, SEs, and CIs for initial occupancy in the baited and unbaited sites. See the examples here:

help(unmarked)

The predictions will just be:

Pr(occupied at t=1 at baited site) = plogis(-1.46)

Pr(occupied at t=1 at unbaited site) = plogis(-1.46+0.65)

Richard

---

From: unma...@googlegroups.com <unma...@googlegroups.com> on behalf of Vishnu r menon <vishnur...@gmail.com>

Sent: Wednesday, February 24, 2021 7:46 PM
To: unmarked <unma...@googlegroups.com>
Subject: Re: [unmarked] occupancy estimates for single species multiseason modelling (colext)

 

To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/7b31a195-8a4d-4726-84d0-61e50dd4ef45n%40googlegroups.com.

[Vishnu r menon]

Hi Richard,

Thanks for your reply. Predict gives me initial occupancies (pre-fire) for baited and unbaited sites which I understand. But if I want to get occupancy estimates for baited and unbaited sites postfire (season 2), I am not sure how to do it because @projected gives me occupancy probabilities for sites that are occupied and unoccupied instead of baited and unbaited. I couldn't find a clear explanation for this in unmarked unfortunately. Or does the 'occupied' and 'unoccupied' refer to something else?

@projected output is :

Thanks again for your help.

Regards,

Vishnu

[Ken Kellner]

Vishnu,

If I'm understanding you correctly what you ultimately want is a matrix of occupancy with two rows (one for baited and one for unbaited) and two columns (one for season 1 and one for season 2). If that's correct I think posterior prediction with Empirical Bayes is a solution. Rather than explain here see this previous thread:

https://groups.google.com/g/unmarked/c/l9CCdAplmAo/m/H9r42ErIBAAJ

The code in that thread is outdated since this feature is now officially part of unmarked. I wrote up a quick example below using the built-in frogs dataset which should get you started. See also ?predict, specifically the final section in the help file related to unmarkedRanef objects.

Ken

-------------------------------

library(unmarked)

data(frogs)

umf <- formatMult(masspcru)[1:200,]

# Make fake categorical covariate "x" with two levels a and b

set.seed(123)

siteCovs(umf) <- data.frame(x = factor(sample(letters[1:2],numSites(umf), replace=T)))

head(siteCovs(umf))

# Fit model with effect of x on initial occupancy

(fit <- colext(~x, ~1, ~1, ~1, umf))

# Get estimates of latent occupancy state (z) at each site using empirical bayes

r <- ranef(fit)

# Now we can generate an empirical bayes posterior predictive distribution

# for mean occupancy in each group a and b

# Function that takes argument zmat, which will be a matrix M x T of latent occupancy

# zmat will be 200 x 7 in this case since there are 200 sites and 7 primary periods

split_z_by_x <- function(zmat){

  T <- ncol(zmat) # primary periods

  x <- siteCovs(umf)$x

  z_a <- zmat[x=="a",] # sites where x == "a"

  z_b <- zmat[x=="b",] # sites where x == "b"

  #output object: 2 rows (for categories a and b) and 7 columns (for 7 primary periods)

  out <- matrix(NA, nrow=2, ncol=T)

  out[1,] <- colMeans(z_a) # mean of posterior of z for sites with x == "a"

  out[2,] <- colMeans(z_b) # same for b

  rownames(out) <- c("a","b")

  colnames(out) <- 1:T

  out

}

# Use predict method for unmarkedRanef, supplying the function, and get 100 samples

pr <- predict(r, func=split_z_by_x, nsim=100)

# Output will be 2 groups x 7 periods x 100 simulations

dim(pr)

#Get mean of posterior for each group and primary period

(post_mean <- apply(pr, c(1,2), mean))

# Get 95% CI for each group and primary period

post_ci <- apply(pr, c(1,2), quantile, c(0.025,0.975))

# Plot means and 95% CIs

plot(x=1:7-0.1, y=post_mean[1,], col='blue', type='b', ylim=c(0.82,1),

    ylab="Proportion occupied", xlab="Primary period")

for (i in 1:7){

  segments(i-0.1, post_ci[1,1,i], i-0.1, post_ci[2,1,i], col='blue')

}

lines(1:7+0.1, post_mean[2,], type='b', col='red')

for (i in 1:7){

  segments(i+0.1, post_ci[1,2,i], i+0.1, post_ci[2,2,i], col='red')

}

legend("bottomright", col=c("blue","red"), legend=c("a","b"),pch=21,

       title="Value of x")

[Richard Chandler]

If you're after the realized values of occupancy, then I agree with Ken.

If you want the expected values, I would use the equations I provided earlier. The expected values for each site and each season are computed and stored in the fitted model object. You can access these like so:

```

library(unmarked)

example(colext)

projected(fm, mean=FALSE) # The second row contains the expected occupancy probs

```

To get SEs and CIs in unmarked, you would need to bootstrap. 

Richard

---

From: unma...@googlegroups.com <unma...@googlegroups.com> on behalf of Ken Kellner <ken.k...@gmail.com>
Sent: Sunday, February 28, 2021 11:55 AM
To: unmarked <unma...@googlegroups.com>

To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/fe03c9fd-35c0-425a-889b-509f31bbecf7n%40googlegroups.com.