Description Usage Arguments Details Value Warning Note Author(s) References See Also Examples
Creates ragged arrays, writes a model file, and generates sensible starting estimates.
1 2 3 4 5 6 7 8 |
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 |
brugs |
compatiblity with OpenBUGS, using the inprod function in place of inprod2, defaults to FALSE on windows and TRUE on unix platforms. |
Consider the following model, where Y_ijk is the number of absences from individual k from class j in school k.
Y_ijk ~ Poisson(mu_i)
log(mu_i) = intercept + age_ijk beta + classSize_ij alpha + schoolCategory_i gamma + U_i + V_ij
U_i ~ N(0, sigma^2)
V_ij ~ 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 ~ N(U_i + classSize_ij alpha , nu^2)
U_i ~ N(intercept + 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.
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
.
You are strongly encouraged to modify the model file and the intial value function file prior to using them.
glmmBUGS uses the inprod2
function, which isn't implimented in OpenBugs, the model file will have to be modified for use with OpenBUGS.
Patrick Brown, patrick.brown@utoronto.ca
"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
winBugsRaggedArray
, glmmPQLstrings
, writeBugsModel
, getStartingValues
, startingFunction
,bugs
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 |
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)
## End(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")
## End(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)
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.