miceadds: ml_mcmc – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

ml_mcmc

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_iter`th 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=\boldmath{X}\boldmath{beta}+ \boldmath{Z}\boldmath{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

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

Some Additional Multiple Imputation Functions, Especially for 'mice'

v3.11-6

GPL (>= 2)

Authors

Alexander Robitzsch [aut,cre] (<https://orcid.org/0000-0002-8226-3132>), Simon Grund [aut] (<https://orcid.org/0000-0002-1290-8986>), Thorsten Henke [ctb]

Initial release

2021-01-21 11:48:47