[R-sig-ME] Does MCMCglmm allow offset?

[Jens Åström]

Setting the prior of "off30" manually to 1 produce results that appear
very similar to glmer with "off30" as an offset.

It would be interesting to know if this is the intended way of
specifying an offset with MCMCglmm.


########################
prior1 = list(B= list (mu = matrix(c(0,0,0,1,0),5),V = diag(5)*(10)))
diag(prior1$B$V)[4]<-1e-9

mcmc.manual.off <- MCMCglmm(rapi ~ gender*time + log(off30), data =
rapi.df,family = "poisson", verbose = TRUE, random = ~ us(1 +
time):id,prior=prior1)


summary(mcmc.manual.off)					
#########################

/Jens



------------------------------

Message: 6
Date: Mon, 20 Dec 2010 14:20:35 -0800
From: David Atkins <datkins at u.washington.edu>
To: "r-sig-mixed-models at r-project.org"
	<r-sig-mixed-models at r-project.org>
Subject: [R-sig-ME] Does MCMCglmm allow offset?
Message-ID: <4D0FD6B3.2070300 at u.washington.edu>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed


Hi all--

I am wondering whether MCMCglmm will take an offset?  Specifically, we
have a study with number of drinks as the outcome, but we also have an
estimate of over how many hours those drinks were consumed.  Thus, might
make sense to model it as a rate of drinks per hour.

As far as I can tell, glmer does allow an offset, but it seems like
maybe MCMCglmm does not.  A search of MCMCglmm's course notes and JSS
paper does not turn up any discussion.

Below is a real dataset, with a "toy" offset term that would seem to
support this.

Any input appreciated.

cheers, Dave

[sessionInfo() at very end]

### MCMCglmm and offset
#
### read data
rapi.df <-
read.csv("http://depts.washington.edu/cshrb/newweb/stats%20documents/RAPI.Final.csv")
head(rapi.df)

### "rapi" is the rutgers alcohol problems index and was measured
### over previous 30 days -- thus, constant exposure, but we can still
### use it as toy example
#
### over-dispersion term
rapi.df$over <- 1:nrow(rapi.df)

### make up "offset" term
rapi.df$off30 <- rpois(n = nrow(rapi.df), lambda = 30)

### GLM (no random-effects) just as baseline comparator
glm.nooff <- glm(rapi ~ gender*time,
						data = rapi.df, family = "poisson")
glm.off <- glm(rapi ~ gender*time + offset(log(off30)),
						data = rapi.df, family = "poisson")
summary(glm.nooff)						
summary(glm.off)						

### over-dispersed Poisson GLMM in glmer
library(lme4)

### without offset
glmer.nooff <- glmer(rapi ~ gender*time + (time | id) + (1 | over),
						data = rapi.df, verbose = TRUE, family = "poisson")

### with offset
glmer.off <- glmer(rapi ~ gender*time + offset(log(off30)) + (time | id)
+ (1 | over),
						data = rapi.df, verbose = TRUE, family = "poisson")
print(glmer.nooff, corr = FALSE)	
print(glmer.off, corr = FALSE)	

### same models using MCMCglmm
library(MCMCglmm)

mcmc.nooff <- MCMCglmm(rapi ~ gender*time, data = rapi.df,
								family = "poisson", verbose = TRUE,
								random = ~ us(1 + time):id)
mcmc.off <- MCMCglmm(rapi ~ gender*time + offset(log(off30)), data =
rapi.df,
								family = "poisson", verbose = TRUE,
								random = ~ us(1 + time):id)
summary(mcmc.nooff)					
summary(mcmc.off)	


 > sessionInfo()
R version 2.12.1 (2010-12-16)
Platform: i386-apple-darwin9.8.0/i386 (32-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods
[7] base

other attached packages:
[1] MCMCglmm_2.10      corpcor_1.5.7      ape_2.6-2
[4] coda_0.14-2        tensorA_0.36       lme4_0.999375-37
[7] Matrix_0.999375-47 lattice_0.19-16

loaded via a namespace (and not attached):
[1] gee_4.13-16   grid_2.12.1   nlme_3.1-97   stats4_2.12.1
[5] tools_2.12.1

-- Dave Atkins, PhD Research Associate Professor Department of
Psychiatry and Behavioral Science University of Washington
datkins at u.washington.edu Center for the Study of Health and Risk
Behaviors (CSHRB) 1100 NE 45th Street, Suite 300 Seattle, WA 98105
206-616-3879 http://depts.washington.edu/cshrb/ (Mon-Wed) Center for
Healthcare Improvement, for Addictions, Mental Illness, Medically
Vulnerable Populations (CHAMMP) 325 9th Avenue, 2HH-15 Box 359911
Seattle, WA 98104 http://www.chammp.org (Thurs)