glmmBUGS (version 1.0)

glmmBUGS: A function to run Generalised Linear Mixed Models in Bugs

Description

Creates ragged arrays, writes a model file, and generates sensible starting estimates.

Usage

glmmBUGS(formula, data, effects, modelFile = "model.bug", initFile = "getInits.R", family = c("bernoulli", "binomial", "poisson", "gaussian"))

Arguments

formula

A formula for the fixed effects portion of the model.

data

A data frame containing the response, covariates, and group membership.

effects

A vector of character strings containing the grouping levels, from most general to most specific

modelFile

File for saving the bugs model

initFile

File for saving the function for generating initial values

family

distribution of responses

Value

• Returns a list with the ragged array, from winBugsRaggedArray, and the list of starting values from getStartingValues. Writes a model file and an initial value function. Note that the initial value function in initFile will look for an object called startingValues, which does not exist as this is part of a list. Either create startingValues <- result$startingValues or edit initFile.

Warning

You are strongly encouraged to modify the model file and the intial value function file prior to using them.

Details

glmmBUGS calls the following functions: getDesignMatrix{to convert factors and interactions to indicator variables and find which covariates apply at which levels} winBugsRaggedArray{ to prepare the ragged array} glmmPQLstrings{ estimate starting values} writeBugsModel{ to create a model file} getStartingValues{to extract starting values from the glmmPQL result} startingFunction{to write a function to generate random starting values} Type glmmBUGS on the R command line to see the source code, it provides a good summary of the roles of the various functions in the glmmBUGS package. Consider the following model, where $Y_{ijk}$ is the number of absences from individual k from class j in school k. $$Y_{ijk} \sim Poisson(\mu_i)$$ $$\log(\mu_i) = \delta age_{ijk} \beta + classSize_{ij} \alpha + schoolCategory_i \gamma + U_i + V_{ij}$$ $$U_i \sim N(0, \sigma^2)$$ $$V_{ij} \sim N(0, \nu^2)$$ Here there are covariates which apply to each of the three levels, and random effects at the school and class level. If data is a data frame with one line per individual, the following would impliment this model: glmmBUGS(data, effects=c("school","class"), covariates = list(school="schoolCategory", class="classSize", observations="age"), observations = "absences"), family="poisson") To aid in convergence, the bugs model is actually the following: $$\log(\mu_i) = age_{ijk} \beta + V_{ij}$$ $$V_{ij} \sim N(U_i + classSize_{ij} \alpha , \nu^2)$$ $$U_i \sim N(\delta + schoolCategory_i \gamma, \sigma^2)$$ and the funciton restoreParams subtracts the means from the random effects to restore the original set of equations.

References

"Handling unbalanced datasets" in the "Tricks: Advanced Use of the BUGS Language" section of the bugs manual, at http://mathstat.helsinki.fi/openbugs/data/Docu/Tricks.html#HandlingUnbalancedDatasets

See Also

winBugsRaggedArray, glmmPQLstrings , writeBugsModel, getStartingValues, startingFunction,bugs

Examples

library(nlme)
data(Muscle)
muscleRagged = glmmBUGS(conc ~ length, data=Muscle, effects="Strip", family="gaussian")
startingValues = muscleRagged$startingValues
source("getInits.R")
# library(R2WinBUGS)
# muscleResult = bugs(muscleRagged$ragged, getInits, parameters.to.save = names(getInits()), model.file="model.bug", n.chain=3, n.iter=1000, n.burnin=100, n.thin=10, program="winbugs")   
data(muscleResult)
muscleParams = restoreParams(muscleResult, muscleRagged$ragged)  
summaryChain(muscleParams)
checkChain(muscleParams)