logitRegRE: Logistic regression with random effects

In erlisR/iLaplaceExamples: Extra Utilities and Examples for iLaplace

Description

The following functions give the posterior distribution and related quantities for a logistic regression model with random effects. Only the case with one observation for each group is considered, e.g. there is no replication within groups and the fixed effect is an intercept. In particular, nlpost_rebin gives the negative log-posterior for the fixed parameter (theta) and random effects (u), grU_rebin give the gradient of nlpost_rebin with respect to u and hessU_rebin gives the Hessian matrix.

Usage

nlpost_rebin(u, theta, y, n)

grU_rebin(u, theta, y, n)

hessU_rebin(u, theta, y, n)

Arguments

u	The vector of random effects. There is one random effect for each observation.
theta	The vector of fixed parameters made by a fixed effect parameter and the logartihm of the precision.
y	The response.
n	The sample size

Details

The model can be written as

Y_i \sim Bernoulli(p_i),

logit(p_i) = β + u_i,

u_i \sim N(0,τ^-1),

β \sim N(0,1),

τ \sim Gamma(1,1).

Note the τ is implemented in log.

Value

double

an n+2 dimensional vector

an (n+2)by(n+2) positive definite matrix

Examples

## Not run: 
#generate the  data
n <- 100
s2 <- 1
b <- 2
# model: logit(p_i) = b + u_i; y_i ~ dbinom(1, p_i)
set.seed(2)
u <- rnorm(n, 0, sqrt(s2))
p  <- 1/(1+exp(-(b + u)))
y <- rbinom(n,1, prob=p)

#MCMCM with JAGS (needs JAGS)
library(rjags)
dataList <- list(y = y, N=n)

# Define the model:
modelString = "
model {
for ( i in 1:N ) {
u[i] ~ dnorm(0.0, tau)
logit(p[i]) <- beta + u[i]
y[i] ~ dbern(p[i])
}
beta ~ dnorm(0.0, 1.0)
tau ~ dgamma(1.0,1.0)
ltau <- log(tau)
}
"

writeLines(modelString, con="JAGSmodel.txt" )
initsList <- list(tau=1, beta=2, u = u)

# Run the chains:
jagsModel = jags.model(file="JAGSmodel.txt",
                       data=dataList, inits=initsList,
                       n.chains=1, n.adapt=5000)
update(jagsModel, n.iter=1e+5)
codaSamples = coda.samples(jagsModel,
                           variable.names=c("beta", "ltau"),
                           n.iter = 1e+5)

plot(coda.samples)

## End(Not run)

erlisR/iLaplaceExamples documentation built on May 16, 2019, 8:48 a.m.