[R] Comparing MCMClogit, glm and BRUGS
francogrex
francogrex at mail.com
Sat Apr 28 23:28:12 CEST 2007
Hello,
I have two "related" questions, one about MCMClogit and the other about
BRUGS:
Given the data on nausea due to diuretic and nsaid below:
nsaid diuretic yes no
0 0 185 6527
0 1 53 1444
1 0 42 1293
1 1 25 253
A logistic regression in glm gives:
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.56335 0.07456 -47.794 < 2e-16 ***
nsaid 0.13630 0.17361 0.785 0.43242
diuretic 0.25847 0.15849 1.631 0.10293
I(nsaid * diuretic) 0.85407 0.30603 2.791 0.00526 **
But in BRUGS:
model
{
for(i in 1:N) {
yes[i] ~ dbin(p[i],no[i])
logit(p[i]) <-
beta0+beta1*nsaid[i]+beta2*diuretic[i]+beta3*(nsaid[i]*diuretic[i])
}
beta0 ~ dnorm(0,0.05)
beta1 ~ dnorm(0,0.05)
beta2 ~ dnorm(0,0.05)
beta3 ~ dnorm(0,0.05)
}
results:
> samplesStats("*")
mean sd MC_error val2.5pc median val97.5pc start sample
beta0 -3.5370 0.07481 0.001134 -3.68800 -3.5370 -3.3910 1001 10000
beta1 0.1332 0.17540 0.003035 -0.21610 0.1354 0.4663 1001 10000
beta2 0.2591 0.15710 0.002757 -0.05212 0.2608 0.5610 1001 10000
beta3 0.9142 0.30900 0.005573 0.30840 0.9176 1.5150 1001 10000
The interaction term beta3 (0.9142) is a little different from the one of
glm, why?
Using the MCMClogit (same burnin and iterations as above) from the MCMC
package gives a closer estimate to glm
Mean SD Naive SE Time-series SE
(Intercept) -3.5612 0.07678 0.0007678 0.003306
nsaid 0.1356 0.17240 0.0017240 0.007093
diuretic 0.2453 0.16045 0.0016045 0.005340
nsaid:diuretic 0.8558 0.30756 0.0030756 0.011460
But the data cannot be entered in a summary like they are above (yes and no
counts), instead they have to be entered as such:
nsaid diuretic nausea
0 0 1
0 0 1
0 0 1
0 0 1
1 0 1
etc... more 9800 rows!
Is there a way to use summary data (yes, no) with MCMClogit? THANKS
--
View this message in context: http://www.nabble.com/Comparing-MCMClogit%2C-glm-and-BRUGS-tf3663589.html#a10236835
Sent from the R help mailing list archive at Nabble.com.
More information about the R-help
mailing list