[R-sig-ME] nlme bug ?
Phillip Alday
ph||||p@@|d@y @end|ng |rom mp|@n|
Tue Aug 11 20:58:25 CEST 2020
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 using mpi.nl>
> *Date: *Monday, 10 August 2020 at 16:44
> *To: *"Cole, Tim" <tim.cole using ucl.ac.uk>,
> "r-sig-mixed-models using r-project.org" <r-sig-mixed-models using 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&data=02%7C01%7C%7Ca6498670eae541e3863a08d83d444fe4%7C1faf88fea9984c5b93c9210a11d9a5c2%7C0%7C1%7C637326710909416681&sdata=9iBQs56VTomb9EAJzND9ku7%2BSgZoohwwdJWxKQ%2Bi0Pc%3D&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 using r-project.org
> <mailto:R-sig-mixed-models using r-project.org> mailing list
>
> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=02%7C01%7C%7Ca6498670eae541e3863a08d83d444fe4%7C1faf88fea9984c5b93c9210a11d9a5c2%7C0%7C1%7C637326710909416681&sdata=te3VRPkQ3lk8Db6%2Frsdsl1py6dbbjWA2LIO2neTsWgc%3D&reserved=0
>
>
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list