# Bayes Factors, Posterior Samples, & DIC Computes Bayes Factors for hierarchical linear models including continuous predictors using the Savage-Dickey density ratio ```r modelrun(data, dv, dat.str, randvar.ia = NULL, corstr = NULL, nadapt = NULL, nburn = NULL, nsteps = NULL, checkconv = NULL, mcmc.save.indiv = NULL, plot.post = NULL, dic = NULL, path = NULL) ``` ## Arguments - `data`: a `data.frame` object with the data to be fitted in the long format. - `dv`: `string` indicating the dependent variable of the model. Has to be normally distributed. - `dat.str`: a `data.frame` object indicating the hierarchical structure in the model with column names "iv" and "type" and an arbitrary number of random variables as the following column names. iv indicates the name of an independent variable as in `data`, type its scale of measurement ("cont" for continuous or "cat" for categorical), and the following entries indicate whether a random effect should be modeled for this variable (1) or not (0). Continuous variables have to be entered before categorical variables. The name for the random variable(s) has to be the same as in `data`. A categorical variable with n levels is entered as n - 1 simple codes into the model with the first level of the variable being the reference category. - `randvar.ia`: a `list` containing n `matrix` objects with n being the number of random variables. In each `matrix` the lower triangle can be used to declare the respective two-way interaction as random within the specific random variable. The row- and column- ordering of independent variables is the same as in `dat.str`. When not specified, interactions are only modeled as fixed effects by default. - `corstr`: a `list` containing n `matrix` objects with n being the number of random variables. In each `matrix` the lower triangle can be used to assign correlations between predictors (including the intercept) for each random effect. The first row and column in each `matrix` object represents the intercept. The following rows and columns represent the independent variables with the same ordering as in `dat.str`. When not specified, no correlations are modeled by default. - `nadapt`: number of MCMC steps to adapt the sampler (2000 by default). - `nburn`: number of MCMC steps to burn in the sampler (2000 by default). - `nsteps`: number of saved MCMC steps in all chains (100'000 by default). - `checkconv`: indicates that convergence statistics of the main model parameters should be returned in the console and that figures of the chains should be plotted when set to 1 (0 by default). - `mcmc.save.indiv`: indicates that the chains should be saved in a `data.frame` object when set to 1 (0 by default). - `plot.post`: indicates that the 95 percent highest-density interval of the posterior of the group parameters should be plotted as a figure with the corresponding Bayes Factors when set to 1 (0 by default). - `dic`: indicates that the deviation information criterion (Spiegelhalter, Best, Carlin, & Linde, 2002) should be computed for a given model when set to 1 (0 by default). - `path`: defines directory where model is saved as .txt file and model name. Is set to file.path(tempdir(), "model.txt") by default. ## Returns returns a list with components: * `bf`: a `data.frame` object with the Bayes Factor estimates of the group parameters (aka fixed effects). * `mcmcdf`: a `data.frame` object with the saved MCMC chains. * `dic`: DIC of the fitted model. ## Details The argument `corstr` can be used to model correlations between (a) pairs of predictors and (b) more than two predictors. When both is done within the same random variable, a predictor can only appear in (a) or (b). `modelrun` z-standardizes the dependent variable and the continuous independent variables. To obtain the posteriors in the original scale they have to be retransformed. ### Savage Dickey Bayes Factors are computed with the Savage-Dickey density ratio. We use the normal approximation (e.g., Wetzels, Raaijmakers, Jakab, & Wagenmakers, 2009) to estimate the density of the posterior. ## Examples ```r data(bayesrsdata) #load data ## ----------------------------------------------------------------- ## Example 1: Estimation of Bayes Factors from a continuous ## independent variable (IV) with random slopes ## - repeated measures for each participant ## - continuous variable with 5 values: x.time ## ------------------------------------------------------------------ ## JAGS Sampler Settings # ----------------- # nr of adaptation, burn-in, and saved mcmc steps only for exemplary use nadapt = 2000 # number of adaptation steps nburn = 2000 # number of burn-in samples mcmcstep = 100000 # number of saved mcmc samples, min. should be 100'000 # Define model structure; dat.str <- data.frame(iv = c("x.time"), type = c("cont"), subject = c(1)) # name of random variable (here 'subject') needs to match data frame # Run modelrun function out <- modelrun(data = bayesrsdata, dv = "y", dat.str = dat.str, nadapt = nadapt, nburn = nburn, nsteps = mcmcstep, checkconv = 0) # Obtain Bayes factor bf <- out[[1]] bf ## ----------------------------------------------------------------- ## Example 2: Estimation of Bayes Factors from a continuous ## independent variable with random slopes that ## are correlated with the random slopes of a categorical variable. ## - Repeated measures for each participant ## - a continuous IV with 5 values: x.time ## - a categorical variable with 2 levels: x.domain ## ------------------------------------------------------------------ ## JAGS Sampler Settings # nr of adaptation, burn-in, and saved mcmc steps only for exemplary use # ----------------- nadapt = 2000 # number of adaptation steps nburn = 2000 # number of burn-in samples mcmcstep = 100000 # number of saved mcmc samples, min. should be 100'000 # Define model structure; # order of IVs: continuous variable(s) needs to go first dat.str <- data.frame(iv = c("x.time", "x.domain"), type = c("cont", "cat"), subject = c(1,1)) # name of random variable (here 'subject') needs to match data frame # Define random effect structure on interaction for each random variable ias.subject <- matrix(0, nrow=nrow(dat.str), ncol = nrow(dat.str)) ias.subject[c(2)] <- 1 randvar.ia <- list(ias.subject) # Define correlation structure between predictors within a random variable cor.subject <- matrix(0, nrow=nrow(dat.str)+1, ncol = nrow(dat.str)+1) cor.subject[c(2,3,6)] <- 1 corstr <- list(cor.subject) # Run modelrun function out <- modelrun(data = bayesrsdata, dv = "y", dat.str = dat.str, randvar.ia = randvar.ia, nadapt = nadapt, nburn = nburn, nsteps = mcmcstep, checkconv = 0, mcmc.save.indiv = 1, corstr = corstr) # Obtain Bayes factors for continous main effect, # categorical main effect, and their interaction bf <- out[[1]] bf ``` ## References Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & van der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(4), 583. Wetzels, R., Raaijmakers, J. G. W., Jakab, E., & Wagenmakers, E.-J. (2009). How to quantify support for and against the null hypothesis: A flexible WinBUGS implementation of a default Bayesian t test. Psychonomic Bulletin & Review, 16(4), 752-760. https://doi.org/10.3758/PBR.16.4.752 ## Author(s) Thalmann, M., Niklaus, M. Part of this package uses code from John Kruschke.