lme4: plotting profile density (not Zeta) manually not by lattice

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