top of page

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 8.5 Simulated random intercept binary logistic data

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

N <- 4000                                                                 # 20 groups, each with 200 observations
NGroups <- 20

â€‹

x1 <- runif(N)
x2 <- runif(N)

Groups <- rep(1:20, each = 200)

â€‹

a <- rnorm(NGroups, mean = 0, sd = 0.5)
eta <- 1 + 0.2 * x1 - 0.75 * x2 + a[Groups]

mu <- 1/(1+exp(-eta))
y <- rbinom(N, prob=mu, size=1)

â€‹

logitr <- data.frame(
y = y,
x1 = x1,
x2 = x2,
Groups = Groups,
RE = a[Groups]
)

â€‹

â€‹

â€‹

â€‹

â€‹

Code 8.6 - Bayesian Random intercept binary model in R

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

library(MCMCglmm)

â€‹

BayLogitRI <- MCMCglmm(y ~ x1 + x2, random= ~Groups,
family="ordinal", data=logitr,
verbose=FALSE, burnin=10000, nitt=20000)

â€‹

summary(BayLogitRI)

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

Anchor 1

Output on screen:

â€‹

Iterations = 10001:19991

Thinning interval = 10

Sample size = 1000

DIC: 4705.982

G-structure:   ~  Groups

post.mean       l-95% CI           u-95% CI          eff.samp

Groups           0.2735        0.06902               0.6496               6.902

R-structure:     ~units

post.mean          l-95%          CI u-95%            CI eff.samp

units               2.128         0.7019                6.548                      3.658

Location effects: y ~ x1 + x2

post.mean       l-95% CI      u-95% CI       eff.samp           pMCMC

(Intercept)          0.9364           0.5490           1.5244            4.614              <0.001 ***

x1                       0.2526         -0.0136            0.5378         53.777                0.054 .

x2                     -0.8450          -1.3951          -0.5057            4.543              <0.001 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

bottom of page