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

Dear Simon,
you could use the "estimate_density()" function from the bayestestR package
(see examples here:
https://easystats.github.io/see/articles/bayestestR.html), which works on
data frames (like the object returned by "profile()") and has an own
plot()-method. Another approach that *might* come close to what you would
like to achieve is simulating parameters from a multivariate normal
distribution, based on MASS:: mvrnorm(). This can be done with the
"simulate_parameters()" function from the parameters package, which also has
a print()-method (see visual examples here:
https://easystats.github.io/see/articles/parameters.html#simulated-model-par
ameters-1).

Here is some working code, I hope I got your point and provided a useful
answer...

library(lme4)
library(bayestestR)
library(parameters)
library(ggplot2)

theme_set(theme_classic())
hsb <-
read.csv('https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')
m31 <- lmer(math ~ ses*meanses + (ses | sch.id), data = hsb)
mp <- profile(m31)

# some profiled densities have very high spikes, so plot is probably less
useful...
d1 <- estimate_density(mp[c("ses", "meanses", "ses:meanses")])
d2 <- estimate_density(mp$ses)
d3 <- estimate_density(mp[c("ses", "ses:meanses")])

plot(d1) + ggplot2::xlim(c(0, 4))
plot(d2) + ggplot2::xlim(c(2.15, 2.25))
plot(d3) + ggplot2::xlim(c(0, 3))

sp <- simulate_parameters(m31)
plot(sp)
plot(sp, stack = FALSE)

Best
Daniel


-----Urspr?ngliche Nachricht-----
Von: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> Im
Auftrag von Simon Harmel
Gesendet: Montag, 26. Oktober 2020 01:55
An: r-sig-mixed-models <r-sig-mixed-models at r-project.org>
Betreff: [R-sig-ME] lme4: plotting profile density (not Zeta) manually not
by lattice

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?

library(lattice)
library(lme4)

hsb <- read.csv('
https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')
m31 <- lmer(math ~ ses*meanses + (ses | sch.id), data = hsb)

p <-  profile(m31)
densityplot(p)
Error in UseMethod("predict") :
  no applicable method for 'predict' applied to an object of class "NULL"


_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

--

_____________________________________________________________________

Universit?tsklinikum Hamburg-Eppendorf; K?rperschaft des ?ffentlichen Rechts; Gerichtsstand: Hamburg | www.uke.de
Vorstandsmitglieder: Prof. Dr. Burkhard G?ke (Vorsitzender), Joachim Pr?l?, Prof. Dr. Blanche Schwappach-Pignataro, Marya Verdel
_____________________________________________________________________

SAVE PAPER - THINK BEFORE PRINTING

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

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?

library(lattice)
library(lme4)

hsb <- read.csv('
https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')
m31 <- lmer(math ~ ses*meanses + (ses | sch.id), data = hsb)

p <-  profile(m31)
densityplot(p)
Error in UseMethod("predict") :
   no applicable method for 'predict' applied to an object of class "NULL"

	[[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

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)

    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?

library(lattice)
library(lme4)

hsb <- read.csv('
https://raw.githubusercontent.com/rnorouzian/e/master/hsb.csv')
m31 <- lmer(math ~ ses*meanses + (ses | sch.id), data = hsb)

p <-  profile(m31)
densityplot(p)
Error in UseMethod("predict") :
   no applicable method for 'predict' applied to an object of class

"NULL"

      [[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

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)

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

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

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

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

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

    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

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

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

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

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