vdorie / blme

Bayesian Linear Mixed Effect Models
40 stars 8 forks source link

Sometimes can't run summary() on blme output #2

Closed longouyang closed 9 years ago

longouyang commented 9 years ago

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?

vdorie commented 9 years ago

Hi, thanks for this. I believe that I've fixed the issue in master, and should be pushing out an update to CRAN in a few days.

The convergence warnings are a product of lme4 changing its default optimizer. The old one is available by adding the argument control = lmerControl(optimizer = "Nelder_Mead").