[R-sig-ME] lme4: plotting profile density (not Zeta) manually not by lattice

[Martin Maechler]

>>>>> 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)

    > On Mon, Oct 26, 2020, 8:19 AM Ben Bolker <bbolker using gmail.com> wrote:

    >> 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
    >> 
    >> On 10/25/20 8:54 PM, Simon Harmel wrote:
    >> > 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