nlme bug ?

Tue, Aug 11, 2020 11:58 AM

For the lme4 model fit with REML, you can get something comparable to
wait nlme calles the log-likelihood by dividing the REML criterion by
-2. (There have been a number of exchanges on this list, even in the
last year, about the problems of REML vs. _the_ likelihood, so I won't
go there, but just note that I can feel the Mathematical Powers That Be
looking disappointingly over my shoulder while I brush these problems
aside.)

On my machine that gives me:

nlme: -315.0308

lme: -313.6868

lmer: -313.6868

These are all quite similar and the corresponding coefficient values are
quite similar, so any differences are related to rounding issues and
machine-level variation in the nonlinear optimizer (which doesn't refer
to nlme here, but rather all of these tools actually fit the model).? (I
don't fully understand all the implementation details of the optimizers
themselves, but I have seen them take a different number of iterations
to fit the same model on different machines with ostensibly identical
software versions.) The slightly better fit returned by the linear
methods on my machine should generally be preferred -- but if your
inferences are sensitive to such small differences, you probably have
bigger problems. ;)

In other words, I wouldn't worry.


Phillip

On 11/8/20 3:35 pm, Cole, Tim wrote:

Thanks Phillip for your thoughts.

?

Thinking about it some more, I realise there is a redundancy which
explains why it fails. The construct

?

a + drop(s1 * splines::ns(x * exp(c)))

?

involves three parameters a, s1 and c, whereas it is a linear model
with intercept and slope and ought to have just two, i.e.

?

a + x * b.

?

This leads to another curiosity. This simpler model fits fine in nlme,
but unexpectedly it is not identical to the equivalent model fitted in
lme or lme4 ? see below. The mean relative difference between fitted
values is 1000 times greater for nlme vs lme than for lme vs lme4.

?

Is this simply a rounding issue, or are the nlme and lme models
genuinely different, and if so, in what way?

?

Thanks,

Tim

?

``` r

# install.packages('sitar')

suppressMessages({library(dplyr)

? library(sitar)

? library(lme4)

})

?

dat <- heights %>%

? mutate(x = scale(age, scale = FALSE)) %>%

? select(x, y = height, id)

?

start <- setNames(lm(y ~ x, data = dat)$coef, c('a', 'b'))

?

fitnlme <- function(x, a, b) {a + x * b}

?

nlme1 <- nlme(y ~ fitnlme(x, a, b), fixed = a + b ~ 1, random = a + b
~ 1 | id, start = start, data = dat)

?

lme1 <- lme(y ~ x, random = ~ x | id, data = dat)

lmer1 <- lmer(y ~ x + (x | id), data = dat)

?

all.equal(fitted(nlme1), fitted(lme1), check.attributes = FALSE)

#> [1] "Mean relative difference: 9.603422e-05"

all.equal(fitted(lme1), fitted(lmer1), check.attributes = FALSE)

#> [1] "Mean relative difference: 4.880377e-08"

```

?

<sup>Created on 2020-08-11 by the [reprex
package](https://reprex.tidyverse.org) (v0.3.0)</sup>

--

Population Policy and Practice

UCL Great Ormond Street Institute of Child Health,

30 Guilford Street, London WC1N 1EH, UK

?

?

*From: *Phillip Alday <phillip.alday at mpi.nl>
*Date: *Monday, 10 August 2020 at 16:44
*To: *"Cole, Tim" <tim.cole at ucl.ac.uk>,
"r-sig-mixed-models at r-project.org" <r-sig-mixed-models at r-project.org>
*Subject: *Re: [R-sig-ME] nlme bug ?

?

This could be a bug, but I would check two related things first:

?

1. Have you tried plotting a profile plot of the likelihood?

?

2. Have you tried? different starting values?

?

In my limited experience, these types of models can be surprisingly

sensitive to starting values.

?

Phillip

?

On 29/7/20 7:01 pm, Cole, Tim wrote:

    My sitar package fits the Lindstrom-Beath growth curve model,
    which varies the slope of the fitted spline curve y by scaling x.

    ?

    The model generally works well, but it fails in the simplest case,
    with a 1 df spline curve corresponding to a linear regression. The
    noddy example below fits a random intercept a and scaling factor
    c, with regression slope s1. The fitnlme function is called 20
    times and the values hardly change, and then it produces nonsense
    numbers and fails.

    ?

    Might this be a bug? It works fine with 2+ df for the spline curve.

    ?

    Thanks for your thoughts.

    ?

    Tim Cole

    ?

    ``` r

    # install.packages('sitar')

    # also uses splines and nlme

    suppressMessages(library(dplyr))

    ?

    dat <- sitar::heights %>%

    ?? select(x = 'age', y = 'height', id = 'id') %>%

    ?? mutate(x = scale(x, TRUE, FALSE))

    ?

    start <- lm(y ~ splines::ns(x), data = dat)$coef

    (start <- c(s1 = start[[2]], a = start[[1]], c = 0))

    #>????????s1???????? a???????? c

    #>??56.82702 128.76694?? 0.00000

    ?

    fitnlme <- function(x, s1, a, c) {

    ?? print(unique(s1))

    ?? print(unique(a))

    ?? print(unique(c))

    ?? a + drop(s1 * splines::ns(x * exp(c)))

    }

    ?

    nlme::nlme(y ~ fitnlme(x, s1, a, c), fixed = s1 + a + c ~ 1,
    random = a + c ~ 1 | id, start = start, data = dat)

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 1.490116e-08

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 1.490116e-08

    #> Warning in nlme.formula(y ~ fitnlme(x, s1, a, c), fixed = s1 +
    a + c ~ 1, :

    #> Iteration 1, LME step: nlminb() did not converge (code = 1).
    PORT message: false

    #> convergence (8)

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 1.490116e-08

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 0

    #> [1] 56.82702

    #> [1] 128.7669

    #> [1] 1.490116e-08

    #> [1] 48.22096

    #>??[1] 134.7411 130.4434 125.1878 134.2582 139.3737 132.7397
    120.4782 128.5690

    #>??[9] 137.1287 125.7633 136.5816 131.8558

    #>??[1]??168940.88??664404.25 1336620.21??242694.18
    -481809.06??256605.51

    #>??[7] 1934638.21??901713.94 -251016.54
    1287590.98??-99830.76??482104.58

    #> Error in qr.default(t(const)): NA/NaN/Inf in foreign function
    call (arg 1)

    ```

    ?

    <sup>Created on 2020-07-29 by the [reprex
    package](https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Freprex.tidyverse.org%2F&amp;data=02%7C01%7C%7Ca6498670eae541e3863a08d83d444fe4%7C1faf88fea9984c5b93c9210a11d9a5c2%7C0%7C1%7C637326710909416681&amp;sdata=9iBQs56VTomb9EAJzND9ku7%2BSgZoohwwdJWxKQ%2Bi0Pc%3D&amp;reserved=0)
    (v0.3.0)</sup>

    ?

    --

    Population Policy and Practice

    UCL Great Ormond Street Institute of Child Health,

    30 Guilford Street, London WC1N 1EH, UK

    ?

    ?

    ?????????? [[alternative HTML version deleted]]

    ?

    _______________________________________________

    R-sig-mixed-models at r-project.org
    <mailto:R-sig-mixed-models at r-project.org> mailing list

    https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&amp;data=02%7C01%7C%7Ca6498670eae541e3863a08d83d444fe4%7C1faf88fea9984c5b93c9210a11d9a5c2%7C0%7C1%7C637326710909416681&amp;sdata=te3VRPkQ3lk8Db6%2Frsdsl1py6dbbjWA2LIO2neTsWgc%3D&amp;reserved=0

?