delta method using RMark

[megpetrie]

Hello,
I need to do exactly what is outlined in the delta method appendix B in the Gentle Intro (B-12 thru B-14) for a large number of parameters. That is, I need to calculate a standard error for the product of three consecutive survival parameters, for 3-4 populations in each of three years. I was hoping there was a way to do it in RMark but I wasn't clear on how to set it up based on the delta method section in the RMark appendix.
I know it should be something like:
deltamethod=(~(x1)*(x2)*(x3),mean=beta$estimate,mymodel$results$beta.vcv) but I don't know how to specify the right parameters in the model to use for the transformation or for the variance-covariance matrix.
Does anyone have some code they could share so that I know I'm doing it right?


Thanks,
Megan

[jlaake]

You should be using get.real which provides the real parameter estimates and their v-c matrix. Also, did you look at deltamethod.special? It handles a product of estimates as you need. Below I modified the example in the help for deltamethod.special so it does 2 products of 3 estimates.

Code: Select all

data(dipper)
# run model
  mod=mark(dipper,model.parameters=list(Phi=list(formula=~time)))
# Get estimates of Phi and their v-c matrix
  rr=get.real(mod,"Phi",se=TRUE,vcv=TRUE)
# Compute se of the product of the first 3 survivals
  deltamethod.special("prod",rr$estimates$estimate[1:3],rr$vcv.real[1:3,1:3])
# Compute se of the product of the last 3 survivals
  deltamethod.special("prod",rr$estimates$estimate[4:6],rr$vcv.real[4:6,4:6])

You need to use get.real to get the vector of real parameter estimates and their v-c matrix and then you just need to subset it to meet your needs.

What the above will NOT do is to give you the v-c matrix of sets of products. Products are extremely easy with the delta method because the derivative with respect to the real parameter is just the product divided by that real parameter -- which is the product of the other 2 parameters.

Code: Select all

#Below I wrote a function in R that will compute the v-c matrix for an abitrary set of products of all the same length. x is the vector of estimates,
#vcv is the variance-covariance matrix of the estimates and indices is a matrix of indices in x for how you want to combine them. The number of rows
#in indices is the number of products and the number of columns is the number of estimates in each product.

dm.prod=function(x,vcv,indices)
{
  D=matrix(0,nrow=nrow(indices),ncol=prod(dim(indices)))
  for(i in 1:nrow(indices))
  {
    estimates=x[indices[i,]]
   D[i,indices[i,]]=prod(estimates)/estimates
  }
  return(D%*%vcv%*%t(D))
}

#Here is an example with the dipper data from above

# create 2 products with first 3 and last 3
indices=matrix(1:6,ncol=3,nrow=2,byrow=T)
# show indices
indices
# compute v-c matrix of the 2 products
vcv=dm.prod(rr$estimates$estimate,rr$vcv,indices=indices)
# show that std errors match the results above
sqrt(diag(vcv))

#You can expand to any number of products of any number of estimates by specifying a different index matrix.
# For example here is one that creates 3 products of each pair of estimates
indices=matrix(1:6,ncol=2,nrow=3,byrow=T)
indices
# Note indices do not have to be in order, just easier to do for the example
vcv=dm.prod(rr$estimates$estimate,rr$vcv,indices=indices)
vcv

I hope this helps. --jeff

[megpetrie]

This is extremely helpful thanks Jeff.
I've started computing my standard errors but they seem really low. Should they be lower than any of the individual parameter standard errors? I'm thinking there's an issue with estimating standard errors when there is fixed detection rates involved. Could that be the issue?

Megan

[jlaake]

If you are multiplying values less than one together, the product is smaller than the individual values and the se will likely be as well. What you should look at is the cv=se/estimate. Also, with covariances it is hard to predict the direction of the se. For example the sum of 2 random variables that are perfectly negatively correlated (R=-1) will be 0.

--jeff

[megpetrie]

OK thanks.
One more question. Is it cause for alarm if the se is larger than the (product) estimate? I guess I'm just wondering what error this new standard error is accounting for. Is it a composite of the se's from all three estimates? It just doesn't seem to correlate with the individual estimate se's. Basically, I can't see where some of the really high se's are coming from.

Megan

[jlaake]

It is certainly not want you want but is possible if your individual estimates are highly positively correlated. If they were independent (covariances=0) then the delta method approximation to the cv is
cv(x1*x2*x3)=sqrt(cv1^2+cv2^2+cv3^2) where cvi is the cv of the ith estimate. But it is unlikely that the covariances are 0. Send me an example offlist of the 3 estimates and their v-c matrix and I'll look at it.

--jeff