Skip to content

Comparing variance components of crossed effects models fit with lme4 and nlme

6 messages · Thierry Onkelinx, Joshua Rosenberg, Ben Bolker

#
Hi all,

I'm trying to fit models with a) crossed random effects and b) a specific
residual structure (auto-correlation). Based on my understanding of what
nlme and lme4 do well, I would normally turn to lme4 to fit a model with
crossed random effects, but because I'm trying to structure the residuals,
I am trying nlme.

In trying to fit and compare the same variance components (no fixed
effects) model using lme4 and nlme, I found the output is similar but a bit
different. Specifically, the standard deviations of the random effects and
the log-likelihood statistics are different. Would you expect the output to
be a bit different?

The models I fit to compare the output are here, though the output is also
here:
https://bookdown.org/jmichaelrosenberg/comparing_crossed_effects_models/


library(lme4)
library(nlme)

m_lme4 <- lmer(diameter ~ 1 + (1 | plate) + (1 | sample), data = Penicillin)
m_lme4

m_nlme <- lme(diameter ~ 1, random = list(plate = ~ 1, sample = ~ 1), data
= Penicillin)
m_nlme


?Thank you for considering this question,
Josh?
#
Dear Joshua,

Crossed random effects are difficult to specify in nlme. I think that you
have to use pdBlocked() in the specification.

When I need correlation I often use INLA (r-inla.org). It allows for
correlated random effects. Crossed random effects are no problem.

Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

2017-08-10 23:05 GMT+02:00 Joshua Rosenberg <jmichaelrosenberg at gmail.com>:

  
  
3 days later
#
Thank you, I will explore use of INLA (or potentially the brms package
because of my familiarity with the [lme4-like] syntax).

I'm curious whether you (or anyone else) has thoughts / advice on using a
package that uses a Bayesian approach for carrying out mixed effects
modeling. In my field / area of research, mixed effects models are new! And
so a Bayesian approach to them would be *very *new. Even though if I
understand (very preliminarily), with some (uniform) prior specification,
results can be comparable to models specified with a maximum likelihood
approach, when possible.

Thank you again!
Josh

On Fri, Aug 11, 2017 at 8:38 AM, Thierry Onkelinx <thierry.onkelinx at inbo.be>
wrote:

  
    
#
A couple of thoughts:

(1) INLA *is* explicitly Bayesian, although I don't know what it
specifies (implicitly or explicitly) for priors or whether it allows
them to be user-adjusted (I'm too lazy to go look at the documentation
or Google "INLA priors" right now ...)
(2) it's worth making a distinction between
    (a) stochastic optimization (as in Bayesian MCMC, or frequentist
Monte Carlo expectation-maximization (MCEM) approaches) and
hill-climbing/deterministic optimization (as in INLA, or lme4, or
glmmTMB -- anything that says it uses the Laplace approximation, or
Gaussian quadrature ...)
   (b) inference based on a maximum point (MLE, or maximum a
posteriori [MAP] estimates in the Bayesian world) and inferences based
on means and distributions (MCMC). Typically the former goes with
deterministic optimization and the latter goes with stochastic
optimization
(3) in addition to INLA, there are a variety of existing Bayesian
machines in R (blme, MCMCglmm, brms, rstanarm ...) -- I think MCMCglmm
and brms implement some flavours of (temporal) autoregression ...

  Depending on the kind of autoregressive structure you want, glmmTMB
is also a possibility.

On Mon, Aug 14, 2017 at 12:23 PM, Joshua Rosenberg
<jmichaelrosenberg at gmail.com> wrote:
1 day later
#
Thank you for helping me build a "model" around these ideas and
distinctions.

Regarding point 2 (below), for software such as brms, I notice that some
kind of (sorry if my language is too loose here) estimate is produced for
fixed effects parameters, as well as their standard errors. Based on brm's
use (I think) of stochastic optimization and inferences based on means and
distributions, do these estimates correspond to something like the maximum
point - or rather to the mean of the posterior distribution?

thank you very much again
Josh
On Mon, Aug 14, 2017 at 6:54 PM, Ben Bolker <bbolker at gmail.com> wrote:

            

  
    
#
Most Bayesian approaches use the mean of the marginal posterior
distribution as the point estimate (although the median would be a
defensible choice as well). (I assume brms is using the mean but haven't
looked at the code/docs to confirm.)  The maximum of the marginal
posterior distribution would also be fairly easy to compute, but doesn't
have as much meaning in the Bayesian case.  The maximum (mode) of the
*multivariate* posterior distribution is the closest analogue of the
maximum likelihood estimate (with a completely uninformative prior it
would be identical), but is hard to compute.
On 17-08-16 01:21 PM, Joshua Rosenberg wrote: