ml_mcmc: MCMC Estimation for Mixed Effects Model

In miceadds: Some Additional Multiple Imputation Functions, Especially for 'mice'

MCMC Estimation for Mixed Effects Model

Description

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

Usage

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 class 'ml_mcmc'
summary(object, digits=4, file=NULL, ...)

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

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

## S3 method for class '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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1177/1471082X0800900301")}

See Also

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

Examples

## Not run: 
#############################################################################
# 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)

## End(Not run)

miceadds documentation built on May 29, 2024, 11:05 a.m.