From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida  

 

you are kindly asked to include the complete citation if you used this material in a publication

Code 5.15 Synthetic data from logistic model in R

====================================================

set.seed(13979)
nobs <- 5000

x1 <- rbinom(nobs, size = 1, 0.6)
x2 <- runif(nobs)

 

xb <- 2 + 0.75*x1 - 5*x2

exb <- 1/(1+exp(-xb))
by <- rbinom(nobs, size = 1, prob = exb)

 

logitmod <- data.frame(by, x1, x2)

====================================================

 

Code 5.16  - Logistic model using R
====================================================

library(MCMCpack)
 

myL <- MCMClogit(by ~ x1 + x2,
                 burnin = 5000,
                 mcmc = 10000,
                 data = logitmod)

 

summary(myL)

====================================================

 

Output on screen:

Iterations = 5001:15000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 10000

 

1. Empirical mean and standard deviation for each variable,

    plus standard error of the mean:

                                   Mean              SD         Naive SE         Time-series SE

(Intercept)                1.9408      0.08314        0.0008314                  0.002651

x1                             0.8506      0.07183        0.0007183                  0.002356

x2                            -5.0118      0.14463        0.0014463                  0.004709

2. Quantiles for each variable:

                                2.5%         25%          50%          75%         97.5%

(Intercept)            1.7772        1.884        1.943       1.9970          2.096

x1                         0.7116        0.803        0.850       0.8963          1.000

x2                        -5.2904      -5.109       -5.016      -4.9131        -4.719