[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