I ran a popan model in RMark, with 2 grouping variables (Var1 and Var2, two levels for each) and 8 occasions. I can get the N (by the 2 grouping variables and the 8 occasions) and the N by occasion (summed across the grouping variables) using popan.derived. However, I want to calculate N by occasion and by Var1 (=summed across Var2). I used deltamethod.special to try and calculate the correct SEs.
A quick test of the code (with the Var1 loop removed, so that there's only by-occasion calculation) against the N by occasion that is provided by popan.derived showed that my code works perfectly for Occasion 1, but fails for further occasions (=the mean estimates are correct, but the SEs are not). I'm guessing that I'm grabbing the wrong bits from the vcv matrix, but can't figure out how to fix that.
The variable data1.results is a list of MARK models; data1.proc is the processed dataset.
- Code: Select all
derived <- popan.derived(x = data1.proc, model = data1.chat.results)
# calculate N by HW and occasion; SEs are calculated using the Delta method
out <- list()
k <- 1
for(var1 in c("A", "B")){
for(i in unique(derived$N$Occasion)){
id <- which(derived$N$Var1 == var1 & derived$N$Occasion == i)
mat <- derived$N.vcv[id, id]
se <- deltamethod.special("sum", derived$N[id,]$N, mat)
out[[k]] <- data.frame(Var1 = var1, Occasion = i, derived$N = sum(derived$N[id,]$N), SE = se)
k <- k + 1 }}
out <- ldply(out, data.frame)
Any help would be appreciated.