Dear list,
I have troubles with using Bayesian logistic model with count data in bayesglm.
If I consider the following artificial data set, with a response y and a
covariate x, in a form of row data (a) and count (b) :
a<-data.frame(y=c(rep(1,10),rep(0,6)),x=c(rep(5,6),rep(4,4),rep(5,1),rep(4,5)))
a$un<-1
b<-aggregate(a$un,list(y=a$y,x=a$x),sum)
names(b)[3]<-"w"
A binomial model with glm running on a and b give the same results :
(M1=glm(y~x,family=binomial,data=a))
(M2=glm(y~x,data=b,family=binomial,weights=w))
However, with the arm and bayesglm function, the equivalent models with
non-informative prior, give different results:
library(arm)
(M3=bayesglm(y~x,data=a,family=binomial,prior.scale=Inf, prior.df=Inf)
)#M3=M2=M1
(M4=bayesglm(y~x,data=b,family=binomial,weights=w,prior.scale=Inf,
prior.df=Inf))#M4 is different
When I try a formulation with y=response rate or y=(response,failure), I get the
rights coefs, but lower standard errors.
c<-data.frame(y=c(4/9,6/7),x=c(4,5),w=c(9,7))
(M5=bayesglm(y~x,data=c,family=binomial,weights=w,prior.scale=Inf,
prior.df=Inf))
y<-matrix(c(4,6,5,1),ncol=2)
x<-matrix(c(4,5),ncol=1)
(M6=bayesglm(y~x,family=binomial,prior.scale=Inf, prior.df=Inf))
Am I missing something?
Thanks,
Edouard Chatignoux, Statistician
Health oservatory of the Paris ?le-de-France region
43 rue Beaubourg- 75003 PARIS
Tel. : 01 77 49 78 54
Fax. : 01 77 49 78 61
I'm running R 2.13.0 under windows XP
> sessionInfo()
R version 2.13.0 (2011-04-13)
Platform: i386-pc-mingw32/i386 (32-bit)
locale:
[1] LC_COLLATE=French_France.1252 LC_CTYPE=French_France.1252
[3] LC_MONETARY=French_France.1252 LC_NUMERIC=C
[5] LC_TIME=French_France.1252
attached base packages:
[1] splines stats graphics grDevices datasets grid utils
[8] methods base
other attached packages:
[1] foreign_0.8-44 arm_1.4-11 abind_1.3-0 R2WinBUGS_2.1-18
[5] coda_0.14-4 lme4_0.999375-39 Matrix_0.999375-50 lattice_0.19-26
[9] MASS_7.3-13 xtable_1.5-6 mgcv_1.7-6 ggplot2_0.8.9
[13] proto_0.3-9.2 reshape_0.8.4 plyr_1.5.2
loaded via a namespace (and not attached):
[1] nlme_3.1-101 stats4_2.13.0 tools_2.13.0
