[R-sig-ME] Does MCMCglmm allow offset?
David Atkins
datkins at u.washington.edu
Wed Dec 22 23:01:12 CET 2010
Thanks to both Gustaf and Jens for their responses. In a nutshell, I
wasn't thinking like a Bayesian! ;)
Let me unpack their proposals, to confirm *I* understand, and also so
that there is a record of why this likely works with MCMCglmm.
With Poisson regression (and it cousins) we assume a constant exposure,
that is, if we are examining number of problems related to drinking, we
assume that the time over which these problems occurred is constant
across people. If this exposure varies across people, then we typically
include the log of the exposure time as an "offset" term in the model
(effectively changing the count to a rate per specified time).
The key issue here is that the offset term is essentially an unmodelled
regression parameter, where the coefficient is set equal to 1.
What Gustaf and Jens both point out is that we can use the prior
distribution of the fixed-effects to force the coefficient of the offset
term to be 1.
Using Gustaf's code:
On 12/21/10 10:25 AM, Gustaf Granath wrote:
[snip]
>
> prior2 = list(B= list (mu = matrix(c(0,1,0,0,0),5),V = diag(5)*1e+6))
The B component of the prior is the prior for the fixed-effects; the
critical piece above is that the term for the offset is set to have a
mean of 1 (whereas everything else has mean 0). The only thing to watch
is that you know where the offset occurs in the list of variables.
> diag(prior2$B$V)[2]<-(1e-6) #change the prior variance to something small
We also want to set a very small variance around the mean of 1; the code
above targets the variance associated with the mean of 1 and sets it to
a really small amount.
At that point, we are ready to call MCMCglmm (eg, the call below).
Thanks again for the input.
cheers, Dave
> mcmc.off2 <- MCMCglmm(rapi ~ log(off30)+gender*time, data =
> rapi.df,family = "poisson", verbose = TRUE,
> random = ~ us(1 + time):id,prior=prior2)
> summary(mcmc.off2$Sol)
>
> Please correct me if Im wrong.
>
> Gustaf
More information about the R-sig-mixed-models
mailing list