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

[Ben Bolker]

    I've pushed an improved version of densityplot to github.  It 
creates density plots for all but .sig02 (the correlation parameter, I 
think), leaving that panel blank, and warns about skipped parameters. 
Give it a try ...

On 10/26/20 12:03 PM, Martin Maechler wrote:
>>>>>> 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
>