Model-averaging RDHuggins

[wijw]

I ran an RDHuggins model to estimate survival and abundance of my study species. When I tried to model-average abundance estimates, I got the following error message.

Code: Select all

#Compute model-averaged estimates of abundance
estimate = matrix(0,ncol=length(robust.results[[1]]$results$derived$estimate),
                     nrow=nrow(robust.results$model.table))

#Create an empty list for vcv matrix
vcv=list(length=nrow(robust.results$model.table))

#Extract model weights
weight=robust.results$model.table$weight

#The actual model number is the row number for the model.table
model.numbers = as.numeric(row.names(robust.results$model.table))

##Loop over each model in model.table
for(i in 1:nrow(robust.results$model.table))
{
  estimate[i,]=robust.results[[model.numbers[i]]]$results$derived$estimate
  vcv[[i]]=robust.results[[model.numbers[i]]]$results$derived.vcv[[1]]
}
model.average(list(estimate=estimate,weight=weight, vcv = vcv))

Code: Select all

Error in model.average.list(list(estimate = estimate, weight = weight,  :
  dimension of one or more vcv matrices does not match dimenstion (columns) of estimate

Any help is appreciated.

[jlaake]

No way I can tell from that. It is telling you that the problem is that the dimension of the vcv does not match the columns in the estimate matrix. Have you checked to see what is wrong? The str function will help. I believe I recently posted some code to do something similar for derived parameters.

[rcscott]

There is also some very helpful text in the help files for model.average.list . I would look at the example code Jeff provides there (particularly the code for establishing the vcv matrix) and compare it to your own.

[wijw]

No way I can tell from that. It is telling you that the problem is that the dimension of the vcv does not match the columns in the estimate matrix. Have you checked to see what is wrong? The str function will help. I believe I recently posted some code to do something similar for derived parameters.

I think the issue is I have NaN in the VCV for each model. The derived N-hat is 0 for one year at one site because no surveys occurred.

[jlaake]

I would drop that estimate and the corresponding row and column in the var-cov matrix and see if that fixes it. I wouldn't think that would be the problem because the dimensions would still be fine. Do the dimensions match?

[wijw]

I would drop that estimate and the corresponding row and column in the var-cov matrix and see if that fixes it. I wouldn't think that would be the problem because the dimensions would still be fine. Do the dimensions match?

That's a good idea. I will try that.
Yes, the dimensions match.

[jlaake]

Actually the dimensions did not match the way you had the code posted.  There were 0 columns in the estimate matrix because the code should have been as shown below.  Note the [[1]] after derived to point to the first list element. If you did as below then you wouldn't get the message about dimensions not matching but you would get the message 

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

because the code was choking on the NaN when it was asking if it was less than 0 which produced NA rather than TRUE or FALSE. I have added code to check for NaN as well and the code works.  Please download the .zip file at https://drive.google.com/file/d/1-ioCyxtkjrPR9prhsnT7N7t58fjbHNlS/view?usp=drive_link and install from local zip in R. It should work for you now. 
--jeff

Code: Select all

#Compute model-averaged estimates of abundance
estimate = matrix(0,ncol=length(robust.results[[1]]$results$derived[[1]]$estimate),
                     nrow=nrow(robust.results$model.table))

#Create an empty list for vcv matrix
vcv=list(length=nrow(robust.results$model.table))

#Extract model weights
weight=robust.results$model.table$weight

#The actual model number is the row number for the model.table
model.numbers = as.numeric(row.names(robust.results$model.table))

##Loop over each model in model.table
for(i in 1:nrow(robust.results$model.table))
{
  estimate[i,]=robust.results[[model.numbers[i]]]$results$derived[[1]]$estimate
  vcv[[i]]=robust.results[[model.numbers[i]]]$results$derived.vcv[[1]]
}
model.average(list(estimate=estimate,weight=weight, vcv = vcv))