blmer: Fit Bayesian Linear and Generalized Linear Mixed-Effects...

Description Usage Arguments Details Value See Also Examples

Description

Maximum a posteriori estimation for linear and generalized linear mixed-effects models in a Bayesian setting. Built off of lmer.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
blmer(formula, data, REML = TRUE,
      control = lmerControl(), start = NULL, verbose = 0L,
      subset, weights, na.action, offset, contrasts = NULL,
      devFunOnly = FALSE, cov.prior = wishart,
      fixef.prior = NULL, resid.prior = NULL, ...)
bglmer(formula, data, family = gaussian,
       control = glmerControl(), start = NULL, verbose = 0L,
       maxit = 100L,
       nAGQ = 1L, subset, weights, na.action, offset,
       contrasts = NULL, mustart, etastart,
       devFunOnly = FALSE, cov.prior = wishart,
       fixef.prior = NULL, ...)

Arguments

cov.prior

a BLME prior or list of priors with allowable distributions: wishart, invwishart, gamma, invgamma, or NULL. Imposes a prior over the covariance of the random effects/modeled coefficients. Default is wishart. The NULL argument imposes flat priors over all relevant parameters.

fixef.prior

a BLME prior of family normal, t, or NULL. Imposes a prior over the fixed effects/modeled coefficients. Default is NULL.

resid.prior

a BLME prior of family gamma, invamma, point or NULL. Imposes a prior over the noise/residual variance, also known as common scale parameter or the conditional variance given the random effects. Default is NULL.

start

like the start arguments for lmer and glmer a numeric vector or named list. Unlike the aforementioned, list members of fixef and sigma are applicable to linear mixed models provided that numeric optimization is required for these parameters.

formula, data, REML, family, control, verbose, maxit, nAGQ, mustart, etastart, devFunOnly, ...

model specification arguments as in lmer and glmer; see there for details.

subset, weights, na.action, offset, contrasts

further model specification arguments as in lm; see there for details.

Details

The bulk of the usage for blmer and bglmer closely follows the functions lmer and glmer. Those help pages provide a good overview of fitting linear and generalized linear mixed models. The primary distinction is that blmer and bglmer allow the user to do Bayesian inference or penalized maximum likelihood, with priors imposed on the different model components. For the specifics of any distribution listed below, see the distributions page.

Covariance Prior

The cov.prior argument applies a prior over the covariance matrix of the random effects/modeled coefficients. As there is one covariance matrix for every named grouping factor - that is every element that appears to the right of a vertical bar ("|") in the model formula - it is possible to apply as many different priors as there are said factors.

The general formats of an argument to blmer or bglmer for such a prior are of the form:

If the “factor.name ~” construct is ommitted, the prior is interpretted as a default and applied to all factors that lack specific priors of their own. Options are not required, but permit fine-tuning of the model.

Supported distributions are gamma, invgamma, wishart, invwishart, NULL, and custom.

The common.scale option, a logical, determines whether or not the prior applies to in the absolute-real world sense (value = FALSE), or if the prior is applied to the random effect covariance divided by the estimated residual variance (TRUE). As a practical matter, when false computation can be slower as the profiled common scale may no longer have a closed-form solution. As such, the default for all cases is TRUE.

Other options are specified along with the specific distributions and defaults are explained in the blme distributions page.

Fixed Effects Prior

Priors on the fixed effects, or unmodeled coefficients, are specified in a fashion similar to that of covariance priors. The general format is

At present, the implemented multivariate distributions are normal, t, and NULL. t priors cannot be used when REML is TRUE, as that integral does not have a closed form solution.

Residual Variance Prior

The general format for a residual variance prior is the same as for a fixed effect prior. The supported distributions are point, gamma, invgamma.

Value

An object of class "bmerMod", for which many methods are available. See there for details.

See Also

lmer, glmer, merMod class, and lm.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
##  covariance prior
data("sleepstudy", package = "lme4")
(fm1 <- blmer(Reaction ~ Days + (0 + Days|Subject), sleepstudy,
              cov.prior = gamma))
(fm2 <- blmer(Reaction ~ Days + (0 + Days|Subject), sleepstudy,
              cov.prior = gamma(shape = 2, rate = 0.5, posterior.scale = 'sd')))
(fm3 <- blmer(Reaction ~ Days + (1 + Days|Subject), sleepstudy,
              cov.prior = wishart))
(fm4 <- blmer(Reaction ~ Days + (1 + Days|Subject), sleepstudy,
              cov.prior = invwishart(df = 5, scale = diag(0.5, 2))))

##  custom prior
penaltyFn <- function(sigma)
  dcauchy(sigma, 0, 10, log = TRUE)
(fm5 <- blmer(Reaction ~ Days + (0 + Days|Subject), sleepstudy,
              cov.prior = custom(penaltyFn, chol = TRUE, scale = "log")))


##  fixed effect prior
(fm6 <- blmer(Reaction ~ Days + (1 + Days|Subject), sleepstudy,
              cov.prior = NULL,
              fixef.prior = normal))
(fm7 <- blmer(Reaction ~ Days + (1 + Days|Subject), sleepstudy,
              cov.prior = NULL,
              fixef.prior = normal(cov = diag(0.5, 2), common.scale = FALSE)))

##  residual variance prior
##  eight schools example
y <- c(28, 8, -3, 7, -1, 1, 18, 12);
sigma <- c(15, 10, 16, 11, 9, 11, 10, 18);
g <- 1:8;
  
(schools <- blmer(y ~ 1 + (1 | g), resid.prior = point,
                  cov.prior = NULL, REML = FALSE,
                  weights = 1 / sigma^2));

Example output

Loading required package: lme4
Loading required package: Matrix
Cov prior  : Subject ~ gamma(shape = 2.5, rate = 0, posterior.scale = sd, common.scale = TRUE)
Prior dev  : 4.2263

Linear mixed model fit by REML ['blmerMod']
Formula: Reaction ~ Days + (0 + Days | Subject)
   Data: sleepstudy
REML criterion at convergence: 1766.54
Random effects:
 Groups   Name Std.Dev.
 Subject  Days  7.108  
 Residual      29.080  
Number of obs: 180, groups:  Subject, 18
Fixed Effects:
(Intercept)         Days  
     251.41        10.47  
convergence code 0; 1 optimizer warnings; 0 lme4 warnings 
Warning message:
In get("checkConv", lme4Namespace)(attr(opt, "derivs"), opt$par,  :
  Model failed to converge with max|grad| = 0.550531 (tol = 0.002, component 1)
Cov prior  : Subject ~ gamma(shape = 2, rate = 0.5, posterior.scale = sd, common.scale = TRUE)
Prior dev  : 5.8346

Linear mixed model fit by REML ['blmerMod']
Formula: Reaction ~ Days + (0 + Days | Subject)
   Data: sleepstudy
REML criterion at convergence: 1766.54
Random effects:
 Groups   Name Std.Dev.
 Subject  Days  7.108  
 Residual      29.080  
Number of obs: 180, groups:  Subject, 18
Fixed Effects:
(Intercept)         Days  
     251.41        10.47  
convergence code 0; 1 optimizer warnings; 0 lme4 warnings 
Warning message:
In get("checkConv", lme4Namespace)(attr(opt, "derivs"), opt$par,  :
  Model failed to converge with max|grad| = 0.394449 (tol = 0.002, component 1)
Cov prior  : Subject ~ wishart(df = 4.5, scale = Inf, posterior.scale = cov, common.scale = TRUE)
Prior dev  : 2.4021

Linear mixed model fit by REML ['blmerMod']
Formula: Reaction ~ Days + (1 + Days | Subject)
   Data: sleepstudy
REML criterion at convergence: 1749.486
Random effects:
 Groups   Name        Std.Dev. Corr 
 Subject  (Intercept) 28.994        
          Days         9.864   -0.30
 Residual             24.631        
Number of obs: 180, groups:  Subject, 18
Fixed Effects:
(Intercept)         Days  
     251.41        10.47  
convergence code 0; 1 optimizer warnings; 0 lme4 warnings 
Warning message:
In get("checkConv", lme4Namespace)(attr(opt, "derivs"), opt$par,  :
  Model failed to converge with max|grad| = 3.39333 (tol = 0.002, component 1)
Cov prior  : Subject ~ invwishart(df = 5, scale = c(0.5, 0, 0, 0.5), posterior.scale = cov, common.scale = TRUE)
Prior dev  : -0.9156

Linear mixed model fit by REML ['blmerMod']
Formula: Reaction ~ Days + (1 + Days | Subject)
   Data: sleepstudy
REML criterion at convergence: 1752.517
Random effects:
 Groups   Name        Std.Dev. Corr 
 Subject  (Intercept) 17.57         
          Days        10.43    -0.19
 Residual             25.64         
Number of obs: 180, groups:  Subject, 18
Fixed Effects:
(Intercept)         Days  
     251.41        10.47  
convergence code 0; 1 optimizer warnings; 0 lme4 warnings 
Warning message:
In get("checkConv", lme4Namespace)(attr(opt, "derivs"), opt$par,  :
  Model failed to converge with max|grad| = 9.83116 (tol = 0.002, component 1)
Cov prior  : Subject ~ custom(fn = penaltyFn, chol = TRUE, scale = log, common.scale = TRUE)
Prior dev  : 6.8959

Linear mixed model fit by REML ['blmerMod']
Formula: Reaction ~ Days + (0 + Days | Subject)
   Data: sleepstudy
REML criterion at convergence: 1766.525
Random effects:
 Groups   Name Std.Dev.
 Subject  Days  7.26   
 Residual      29.02   
Number of obs: 180, groups:  Subject, 18
Fixed Effects:
(Intercept)         Days  
     251.41        10.47  
Fixef prior: normal(sd = c(10, 2.5), corr = 0, common.scale = TRUE)
Prior dev  : 24.0661

Linear mixed model fit by REML ['blmerMod']
Formula: Reaction ~ Days + (1 + Days | Subject)
   Data: sleepstudy
REML criterion at convergence: 1743.634
Random effects:
 Groups   Name        Std.Dev. Corr
 Subject  (Intercept) 24.772       
          Days         5.927   0.06
 Residual             25.503       
Number of obs: 180, groups:  Subject, 18
Fixed Effects:
(Intercept)         Days  
     251.23        10.47  
Fixef prior: normal(sd = c(0.7071, 0.7071), corr = 0, common.scale = FALSE)
Prior dev  : 2.2959

Linear mixed model fit by REML ['blmerMod']
Formula: Reaction ~ Days + (1 + Days | Subject)
   Data: sleepstudy
REML criterion at convergence: 1830.665
Random effects:
 Groups   Name        Std.Dev. Corr
 Subject  (Intercept) 252.50       
          Days         11.90   0.88
 Residual              25.59       
Number of obs: 180, groups:  Subject, 18
Fixed Effects:
(Intercept)         Days  
    0.03350      0.04581  
Resid prior: point(value = 1)
Prior dev  : 0

Linear mixed model fit by maximum likelihood  ['blmerMod']
Formula: y ~ 1 + (1 | g)
Weights: 1/sigma^2
     AIC      BIC   logLik deviance df.resid 
 65.3485  65.5868 -29.6742  59.3485        5 
Random effects:
 Groups   Name        Std.Dev.
 g        (Intercept) 0       
 Residual             1       
Number of obs: 8, groups:  g, 8
Fixed Effects:
(Intercept)  
      7.686  

blme documentation built on May 2, 2019, 10 a.m.