Model averaging derived parameter estimates in RMark

[jlaufenb]

I recently discovered that the function model.average was changed to allow more flexibility for model averaging real parameter estimates. ��

The function model.average has been changed to a generic function. Currently it supports 2 classes: 1) list, and 2) marklist. The latter was the original model.average which has been renamed model.average.marklist and the first argument has been renamed x instead of model.list to match the standard generic function approach. The previous syntax model.average(...) will work as long as the usage does not name the first agument as in the example model.average(model.list=dipper.results,...). The list formulation (model.average.list) was created to enable a generic model averaging of estimates instead of just real parameter estimates from a mark model. It could be used with any set of estimates, model weights and estimates of precision.

Can this be used for derived estimates? So far, I have not been able to figure out how to create the vectors for N-hat and se ("Huggins" data type)that will work with model.average.list. I know that I can just as easily compute, by hand, the model average estimates using formulas from B&A and weight values. However, my candidate model set contains 20 models, which will be run on 100 random subsets of my full data set for each of 20 different subsampling scenarios (i.e., 2000 total subsets). If model.average.list is not the way to go, how can this be accomplished? I have tried every idea I can come up with to extract N-hat, se's, and weights and combine them into a single dataframe, but none have worked. Also, what is the easiest way to export dataframes, matrices, vectors, etc from R to Excel?

[dhewitt]

Also, what is the easiest way to export dataframes, matrices, vectors, etc from R to Excel?

Have a look at ?write.table (write.csv)

[jlaake]

Here is an example that is created with the edwards.eberhardt data with RMark. The first part defines and runs some Huggins models with those data. The second part shows how to use model.average.list to get a model averaged estimate of N and its se. I took code from the example with model.average.lst and simply modified it to work with this example (note: I didn't check the results. Please do so). I use a loop because it is easier for most folks to comprehend but it is possible to use sapply to do the same thing. Note that you are constructing a list for model.average.list and not a dataframe. In some cases the list could be a dataframe because each element is a vector (as in this example), but that is not always the case and in particular when you have more than one estimate per model and you want to specify the vcv matrix. If you have more than a single Nhat estimate per model then you'll need to adjust the code that I wrote to use matrices for estimate and se. See the model.average.list example to see how that is done. Let me know if this is still not clear.

With regard to the Excel question, in addition to Dave's suggestion there are also packages in R that let you read/write Excel sheets directly. See CRAN for a list of packages and their descriptions. Also, see the R Data Import/Export manual under Help/Manuals in R.

--jeff



data(edwards.eberhardt)
#
# create function that defines and runs the analyses as defined in MARK example dbf file
#
run.edwards.eberhardt=function()
{
#
# Define parameter models
#
pdotshared=list(formula=~1,share=TRUE)
ptimeshared=list(formula=~time,share=TRUE)
ptime.c=list(formula=~time+c,share=TRUE)
ptimemixtureshared=list(formula=~time+mixture,share=TRUE)
pmixture=list(formula=~mixture)
#
#
# Huggins models
#
# p=c constant over time
ee.huggins.m0 =mark(edwards.eberhardt,model="Huggins",model.parameters=list(p=pdotshared))
# p constant c constant but different; this is default model for Huggins
ee.huggins.m0.c =mark(edwards.eberhardt,model="Huggins")
# Huggins Mt
ee.huggins.Mt =mark(edwards.eberhardt,model="Huggins",model.parameters=list(p=ptimeshared),adjust=TRUE)
#
# Huggins heterogeneity models
#
# Mh2 - p different for mixture
ee.huggins.Mh2 =mark(edwards.eberhardt,model="HugHet",model.parameters=list(p=pmixture))
# Huggins Mth2 - p different for time; mixture additive
ee.huggins.Mth2.additive =mark(edwards.eberhardt,model="HugFullHet",model.parameters=list(p=ptimemixtureshared),adjust=TRUE)
#
# Return model table and list of models
#
return(collect.models() )
}
#
# fit models in mark by calling function created above
#
ee.results=run.edwards.eberhardt()


# model average Nhat value
num.models=nrow(ee.results$model.table)
estimate=vector("numeric",length=num.models)
se=vector("numeric",length=num.models)
for(i in 1:num.models)
{
# The actual model number is the row number for the model.table
model.numbers= as.numeric(row.names(ee.results$model.table))
# For each model extract the derived parameter valuess and their se
x=ee.results[[model.numbers[i]]]$results$derived
estimate[i]=x$estimate
se[i]=x$se
}
# Call model.average using the list structure which includes estimate, weight and vcv list in this case
model.average(list(estimate=estimate,weight=ee.results$model.table$weight,se=se))