[R-sig-ME] lme4: plotting profile density (not Zeta) manually not by lattice
Martin Maechler
m@ech|er @end|ng |rom @t@t@m@th@ethz@ch
Mon Oct 26 17:03:14 CET 2020
>>>>> 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
More information about the R-sig-mixed-models
mailing list