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

[Martin Maechler]

>>>>> Ben Bolker 
>>>>>     on Mon, 26 Oct 2020 22:38:10 -0400 writes:

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

That may be the best solution.

Note that xyplot()  which I'd prefer over densityplot() almost
always *does* deal with the situation,
using linear interpolation instead of the splines (which failed
here for .sig02).

We may also think of an alternative for parameters where the
splines had failed.  Simply skipping it for the plot seems "too
dangerous".

In your case the good model is really the m30 I used in my
previous e-mail, which indeed is the one where you set
the ".sig02"  \sigma_2 = 0

Martin


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