Standard errors from recursive equation

[mFidino]

Hi all,

I am using the recursive equation to calculate future occupancy rates as a function of colonization and extinction probabilities (in the same way as done with the House finch expansion example from the single-species, multiple season chapter of MacKenzie et al. 2006). However, I am uncertain how to calculate the standard errors associated to these estimated occupancy probabilities as they are related to initial occupancy, colonization, and extinction (as well as the individual covariates within those parameters).

I would like to generate a table similar to table 7.3 (pg 204) of MacKenzie et al. 2006 using two different individual covariates that are related to colonization and extinction. Any help would be greatly appreciated.

[jhines]

PRESENCE will compute these 'recursive' occupancy rates (seasons 2-?) for each site in the data as a 'derived' parameter. If you only want those estimates for values of the covariates which are in the data, then you can just use the PRESENCE output. If you want the estimates for covariate values which aren't currently in the data, you can create 'dummy' sites (sites with all detection-history data equal to '-') which have all of the combinations of the covariate values you want and re-run the model. This would require modifying the input pao file, so I recommend creating a new project and adding the original and dummy data to it.

Alternatively, you could compute the variances using the 'Delta method'. This method is described on pg 66 of the occupancy book. Here is some sample R code I wrote for computing the variance of another 'derived' parameter:

Code: Select all

#  delta_method_example1.r
#    compute std.err of conditional psi using delta method
#    on blue ridge salamander data (model psi,p())

bta=c(0.383108,-1.052585)   #  beta's from output

#  var cov matrix of beta's from output
vc=matrix(c( 0.258663,  -0.083279,  -0.083279,  0.090511),2,2,byrow=T)
cat('var-cov(betas):\n'); print(vc)
     
psi=plogis(bta[1]); p=plogis(bta[2]); print(cbind(psi,p))

e1=.00000001; g=matrix(0,2,2)  #  g = derivative of cpsi wrt betas
g[1,1]=(plogis(bta[1]+e1)-psi)/e1;
g[2,2]=(plogis(bta[2]+e1)-p)/e1

v2= g %*% vc %*% t(g)   #  v2 = var cov matrix of real parms (psi,p)

cat('var-cov(real parms):\n');  print(v2)
   #  cpsi function to compute cond. psi from real psi,p
cpsi <- function(psi,p) {psi*(1-p)^5/(psi*(1-p)^5+1-psi)}

#  gg = derivative of cpsi wrt real parms (psi,p)
gg = (cpsi(psi+e1,p)-cpsi(psi,p))/e1
gg= c(gg,(cpsi(psi,p+e1)-cpsi(psi,p))/e1); dim(gg)=c(1,2)

vv=gg %*% v2 %*% t(gg)  #  vv = variance of cpsi

cat('var-cov(cond psi):\n'); print(vv)
cat('cond-psi=',round(cpsi(psi,p),4),' se=',round(sqrt(vv),4),'\n')

[mFidino]

Thanks!