APPENDIX C. Description of the procedure for implementing the bma_bic function for R and detailed explanation of the function.

1. Procedure

a) In R, first load Hmisc library and load the dataset that will be used to conduct the model averaging.

b) Modify the bma_bic function provided in Appendix D with the prior model weight of your choice. The prior is currently set to Occam, but the three other prior model weights used in the manuscript are available in the code. The user interested in using one of these would simply have to comment the line that contains the equation for the Occam’s prior and uncomment the line of the prior of interest (i.e., Uniform, Complex, of Kullback-Leibler).

c) Paste bma_bic function in R console.

d) Run command following the example at the end of this appendix.

2. Details about the bma_bic function for R.

Description

This function performs Bayesian Model Averaging with the BIC approximation to Bayes factors from a list of models provided by the user. The list can be all possible combinations of a set of variables, or fewer a priori selected models. The type of prior model weight must be modified by the user directly in the code. The code is set to run with the Occam’s prior of parsimony, but other priors are available as well (Uniform, Complex, Kullback-Leibler).

Usage

bma_bic(modcall,yvar,xvar,species,...)

Arguments

modcallUser-uploaded vector that lists all models that will be used to update the intercept model (basemod). See details below.

yvarThe response variable as a text string in quotations (e.g., "richness", or "sqrt(abundance)”). The text string will be used to create a statement that will update the intercept only model. If, for example, the variable name in the table is abundance, but the user is interested in fitting a linear model where the square root of the abundance is taken, then yvar should be “sqrt(abundance)”.

xvarA single string listing of all explanatory variables that are considered in the model averaging (including the intercept). The following syntax, including the quotation marks, must be used in order for the string to be correctly parsed within the function:

"Cs(Intercept,var1,var2,...)"

Note: Cs is a function of the library Hmisc that creates a vector of character strings from unquoted names.

ORThe Occam's window criteria (OR = 20 was used in the manuscript).

mydataThe name of the dataset to be used.

Details

Each model listed in modcall must be a text string of the following format without the intercept:

"update(basemod,. ~ var1 + var2 + var 3 + ...)"

The user may upload a vector of strings stored in a text file where each line corresponds to a model with variable 1, variable 2, etc. The intercept only model (basemod in the function) will automatically be updated with each new set of variables within the function.

Value

ModelsReturns the list of models

bic Returns the BIC value of the models selected for the averaging

FvalReturns the F-value of the likelihood ratio test

PvalueThe significance of the models against the null model

PostprobPosterior probability of each model

AdjR2Adjusted coefficient of determination of each model selected in the reduced set

modsizeSize of each model

varnamesList of input variables

postprobcoefPosterior probability of the model averaged coefficients

postmeanPosterior mean of the model averaged coefficients

postsdPosterior standard error of the model averaged coefficients

Reference

Hoeting, J. A., D. Madigan, A. E. Raftery, and C. T. Volinsky. 1999. Bayesian Model Averaging: A tutorial. Statistical Science 14:382-417

Authors

This function was inspired by the bic.glm and bicreg functions of the BMA package and written by Véronique St-Louis and James D. Forester.

Example of implementation

install.packages(“Hmisc”)

library(Hmisc)

#read in dataset

mybirddat<-read.csv(‘mybirddat.csv’)

#read in file that contains all the models

mymods<-read.csv(‘mymods.csv’)

#fit function using 20 as Occam’s window criteria. The explanatory variables are grass,

# shrubs, trees, urban, and forest. The response variable is the square root of the abundance.

mybma<-bma_bic(mymods,“sqrt(abundance)”,“Cs(Intercept,grass,shrubs,trees,urban,forest)”,20,mybirddat)

#To view the results, we can simply type:

mybma$postprobcoef

mybma$postmean

…