MCMC Estimation for Mixed Effects Model

Fits a mixed effects model via MCMC. The outcome can be normally distributed or ordinal (Goldstein, 2011; Goldstein, Carpenter, Kenward & Levin, 2009).

ml_mcmc( formula, data, iter=3000, burnin=500, print_iter=100, outcome="normal",
     nu0=NULL, s0=1, psi_nu0_list=NULL, psi_S0_list=NULL, inits_lme4=FALSE,
     thresh_fac=5.8, ridge=1e-5)

# S3 method for ml_mcmc
summary(object, digits=4, file=NULL, ...)

# S3 method for ml_mcmc
plot(x, ask=TRUE, ...)

# S3 method for ml_mcmc
coef(object, ...)

# S3 method for ml_mcmc
vcov(object, ...)

ml_mcmc_fit(y, X, Z_list, beta, Psi_list, sigma2, alpha, u_list, idcluster_list,
    onlyintercept_list, ncluster_list, sigma2_nu0, sigma2_sigma2_0, psi_nu0_list,
    psi_S0_list, est_sigma2, est_probit, parameter_index, est_parameter, npar, iter,
    save_iter, verbose=TRUE, print_iter=500, parnames0=NULL, K=9999, est_thresh=FALSE,
    thresh_fac=5.8, ridge=1e-5, parm_summary=TRUE )

## exported Rcpp functions
miceadds_rcpp_ml_mcmc_sample_beta(xtx_inv, X, Z_list, y, u_list, idcluster_list, sigma2,
     onlyintercept_list, NR, ridge)
miceadds_rcpp_ml_mcmc_sample_u(X, beta, Z_list, y, ztz_list, idcluster_list,
    ncluster_list, sigma2, Psi_list, onlyintercept_list, NR, u0_list, ridge)
miceadds_rcpp_ml_mcmc_sample_psi(u_list, nu0_list, S0_list, NR, ridge)
miceadds_rcpp_ml_mcmc_sample_sigma2(y, X, beta, Z_list, u_list, idcluster_list,
     onlyintercept_list, nu0, sigma2_0, NR, ridge)
miceadds_rcpp_ml_mcmc_sample_latent_probit(X, beta, Z_list, u_list, idcluster_list, NR,
     y_int, alpha, minval, maxval)
miceadds_rcpp_ml_mcmc_sample_thresholds(X, beta, Z_list, u_list, idcluster_list, NR, K,
     alpha, sd_proposal, y_int)
miceadds_rcpp_ml_mcmc_predict_fixed_random(X, beta, Z_list, u_list, idcluster_list, NR)
miceadds_rcpp_ml_mcmc_predict_random_list(Z_list, u_list, idcluster_list, NR, N)
miceadds_rcpp_ml_mcmc_predict_random(Z, u, idcluster)
miceadds_rcpp_ml_mcmc_predict_fixed(X, beta)
miceadds_rcpp_ml_mcmc_subtract_fixed(y, X, beta)
miceadds_rcpp_ml_mcmc_subtract_random(y, Z, u, idcluster, onlyintercept)
miceadds_rcpp_ml_mcmc_compute_ztz(Z, idcluster, ncluster)
miceadds_rcpp_ml_mcmc_compute_xtx(X)
miceadds_rcpp_ml_mcmc_probit_category_prob(y_int, alpha, mu1, use_log)
miceadds_rcpp_pnorm(x, mu, sigma)
miceadds_rcpp_qnorm(x, mu, sigma)
miceadds_rcpp_rtnorm(mu, sigma, lower, upper)

Arguments

formula

An R formula in lme4-like specification

data

Data frame

iter

Number of iterations

burnin

Number of burnin iterations

print_iter

Integer indicating that every print_iterth iteration progress should be displayed

outcome

Outcome distribution: "normal" or "probit"

nu0

Prior sample size

s0

Prior guess for variance

inits_lme4

Logical indicating whether initial values should be obtained from fitting the model in the lme4 package

thresh_fac

Factor for proposal variance for estimating thresholds which is determined as thresh_fac\(/N\) (\(5.8/N\) as default).

ridge

Ridge parameter for covariance matrices in sampling

object

Object of class ml_mcmc

digits

Number of digits after decimal used for printing

file

Optional file name

...

Further arguments to be passed

x

Object of class ml_mcmc

ask

Logical indicating whether display of the next plot should be requested via clicking

y

Outcome vector

X

Design matrix fixed effects

Z_list

Design matrices random effects

beta

Initial vector fixed coefficients

Psi_list

Initial covariance matrices random effects

sigma2

Initial residual covariance matrix

alpha

Vector of thresholds

u_list

List with initial values for random effects

idcluster_list

List with cluster identifiers for random effects

onlyintercept_list

List of logicals indicating whether only random intercepts are used for a corresponding random effect

ncluster_list

List containing number of clusters for each random effect

sigma2_nu0

Prior sample size residual variance

sigma2_sigma2_0

Prior guess residual variance

psi_nu0_list

List of prior sample sizes for covariance matrices of random effects

psi_S0_list

List of prior guesses for covariance matrices of random effects

est_sigma2

Logical indicating whether residual variance should be estimated

est_probit

Logical indicating whether probit model for ordinal outcomes should be estimated

parameter_index

List containing integers for saving parameters

est_parameter

List of logicals indicating which parameter type should be estimated

npar

Number of parameters

save_iter

Vector indicating which iterations should be used

verbose

Logical indicating whether progress should be displayed

parnames0

Optional vector of parameter names

K

Number of categories

est_thresh

Logical indicating whether thresholds should be estimated

parm_summary

Logical indicating whether a parameter summary table should be computed

xtx_inv

Matrix

NR

Integer

ztz_list

List containing design matrices for random effects

u0_list

List containing random effects

nu0_list

List with prior sample sizes

S0_list

List with prior guesses

sigma2_0

Numeric

y_int

Integer vector

minval

Numeric

maxval

Numeric

sd_proposal

Numeric vector

N

Integer

Z

Matrix

u

Matrix containing random effects

idcluster

Integer vector

onlyintercept

Logical

ncluster

Integer

mu1

Vector

use_log

Logical

mu

Vector

sigma

Numeric

lower

Vector

upper

Vector

Details

Fits a linear mixed effects model \(y=\bm{X}\bm{beta}+ \bm{Z}\bm{u}+e\) with MCMC sampling. In case of ordinal data, the ordinal variable \(y\) is replaced by an underlying latent normally distributed variable \(y^\ast\) and the residual variance is fixed to 1.

Value

List with following entries (selection)

sampled_values

Sampled values

par_summary

Parameter summary

References

Goldstein, H. (2011). Multilevel statistical models. New York: Wiley. doi:10.1002/9780470973394

Goldstein, H., Carpenter, J., Kenward, M., & Levin, K. (2009). Multilevel models with multivariate mixed response types. Statistical Modelling, 9(3), 173-197. doi:10.1177/1471082X0800900301

See also

See also the MCMCglmm package for MCMC estimation and the lme4 package for maximum likelihood estimation.

Examples

if (FALSE) {
#############################################################################
# EXAMPLE 1: Multilevel model continuous data
#############################################################################

library(lme4)

#*** simulate data
set.seed(9097)

# number of clusters and units within clusters
K <- 150
n <- 15
iccx <- .2
idcluster <- rep( 1:K, each=n )
w <- stats::rnorm( K )
x <- rep( stats::rnorm( K, sd=sqrt(iccx) ), each=n) +
               stats::rnorm( n*K, sd=sqrt( 1 - iccx ))
X <- data.frame(int=1, "x"=x, xaggr=miceadds::gm(x, idcluster),
        w=rep( w, each=n ) )
X <- as.matrix(X)
Sigma <- diag( c(2, .5 ) )
u <- MASS::mvrnorm( K, mu=c(0,0), Sigma=Sigma )
beta <- c( 0, .3, .7, 1 )
Z <- X[, c("int", "x") ]
ypred <- as.matrix(X) %*% beta + rowSums( Z * u[ idcluster, ] )
y <- ypred[,1] + stats::rnorm( n*K, sd=1 )
data <- as.data.frame(X)
data$idcluster <- idcluster
data$y <- y

#*** estimate mixed effects model with miceadds::ml_mcmc() function
formula <- y ~ x + miceadds::gm(x, idcluster) + w + ( 1 + x | idcluster)
mod1 <- miceadds::ml_mcmc( formula=formula, data=data)
plot(mod1)
summary(mod1)

#*** compare results with lme4 package
mod2 <- lme4::lmer(formula=formula, data=data)
summary(mod2)

#############################################################################
# EXAMPLE 2: Multilevel model for ordinal outcome
#############################################################################

#*** simulate data
set.seed(456)
# number of clusters and units within cluster
K <- 500
n <- 10
iccx <- .2
idcluster <- rep( 1:K, each=n )
w <- rnorm( K )
x <- rep( stats::rnorm( K, sd=sqrt(iccx)), each=n) +
                 stats::rnorm( n*K, sd=sqrt( 1 - iccx ))
X <- data.frame("int"=1, "x"=x, "xaggr"=miceadds::gm(x, idcluster),
        w=rep( w, each=n ) )
X <- as.matrix(X)
u <- matrix( stats::rnorm(K, sd=sqrt(.5) ), ncol=1)
beta <- c( 0, .3, .7, 1 )
Z <- X[, c("int") ]
ypred <- as.matrix(X) %*% beta + Z * u[ idcluster, ]
y <- ypred[,1] + stats::rnorm( n*K, sd=1 )
data <- as.data.frame(X)
data$idcluster <- idcluster
alpha <- c(-Inf, -.4, 0, 1.7,  Inf)
data$y <- cut( y, breaks=alpha, labels=FALSE ) - 1

#*** estimate model
formula <- y ~ miceadds::cwc(x, idcluster) + miceadds::gm(x,idcluster) + w + ( 1 | idcluster)
mod <- miceadds::ml_mcmc( formula=formula, data=data, iter=2000, burnin=500,
                outcome="probit", inits_lme4=FALSE)
summary(mod)
plot(mod)
}