Gibbs Sampler for Hierarchical Linear Model with Normal Heterogeneity
rhierLinearModel implements a Gibbs Sampler for hierarchical linear models with a normal prior.
rhierLinearModel(Data, Prior, Mcmc)
Data |
list(regdata, Z) |
Prior |
list(Deltabar, A, nu.e, ssq, nu, V) |
Mcmc |
list(R, keep, nprint) |
nreg regression equations with nvar X variables in each equation
y_i = X_iβ_i + e_i with e_i ~ N(0, τ_i)
τ_i ~ nu.e*ssq_i/χ^2_{nu.e} where τ_i is the variance of e_i
β_i ~ N(ZΔ[i,], V_{β})
Note: ZΔ is the matrix Z * Δ and [i,] refers to ith row of this product
vec(Δ) given V_{β} ~ N(vec(Deltabar), V_{β}(x) A^{-1})
V_{β} ~ IW(nu,V)
Delta, Deltabar are nz x nvar; A is nz x nz; V_{β} is nvar x nvar.
Note: if you don't have any Z variables, omit Z in the Data argument and
a vector of ones will be inserted; the matrix Δ will be 1 x nvar and should
be interpreted as the mean of all unit βs.
Data = list(regdata, Z) [Z optional]
regdata: |
list of lists with X and y matrices for each of nreg=length(regdata) regressions |
regdata[[i]]$X: |
n_i x nvar design matrix for equation i |
regdata[[i]]$y: |
n_i x 1 vector of observations for equation i |
Z: |
nreg x nz matrix of unit characteristics (def: vector of ones) |
Prior = list(Deltabar, A, nu.e, ssq, nu, V) [optional]
Deltabar: |
nz x nvar matrix of prior means (def: 0) |
A: |
nz x nz matrix for prior precision (def: 0.01I) |
nu.e: |
d.f. parameter for regression error variance prior (def: 3) |
ssq: |
scale parameter for regression error var prior (def: var(y_i)) |
nu: |
d.f. parameter for Vbeta prior (def: nvar+3) |
V: |
Scale location matrix for Vbeta prior (def: nu*I) |
Mcmc = list(R, keep, nprint) [only R required]
R: |
number of MCMC draws |
keep: |
MCMC thinning parm -- keep every keepth draw (def: 1) |
nprint: |
print the estimated time remaining for every nprint'th draw (def: 100, set to 0 for no print)
|
A list containing:
betadraw |
nreg x nvar x R/keep array of individual regression coef draws |
taudraw |
R/keep x nreg matrix of error variance draws |
Deltadraw |
R/keep x nz*nvar matrix of Deltadraws |
Vbetadraw |
R/keep x nvar*nvar matrix of Vbeta draws |
Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.
For further discussion, see Chapter 3, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.
http://www.perossi.org/home/bsm-1
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
set.seed(66)
nreg = 100
nobs = 100
nvar = 3
Vbeta = matrix(c(1, 0.5, 0, 0.5, 2, 0.7, 0, 0.7, 1), ncol=3)
Z = cbind(c(rep(1,nreg)), 3*runif(nreg))
Z[,2] = Z[,2] - mean(Z[,2])
nz = ncol(Z)
Delta = matrix(c(1,-1,2,0,1,0), ncol=2)
Delta = t(Delta) # first row of Delta is means of betas
Beta = matrix(rnorm(nreg*nvar),nrow=nreg)%*%chol(Vbeta) + Z%*%Delta
tau = 0.1
iota = c(rep(1,nobs))
regdata = NULL
for (reg in 1:nreg) {
X = cbind(iota, matrix(runif(nobs*(nvar-1)),ncol=(nvar-1)))
y = X%*%Beta[reg,] + sqrt(tau)*rnorm(nobs)
regdata[[reg]] = list(y=y, X=X)
}
Data1 = list(regdata=regdata, Z=Z)
Mcmc1 = list(R=R, keep=1)
out = rhierLinearModel(Data=Data1, Mcmc=Mcmc1)
cat("Summary of Delta draws", fill=TRUE)
summary(out$Deltadraw, tvalues=as.vector(Delta))
cat("Summary of Vbeta draws", fill=TRUE)
summary(out$Vbetadraw, tvalues=as.vector(Vbeta[upper.tri(Vbeta,diag=TRUE)]))
## plotting examples
if(0){
plot(out$betadraw)
plot(out$Deltadraw)
}Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.