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.txt", 
initFile = "getInits.R", 
family = c("bernoulli", "binomial", "poisson", "gaussian"), 
spatial=NULL, spatialEffect = NULL,
reparam=NULL, prefix=NULL, priors=NULL,
brugs=length(grep("unix|linux",
					.Platform$OS.type,
					ignore.case=TRUE)))

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

spatial

For Markov Random Field models, a polygons or adjacency matrix. For Geostatistical models, a SpatialPoints objects, a matrix or data frame with columns "x" and "y", or a vector of complex numbers.

spatialEffect

spatial variable from data

reparam

vector of random effect names, subtract covariates at this level from the intercept.

prefix

string to append to object names

priors

List or vector where names refer to parameters and elements are prior distributions, for example list(SDsite="dunif(0,10)").

brugs

compatiblity with OpenBUGS, using the inprod function in place of inprod2, defaults to FALSE on windows and TRUE on unix platforms.

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

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.

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.

References

"Handling unbalanced datasets" in the "Tricks: Advanced Use of the BUGS Language" section of the bugs manual, at http://www.openbugs.net/Manuals/Tricks.html

Examples

Run this code

# NOT RUN {
library(nlme)
data(Muscle)

muscleRagged = glmmBUGS(conc ~ length, data=Muscle, 
	effects="Strip", family="gaussian",
	modelFile='model.bug', reparam='Strip')
startingValues = muscleRagged$startingValues

# }
# NOT RUN {
# run with winbugs
source("getInits.R")
require('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", 
	working.directory=getwd()
) 

# a jags example
require('R2jags')
muscleResultJags = jags(
muscleRagged$ragged, getInits, parameters.to.save = names(getInits()),
                model.file="model.bug", n.chain=3, n.iter=1000, 
                n.burnin=100, n.thin=10,
                working.directory=getwd()) 

muscleParamsJags = restoreParams(
	muscleResultJags$BUGSoutput, 
	muscleRagged$ragged) 
checkChain(muscleParamsJags)
# }
# NOT RUN {
 
data(muscleResult)

muscleParams = restoreParams(muscleResult, muscleRagged$ragged)  
summaryChain(muscleParams)
checkChain(muscleParams)



# a spatial example
# }
# NOT RUN {
library(diseasemapping)

data('popdata')
data('casedata')

model = getRates(casedata, popdata, ~age*sex)
ontario = getSMR(popdata, model, casedata)
ontario = ontario@data[,c("CSDUID","observed","logExpected")]

popDataAdjMat = spdep::poly2nb(popdata,
	row.names=as.character(popdata[["CSDUID"]]))

data('popDataAdjMat')
data('ontario')

forBugs = glmmBUGS(formula=observed + logExpected ~ 1,
  effects="CSDUID",   family="poisson", spatial=popDataAdjMat,
  spatialEffect="CSDUID",
  data=ontario)

startingValues = forBugs$startingValues

source("getInits.R")
  # find patrick's OpenBUGS executable file
  if(Sys.info()['user'] =='patrick') {	 
    obExec = system(
      "find /store/patrick/ -name OpenBUGS",
    TRUE)
    obExec = obExec[length(obExec)]
  } else {
    obExec = NULL
  }

bugsResult = bugs(forBugs$ragged, getInits, 
  parameters.to.save = names(getInits()),
    model.file="model.bug", n.chain=3, n.iter=50, n.burnin=10, 
    n.thin=2,
      program="winbugs", debug=T,working.directory=getwd())


data('ontarioResult')

ontarioParams = restoreParams(ontarioResult, forBugs$ragged)

ontarioSummary = summaryChain(ontarioParams)

# posterior probability of having 10x excess risk
postProb = apply(ontarioParams$FittedRCSDUID, 3, function(x) mean(x>log(10)) )
hist(postProb)


ontario = mergeBugsData(popdata, ontarioSummary)

spplot(ontario, "FittedRateCSDUID.mean")


ontario = mergeBugsData(ontario, postProb, newcol="postProb", by.x="CSDUID")
spplot(ontario, "postProb")



# }
# NOT RUN {
# geostatistical example

# }
# NOT RUN {
rongelap= read.table(url(
	paste("http://www.leg.ufpr.br/lib/exe/fetch.php/",
	"pessoais:paulojus:mbgbook:datasets:rongelap.txt",
	sep="")
	),header=TRUE
)

library('spdep')
coordinates(rongelap) = ~cX+cY


rongelap$logOffset = log(rongelap$time)
rongelap$site = seq(1, length(rongelap$time)) 
  
forBugs = glmmBUGS(
formula=counts + logOffset ~ 1, family="poisson",
    data=rongelap@data, effects="site", spatial=rongelap,
    priors=list(phisite="dgamma(100,1)"))
    
startingValues = forBugs$startingValues
startingValues$phi$site = 100

source("getInits.R")

rongelapResult = bugs(forBugs$ragged, getInits, 
  parameters.to.save = names(getInits()),
    model.file="model.bug", n.chain=2, n.iter=20, n.burnin=4, n.thin=2,
      program="winbugs", debug=TRUE,
      working.directory=getwd())

data('rongelapResult')

rongelapParams = restoreParams(rongelapResult, forBugs$ragged)      
      
checkChain(rongelapParams)

rongelapParams$siteGrid = CondSimuPosterior(rongelapParams, rongelap,
	gridSize=100) 

rongelapSummary=summaryChain(rongelapParams)

# plot posterior probabilities of being above average
image(rongelapSummary$siteGrid$pgt0)
# }
# NOT RUN {
# }

Run the code above in your browser using DataLab