Averaging derived parameters with NaN in the v-c matrix

[jCeradini]

Hi Folks,

I'm modeling small mammal data in RMark with RDHuggins.
I've model averaged the real parameters for the AICc confidence set but am having trouble averaging the derived parameter, N-hat. There are many posts on phidot regarding this topic, which I've found very helpful, but I have not found anything addressing this issue...

To average the derived parameter, I modified a loop I got off this forum:

Code: Select all

# create matrix for estimates
      estm <- matrix(0, ncol=length(test.results[[1]]$results$derived$estimate),
                                    nrow=nrow(test.results$model.table))
      # Extract models weights
      wt <- test.results$model.table$weight
      # Create empty list for vcv matrices
      vcv <- vector("list", length=nrow(test.results$model.table))
      # Loop
      for(i in 1:nrow(test.results$model.table)){
        mod.num <- as.numeric(row.names(test.results$model.table))
        x <- test.results[[mod.num[i]]]$results
        estm[i, ] <- x$derived$estimate
        vcv[[i]] <- x$derived.vcv
      }

model.average(list(estimate=estm, weight=wt, vcv=vcv))

Error in if (any(diag(x$vcv[[i]]) < 0)) { :
  missing value where TRUE/FALSE needed

I tried adding drop=FALSE argument (assuming my syntax is correct) in case there are negative variances in the v-c matrix (although my parameter estimates are not at boundaries), which returns the same error. I'm not sure if that argument can be used with anything except model.average.marklist.
model.average(list(estimate=estm, weight=wt, vcv=vcv), drop=FALSE)

I have p and c fixed to zero for many groups, which represent structural zeros (not all sites were sampled in all primary periods). This results in "NaN" in the v-c matrix for the derived parameter for the groups fixed to zero, which I'm guessing is the problem.

If the NaN's in the v-c matrix are the problem, is there an argument that can be added so they're ignored, or some other solution?

Is there an efficient way to relate negative variance values in the v-c matrix back to the parameters (my v-c matrix is fairly large)?

I appreciate any advice.
Thanks!
Joe

[jCeradini]

I talked to Jeff off list - it's a nice simple solution (that I should have thought of...). Change NaN's in the v-c matrix to zeros:

Code: Select all

model$results$derived.vcv[is.na(model$results$derived.vcv)] <- 0

After doing that, I used the loop shown in my original post to extract the values, then model.average worked fine.

Joe

[jCeradini]

Just wanted to update this post so it is compatible with RMark 2.1.12 (the for loop and replacing NaN's code above was not working for me in 2.1.12).

Model averaging abundance in Huggins robust design
- returns SE but no CI

Code: Select all

# Create matrix for estimates
estm <- matrix(0, ncol=length(ModSet[[1]]$results$derived$`N Population Size`$estimate), nrow=nrow(ModSet$model.table))

# Extract models weights
wt <- ModSet$model.table$weight

# Create empty list for vcv matrices
vcv <- vector("list", length=nrow(ModSet$model.table))

# Loop over each model
      for(i in 1:nrow(ModSet$model.table)){
        mod.num <- as.numeric(row.names(ModSet$model.table))
        x <- ModSet[[mod.num[i]]]$results
        estm[i, ] <- x$derived$`N Population Size`$estimate
        temp <- x$derived.vcv
        vcv[[i]] <- as.matrix(temp$'N Population Size')
      }

# If have NaN in vcv matrix, model.average will error.
# Can change NaN's to zeros using rapply
vcv <- rapply(vcv, f=function(x) ifelse(is.nan(x), 0, x), how="replace")

# model.average function with what extracted in loop
mod.ave <- model.average(list(estimate=estm, weight=wt, vcv=vcv))

Joe