Hello! For fixef.prior , it is in the format normal(sd = c(10, 2.5), cov, common.scale = TRUE) and it says Normal priors are constrained to have a mean of 0 - non-zero priors are equivalent to shifting covariates. Could you explain why non-zero priors are equivalent to shifting covariates?

Or how hard is it to allow non-zero fixed-effects priors? Do you have any clue how to incorporate the prior mean for fixed-effects? Thanks!

Hello! I am trying to fit a blme where the variance of the random effect is fixed and I found a nice solution for the residual variance for the linear model in blme (e.g. https://stat.ethz.ch/pipermail/r-sig-mixed-models/2014q1/021740.html). It seems that there isn't a way to specify point for the variance so I thought about writing a custom function (e.g. strongly penalizing any divergence from the fixed value). I encountered some issues, as documented below.

library(blme)
#> Loading required package: lme4
#> Loading required package: Matrix
y <- rbinom(100, 1, prob = 0.5)
x <- rnorm(100)
g <- sample(1:10, 100, replace = T)
fn_custom <- function(x){log(x)}

#Fails
bglmer(y ~ x + (1 | g), 
       cov.prior = custom(fn_custom),
       family = binomial())
#> Error in eval(parse(text = common.scale)[[1L]]): object 'none' not found
#Fails
bglmer(y ~ x + (1 | g), 
       cov.prior = custom(fn_custom, scale = "log"),
       family = binomial())
#> Error in validObject(.Object): invalid class "bmerCustomDist" object: invalid object for slot "commonScale" in class "bmerCustomDist": got class "function", should be or extend class "logical"
# Works!
blmer(y ~ x + (1 | g), 
       cov.prior = custom(fn_custom, scale = 'log'))
#> Cov prior  : g ~ custom(fn = fn_custom, chol = FALSE, scale = log, common.scale = TRUE)
#> Prior dev  : 4.0957
#> 
#> Linear mixed model fit by REML ['blmerMod']
#> Formula: y ~ x + (1 | g)
#> REML criterion at convergence: 153.9388
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  g        (Intercept) 0.1751  
#>  Residual             0.4875  
#> Number of obs: 100, groups:  g, 10
#> Fixed Effects:
#> (Intercept)            x  
#>    0.500282    -0.003755

^{Created on 2021-02-24 by the reprex package (v0.3.0)}

First off, thanks so much for writing this package - it's become indispensable in my work.

Both blmer() and bglmer() allow a list to be passed to cov.prior, i.e. bglmer(...,cov.prior=list(fc.nm ~ dist1, fc.nm ~ dist2)), but if we were to set the list as a variable:

cov.prior <-  list(fc.nm ~ dist1, fc.nm ~ dist2)
bglmer(...,cov.prior=cov.prior)

An error is thrown. Some cursory debugging indicates that when a variable that references a list is passed, it gets wrapped in a list twice, which causes problems in evaluateCovPriors(). I was able to hack together a solution but I wanted to raise the issue since there might need to be a more general solution.

Hello! I'm trying to play around with using level.dim to give programmatic arguments to a wishart prior (e.g. level.dim + 2). I find that for bglmer (not blmer), it fails to run. A reprex is below. I will fork the package and see if I can figure out the issue. It appears in both the CRAN and GitHub versions. Thanks so much!

library(lme4)
#> Loading required package: Matrix
library(blme)
library(lattice)
#level.dim works as a complex argument
gm1a <- blmer(incidence ~ period + (1 | herd),
              cbpp, cov.prior = wishart(df = level.dim))
#> Warning in log(diag(Lambda.t)): NaNs produced
#level dim does not work for binomial
gm1a <- bglmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
               cbpp, binomial, 
               cov.prior = wishart(df = level.dim))
#> Error in eval(matchedCall$df): object 'level.dim' not found
#level dim does not work for poisson
gm1a <- bglmer(incidence ~ period + (1 | herd),
               cbpp, poisson, 
               cov.prior = wishart(df = level.dim))
#> Error in eval(matchedCall$df): object 'level.dim' not found

^{Created on 2020-11-02 by the reprex package (v0.3.0)}

Hello! I'm a newer to GitHub.
I know this blme package from a paper：
Chung, Y. , Rabe-Hesketh, S. , Dorie, V. , & Andrew Gelman(2013). A nondegenerate penalized likelihood estimator for variance parameters in multilevel models. Psychometrika, 78(4), 685-709.
which suggests that it is helpful to choose a log-gamma penalty(or a weak information prior) for group-level variance.
But I wonder how can I get estimated standard error and confidence interval from this package? (for group-level variance)
Would you please help me? Thanks a lot.
Best wishes!

I've come across an example where blme seems to give weird results with nAGQ>1 (i.e., non-trivial Gauss-Hermite quadrature rather than Laplace approximation). I haven't thoroughly explored/gotten together a reproducible example yet, but thought I would post a comment here in the meantime (in case this has been reported before, or someone else finds it later before I have a chance to replicate/go into greater detail).

Hi there:

I'm trying to blme on a subset of the data from the radon example in the Gelman and Hill text on hierarchical models.

The data are the measurements from Minnesota

I tried running:

library(lme4)
library(blme)
d <- read.csv("Measurements-counties.csv")
blmer(activity ~ floor + Uppm + (1|county),
      data = d,
      cov.prior = gamma(shape = 1, rate = 1e-5),
      fixef.prior = normal(sd=sqrt(1e5), common.scale = FALSE),
      resid.prior = gamma(shape = 1, rate = 1e-5),
      REML = FALSE)

(Here, activity is the measured radiation levels, floor is which level of the house the measurements were taken on, county is the county that the house is in, and Uppm is the average soil uranium concentration in that county).

This runs okay, although it fails to converge (FWIW, I also tried reproducing the example from Chapter 10 of your dissertation and also got failure to converge). However, I can't run summary() on the blme result to examine standard errors of coefficient estimates. I get:

Error in diag(vcov(object, use.hessian = use.hessian)) : 
  error in evaluating the argument 'x' in selecting a method for function 'diag': Error in `dimnames<-`(`*tmp*`, value = list(c("(Intercept)", "floor",  : 
  invalid dimnames given for “dpoMatrix” object

(Some googling suggested that this problem might be due to a conflict with lme4. I tried it without loading lme4 and got an error about forceSymmetric being undefined)

summary() is not entirely broken, though, as it works if I don't specify priors:

summary(blmer(activity ~ floor + Uppm + (1|county),
      data = d,
      REML = FALSE))

(though this also warns about non-convergence)

Why can't I run summary() on the version with priors? Or, more germane to my current purposes, how would I go about extracting the coefficient standard errors for the priors version?

Is the resid.prior option intended to only be supported for blmer models? I was trying to set a gamma (1,1) prior on the residual variance (essentially, a unit exponential prior) which of course tends to be very useful, but bglmer reports that resid.prior is not an available option. It does appear to be a viable option for blmer however.

Hey, thanks for your work on this really useful software!

I was wondering if blme models are compatible with lme4::bootMer?

• lme4::bootMer(my_bmerMod, my_icc_function, nsim = 1000, use.u = FALSE, type = "parametric") gives no obvious errors.
• I get sensible looking results for my purpose (bootstrapping ICC).
• lme4::bootMer uses a generic refit method (which blme defines for bmerMods) so it isn't just running lmers.

Thanks in advance!

I want to fix the cov.prior to 10 and residual var to 5. But the cov.prior can not be passed to the model.
How to fix the cov.prior? Here I only have an intercept for one random effect 'rep'.
'''

Define the model formula

formula <- y ~ trt + (1|rep)

Define the control parameters

control = lmerControl(optimizer ='bobyqa')

Fit the model with a prior variance of 10 for the random intercepts

model <- blmer(formula, data=data, cov.prior = (diag(1)*10),
resid.prior = point(value=5, posterior.scale='var'),
control=control)
'''

Hi,

I would like to use a fixed effects prior in blme::bglmer() for one of my covariates showing complete separation, following the guidelines of a previous post: Binomial glmm with a categorical variable with full successes. However, in my case, there is an additional covariate, not showing complete separation.
Below is an overview of my dataset. ID identifies each observed fish, haul identifies each sampled haul (on a fishing vessel, with several fish from the same haul), death is my response variable with 2 levels (mortality at asymptote, 0 for survival and 1 for dead), R1 and R2 are my covariates each with 2 levels (reflex for each fish, 0 when reflex present and 1 when reflex absent):

ID   haul death R1 R2
2704    9     1  1  1
2705   14     1  1  0
2715   13     1  1  1
2718   13     1  1  0
2719    8     1  1  0
2722    9     1  1  1
...

Let's say that R1 is the covariate with complete separation for which I would like to add a fixed effects prior. My model would therefore be as follow:

bglmer(death ~ R1 + (1|haul) + (1|ID), data=RAMP_Subset, family=binomial, fixef.prior = normal(cov = diag(9,2)))

Now let's say that my second covariate R2 does not show complete separation. The reference manual for the R package "blme" states that priors on the fixed effects are specified in a fashion similar to that of covariance priors, so I tried the following:

bglmer(death ~ R1 + R2 + (1|haul) + (1|ID), data=RAMP_Subset, family=binomial, fixef.prior = R1 ~ normal(cov = diag(9,2)))

However, my attempt received an error message:

Error in evaluateFixefPrior(fixefPrior, defnEnv, evalEnv) : unrecognized fixef distribution: '~'

I could not find any other example online, and would be very grateful for advice on proper formatting.

Thanks,
Esther

As a possible feature, would it be possible to add a horseshoe prior option for the non-grouped variables? It is very useful in Stan but computationally time-consuming on large data sets. It may be much more effectively explored through an optimization approach. Ideally, this would be available through bayesglm as well to explore non-hierarchical versions of the same models.