Probability of detection after k surveys (standard errors)

[carola]

Dear list members,

I have estimated the single-visit probability of detection for several amphibian species. For most of them, I have found that it changed among surveys within the same breeding season. For each species, I have computed the estimated probability of detecting the species, at least once, during the breeding season after k surveys [p* = 1 – Π(1-pk)] (formula in page 97, Mackenzie et al. (2006). Occupancy Estimation and Modeling).

How can I compute the standard error of p*? Can it be derived from the standard errors of the survey-specific probabilities of detection or by any other procedure?

Thank you!

Carola

[darryl]

Hi Carola,
You need to use something called the delta method (which is also in the book). This will require some calculus; you need to differentiate p* with respect to each of the p's, and also the portion of the variance-covariance matrix associated with all the p's.
Cheers
Darryl

[cooch]

Hi Carola,
You need to use something called the delta method (which is also in the book). This will require some calculus; you need to differentiate p* with respect to each of the p's, and also the portion of the variance-covariance matrix associated with all the p's.
Cheers
Darryl

The book Darryl is referring to of course is the MacKenzie et al. occupancy tome. The Delta method is also covered in some detail in Appendix B of the MARK book

http://www.phidot.org/software/mark/doc ... /app_2.pdf

One or the other reference should help you get started.

[jlaake]

It is not a bad thing to learn the delta method but it can be a bit of a struggle. The material in the Occupancy book on the delta method is rather limited. I suggest that you start with the Appendix in Cooch and White that Evan referred to. In a separate email I'll send you an Excel spreadsheet that shows an example of the calculation for the p(t) psi(.) model for the salamander data that are with RMark and I believe MARK. If any one else is interested in it you can email me directly (jeff.laake@noaa.gov). Evan and I have ben talking about putting some example spreadsheets on the phidot site.

You'll have to modify the EXCEL sheet for the number of occasions, the design matrix that you used and the values themselves. I've not used PRESENCE but I've been told that it does not provide the v-c matrix of the real parameters (correct me if I'm wrong on that Jim or Darryl)and that complicates the calculation somewhat.

If you had V (the v-c matrix of the p estimates) then all you need to construct is the
derivative matrix of p* with respect to p. This is simply (1-p*)/(1-pi) where pi is the
probability for occasion i and p*=1-(1-p1)*(1-p2)*...(1-pk) for k occasions. If I call that vector dp then the se of p* is t(dp)%*%V%*%dp where t(dp) is the transpose of dp and the %*% means matrix multiplication. Note that I'm assuming here model p(t) but if it was p(.) the same holds if V is a matrix with each element containing the var(p.)) although it needn't be done with matrix multiplication for the p(.) model.

Now if you don't have V then you need to construct it from the betas, the design matrix (dm) which relates the betas to the reals, link function and the v-c matrix of the betas. If you use the typical logistic link function then the derivative of the reals with respect to beta is simply the dm*p*(1-p). The dm should be the rows and columns for the p's only. Each row in the dm is for a single p. You multiply each element in row i by pi*(1-pi). Then you take that matrix and pre and post multiply the v-c matrix of the betas to get the v-c matrix for the p's. Then you can do the above step after you get V. In the EXCEL sheet I combine the derivative calculations (the chain rule for those that remember their calculus) and that may make it a little more difficult to follow. Then on the bottom of the sheet I show the calculation that is made when you have V.

For those of you that may use RMark, the calculation can be done quite simply with the following steps in R after you have run the model (called model below). You'll see that the results shown below match the value in the EXCEL sheet to about 7/8 digits. The math in R is more precise than in EXCEL.


> # extract a list containing the p estimates (ps) and their v-c matrix (V)
> plist=get.real(model,"p",vcv=T)
> ps=plist$estimates$estimate
> V=plist$vcv.real
> # compute p* = 1-(1-p1)*(1-p2)...
> pstar=1-prod(1-ps)
> # compute derivative vector of p* with respect to each p
> dp=(1-pstar)/(1-ps)
> # compute se of p*
> sepstar=sqrt(t(dp)%*%V%*%dp)
> # show values
> pstar
[1] 0.7959225
> sepstar
[,1]
[1,] 0.08419747

[jlaake]

Carola-

Send me a note so I can reply with the sheet. I had deleted your message and your email address is not posted in your profile.

Sorry to others for this message.

--jeff

[darryl]

Re: Jeff's query about VC matrix of real parameters in PRESENCE.

For the single season, predefined models you do get one for the real parameters, otherwise it's all in terms of beta-parameters.

Darryl

[cooch]

t(dp)%*%V%*%dp where t(dp) is the transpose of dp and the %*% means matrix multiplication.

Jeff is 'speakng' R here (and has kindly provided the translation). One of the bits of the learning curve in learning R is that A*B in MATLAB is the multiplication of matrices A and B, whereas the same syntax in R performs an elementwise multiplication, which is altogether a different beast. It would have been nice if R could have simply used MATLABs syntax (like Octave did, for example, but...).

On the other hand, R is faster than MATLAB for a lot of linear algebra operations, and you clearly can't argue with the price.

[jlaake]

With regard to Darryl's comment

Thanks. This must have been added in 2.2. I tried with the older version I had and didn't see it but I downloaded 2.2 and see that it is there. Good reason to upgrade.

Should be possible to add this for any specified design matrix since the formulation is the same regardless of the design matrix.

regards --jeff