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

[Simon Harmel]

Awesome, thanks! That was what I was looking for.

On Mon, Oct 26, 2020 at 9:38 PM Ben Bolker <bbolker using gmail.com> wrote:

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

	[[alternative HTML version deleted]]