Power analysis with Open-Occupancy Models

[Jesse Whittington]

I am conducting power analysis using data simulation and open-occupancy (colext) models.  In traditional power analyses using glm or glmm's I have used for each iteration a likelihood ratio-test to compare models with and without the explanatory variable "Year".  I have then estimated power as the percentage of iterations that detected a significant decline.  Because colext models automatically include covariates for colonization and extinction, I do not have a null model for comparing models with and without a "Year".  How I should test for significant declines?

Options are:

1.  Use p-values associated with the quasi-extinction covariate.  The problem with this approach is that p-values do not appear to be correlated with the level of simulated decline (see example code below).  

2.  Compare the effect of year in a simple occupancy model.  The problem with this approach is that it doesn't track the repeated measures nature of the data and will underestimate power.

3.  Determine if the 95% confidence intervals for quasi-extinction estimates fall above a threshold of decline (say 0.01, which would reflect greater than 1% decline in the population).

Any help would be greatly appreciated,

Jesse Whittington

Wildlife Biologist

Banff National Park

# Code to simulate data and to estimate extinction rates.  Based on colext vignette.

library(unmarked)

library(plyr)

fnc.simulate.occu <- function(M=50,J=10, T=2, psi=0.7, p.detect=0.9, survival=0.8, colonize=0.0) {

  muZ <- array(dim = c(M, T))

  z <- array(dim = c(M, T)) # Expected and realized occurrence

  y <- array(NA, dim = c(M, J, T)) # Detection histories

  p <- p.detect

  phi <- survival  

  gamma <- colonize 

  # Generate latent states of occurrence

  # First year

  z[,1] <- rbinom(M, 1, psi) # Initial occupancy state

  # Later years

    for(i in 1:M){ # Loop over sites

      for(k in 2:T){ # Loop over years

        muZ[i,k] <- z[i, k-1]*phi + (1-z[i, k-1])*gamma

        z[i,k] <- rbinom(1, 1, muZ[i,k])

      }

    }

  # Generate detection/non-detection data

    for(i in 1:M){

      for(k in 1:T){

        prob <- z[i,k] * p

        for(j in 1:J){

          y[i,j,k] <- rbinom(1, 1, prob)

        }

      }

    }

  yy <- matrix(y, M, J*T)

  year <- matrix(c('01','02'), nrow(yy), T, byrow=TRUE)

  data <- unmarkedMultFrame(y = yy, yearlySiteCovs=list(year=year),  numPrimary=2)

  m <- colext(psiformula= ~1, gammaformula = ~ 1, epsilonformula = ~ 1, pformula = ~ 1, data = data, method="BFGS")

  B <- coef(m)[3]

  SE <- SE(m)[3]  

  p <-   1 - pchisq((B/SE)^2, 1)  # Note: this p-value estimate matches p-values from print(m)

 # back Transform results

  prob.ext <- plogis(B)

  CI <- confint(backTransform(m, type='ext'))

  result <- c(survival, prob.ext, CI,p)

  result <- round(result,3)

  names(result) <- c('survival', 'prob.ext','prob.lower','prob.upper','p.value')

  return(result)

} 

 surv.90 <- ldply(1:100, function(x) fnc.simulate.occu(survival=0.9))            

 surv.50 <- ldply(1:100, function(x) fnc.simulate.occu(survival=0.5))

# Plot extinction probability and p-values

par(mfrow=c(2,1))

  # Extinction Probabilities             

  hist(surv.90$prob.ext, col='gray40', xlim=c(0,1), main='Estimated Extinction Probability')

  hist(surv.50$prob.ext, col='gray80', xlim=c(0,1), add=T)

  legend('topright', c('survival=0.9; mortality=0.1', 'survival=0.5, mortality=0.5'), col=c('gray40','gray80'), bty='n')

  # P-values

  hist(surv.90$p.value, col='gray40', xlim=c(0,1), main='P-Value for Coef. Extinction')

  hist(surv.50$p.value, col='gray80', xlim=c(0,1), add=T)

  legend('topright', c('survival=0.9; mortality=0.1', 'survival=0.5, mortality=0.5'), col=c('gray40','gray80'), bty='n')

[Richard Chandler]

Hi Jesse,

This paper might be useful:

@article {guillera-arroita_lahoz-monfort:2012,
  author = {Guillera-Arroita, Gurutzeta and Lahoz-Monfort, José J.},
  title = {Designing studies to detect differences in species occupancy: power analysis under imperfect detection},
  journal = {Methods in Ecology and Evolution},
  url = {http://dx.doi.org/10.1111/j.2041-210X.2012.00225.x},
  pages = {no--no},
  year = {2012},
}


Richard
_____________________________________
Richard Chandler, post-doc
USGS Patuxent Wildlife Research Center
301-497-5696

From:	Jesse Whittington <jesse.whi...@gmail.com>
To:	unma...@googlegroups.com
Date:	07/16/2012 02:35 PM
Subject:	[unmarked] Power analysis with Open-Occupancy Models
Sent by:	unma...@googlegroups.com

---

[Jesse Whittington]

Thank-you for the reference Richard.  It will be very helpful,

Jesse