Profiling fails with prior weights

[stevencarlislewalker]

Here is an example,

library(lme4)
set.seed(1)
nGroups <- 100
nObs <- 1000
y <- rnorm(nObs)
g <- as.factor(rep(1:nGroups,each=nObs/nGroups))
w <- 1/(rpois(nObs,1) + 1)
mWeights <- lmer(y ~ 1 + (1|g), weights = w, REML = FALSE)
m <- lmer(y ~ 1 + (1|g), REML = FALSE)

The profile for the model with weights fails,

> profile(mWeights)
Error in profile.merMod(mWeights) : 
  Profiling over both the residual variance and
fixed effects is not numerically consistent with
profiling over the fixed effects only

while that for the model without weights works,

> profile(m)
         .zeta       .sig01    .sigma  (Intercept)        .par
1   0.00000000 0.000000e+00 1.0343983 -0.011648142      .sig01
2   0.01718913 1.000000e-03 1.0343981 -0.011648142      .sig01
3   0.18913968 1.100000e-02 1.0343768 -0.011648142      .sig01
4   0.61192660 3.548534e-02 1.0341815 -0.011648142      .sig01
5   1.03815225 5.986874e-02 1.0338163 -0.011648142      .sig01
...

I'm pretty sure this is an easy fix, arising because the log determinant of the prior weights matrix is not accounted for in devfun2. Here's the key line from devfun2,

pp$ldL2() + (resp$wrss() + pp$sqrL(1))/sigsq + n * log(2 * pi * sigsq)

[bbolker]

at least it fails rather than giving the wrong answer ...

[stevencarlislewalker]

Ya...I thought that was pretty cool too.

[stevencarlislewalker]

Almost positive my diagnosis is correct. Debugging through the example, the error is thrown because base (the devfun deviance) and dd(opt[seq(npar1)]) (the devfun2 deviance) differ by this non-zero amount,

Browse[2]> dd(opt[seq(npar1)]) - base
   .sigma 
-1215.922 

which is identical to the log determinant of the weights matrix,

Browse[2]> sum(log(environment(dd)$resp$weights))
[1] -1215.922

[stevencarlislewalker]

Not getting the error anymore after this commit but,

> confint(mWeights)
Computing profile confidence intervals ...
                  2.5 %    97.5 %
.sig01       0.00000000       Inf
.sigma       0.78951354 0.8618559
(Intercept) -0.09506946 0.0324549
Warning messages:
1: In nextpar(mat, cc, i, delta, lowcut, upcut) :
  Last two rows have identical or NA .zeta values: using minstep
2: In profile.merMod(object, signames = oldNames, ...) :
  non-monotonic profile
3: In optwrap(optimizer, par = start, fn = function(x) dd(mkpar(npar1,  :
  convergence code 3 from bobyqa: bobyqa -- a trust region step failed to reduce q
4: In optwrap(optimizer, par = start, fn = function(x) dd(mkpar(npar1,  :
  convergence code 3 from bobyqa: bobyqa -- a trust region step failed to reduce q
5: In optwrap(optimizer, par = start, fn = function(x) dd(mkpar(npar1,  :
  convergence code 3 from bobyqa: bobyqa -- a trust region step failed to reduce q
6: In optwrap(optimizer, par = thopt, fn = mkdevfun(rho, 0L), lower = fitted@lower) :
  convergence code 3 from bobyqa: bobyqa -- a trust region step failed to reduce q

This might be the right answer ... ?

[stevencarlislewalker]

I think these profiling warnings are not good. tl;dr: very wonky profiles for fixed effects when prior weights are present.

Here's a modification of the first example in tests/priorWeights.R,

library(lme4)
library(lattice)

set.seed(1)
nGroups <- 100
nObs <- 1000

# explanatory variable with a fixed effect
explVar1 <- rnorm(nObs)
explVar2 <- rnorm(nObs)

# random intercept among levels of a grouping factor
groupFac <- as.factor(rep(1:nGroups,each=nObs/nGroups))
randEff0 <- rep(rnorm(nGroups),each=nObs/nGroups)
randEff1 <- rep(rnorm(nGroups),each=nObs/nGroups)
randEff2 <- rep(rnorm(nGroups),each=nObs/nGroups)

# residuals with heterogeneous variance
residSD <- rpois(nObs,1) + 1
residError <- rnorm(nObs,sd=residSD)

# response variable
respVar <- randEff0 + (1+randEff1)*explVar1 + (1+randEff2)*explVar2 + residError

# rename to fit models on one line
y <- respVar
x <- explVar1
z <- explVar2
g <- groupFac
v <- residSD^2
w <- 1/v
w <- 0.01*w + 0.99  # very similar to unit weights

m <- lmer(y ~ x + (1|g), weights = w, REML = FALSE)
mProf <- profile(m)
densityplot(mProf)

These profiles can't be right given how similar this model is to one without any weights,

m <- lmer(y ~ x + (1|g), REML = FALSE)
mProf <- profile(m)
densityplot(mProf)

Could it be a V versus X thing?

[stevencarlislewalker]

I little more digging strongly suggests prior weights aren't being used in the profiled pred module, even though they are in the response module. Could be the source of the problem.

[stevencarlislewalker]

I now think this last comment is off base. I'm now intrigued by how setting delta = 1 in the call to profile fixes the wonky densityplots. Perhaps cutoff isn't being computed properly when prior weights are present? Nope...weights have absolutely nothing to do with cutoff.

[stevencarlislewalker]

Closing now because the original problem is fixed. There are still issues with the default options for profile.merMod (especially it seems with prior weights).

[stevencarlislewalker]

The problem persists. Here's a simple way to reproduce,

source("tests/priorWeights.R")
profile(lmerMods$ML1)