lme4: plotting profile density (not Zeta) manually not by lattice

Simon Harmel 
    on Mon, 26 Oct 2020 09:51:26 -0500 writes:

> Ben,
    > I expect the exact same plots that densityplot(profile(fitted_model)) from
    > lattice produces?

    > again, densityplot(profile(fitted_model)) throws an error for the model in
    > my original question (and generally when any parameter's likelihood
    > distribution is highly spiked or funny-looking)

Hmm.. interesting.
As I'm coauthor of lme4  and have been doing nonparametric curve
estimation during my ph.d. years ("yesterday, .."),
I'm interested to rather fix the problem than try other
packages.

From your error message, there must be a buglet in either lattice

or lme4 ...

*BUT* (see below)

>> Can you clarify a bit what you want to plot?
    >> as.data.frame(p) is a good way to retrieve a simple data frame from
    >> profile objects that you can then transform/use to plot as you see fit.
    >> 
    >> Ben Bolker
    >>

>> > Dear All,
    >> >
    >> > I'm trying to plot the sampling distributions of my model parameters
    >> using `
    >> > densityplot()` from the `lattice` package but lattice often throws an
    >> error
    >> > even if one of the estimate's density distribution is highly skewed or
    >> > funny-looking.
    >> >
    >> > Is there a better package or a better way (even manually) to plot the
    >> > densities (not Zeta) from a `profile()` call?

   hsb <- read.csv('https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')

   library(lme4) # gets 'lattice' 
   m31 <- lmer(math ~ ses*meanses + (ses | sch.id), data = hsb)

   p <-  profile(m31)

## the profiling above gives *TONS and TONS* of warnings !

## so I guess now wonder you cannot easily plot it ..
## still you should at least get a better  error message

densityplot(p)

Error in UseMethod("predict") :
 no applicable method for 'predict' applied to an object of class
  NULL

Here's what I do to "summarize" ... and show "the solution" (?)

options(nwarnings=2^12) # so we store all the warnings !
system.time( p <- profile(m31) )
##   user  system elapsed 
## 19.007   0.002  19.111 
## There were 92 warnings (use warnings() to see them)

## MM: the cool thing is I wrote a summary() method for warnings in R 
##     a while ago, so use it:
summary( warnings() )
## Summary of (a total of 92) warning messages:
##  3x : In nextpar(mat, cc, i, delta, lowcut, upcut) :
##   unexpected decrease in profile: using minstep
## 88x : In nextpar(mat, cc, i, delta, lowcut, upcut) :
##   Last two rows have identical or NA .zeta values: using minstep
##  1x : In FUN(X[[i]], ...) : non-monotonic profile for .sig02

confint(p)

                 2.5 %     97.5 %
.sig01       1.4034755  1.8925324
.sig02      -0.9025412  0.2035804
.sig03       0.1824510  0.9800896
.sigma       5.9659398  6.1688744
(Intercept) 12.3231883 12.9337235
ses          1.9545565  2.4326048
meanses      3.0178000  4.5260869
ses:meanses -0.4044241  0.7279685
Warning messages:
1: In confint.thpr(p) :
  bad spline fit for .sig02: falling back to linear interpolation
2: In regularize.values(x, y, ties, missing(ties), na.rm = na.rm) :
  collapsing to unique 'x' values

so you see indeed, that  sig02 should probably be omitted from
the model

which I can "easily" confirm :

m30 <- lmer(math ~ ses * meanses + (1|sch.id) + (0+ ses | sch.id), data= hsb)
m20 <- lmer(math ~ ses * meanses + (1|sch.id), data= hsb)

anova(m31,m30,m20)

refitting model(s) with ML (instead of REML)
Data: hsb
Models:
m20: math ~ ses * meanses + (1 | sch.id)
m30: math ~ ses * meanses + (1 | sch.id) + (0 + ses | sch.id)
m31: math ~ ses * meanses + (ses | sch.id)
    npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)  
m20    6 46575 46616 -23282    46563                       
m30    7 46572 46620 -23279    46558 5.5415  1    0.01857 *
m31    8 46573 46628 -23278    46557 0.9669  1    0.32546  
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

So it seems m30, the model with no correlation between intercept
and slope fits well enough

and indeed,

system.time( p30 <- profile(m30 )
## ends in 5 sec, without any warnings,

and then

xyplot(p30)  # <-- more useful I think than
densityplot(p30) # both work fine

-- still I agree there's something we should do to fix the
   buglet !!

Martin Maechler
ETH Zurich