[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