[R-sig-ME] nlme bug ?
Cole, Tim
t|m@co|e @end|ng |rom uc|@@c@uk
Wed Sep 2 14:32:10 CEST 2020
A bit ago I posted some code suggesting an nlme bug. I’ve now realised it was my error, and for info I’m recording the details here.
The code, at the end of the thread, actually contains two errors. The first is hinted at in my exchange with Phillip, that there are three fixed effects for only two parameters. So the c fixed effect should be omitted.
The second bug is more subtle – the fitnlme function alters range(x), which is the default value for Boundary.knots in ns(.). So for fitnlme to work Boundary.knots needs to be constrained to be constant, by explicitly setting it to the initial range(x).
With these two changes the code works.
Best wishes,
Tim
``` r
suppressMessages(library(dplyr))
dat <- sitar::heights %>%
mutate(x = scale(age, scale = FALSE),
y = scale(height, scale = FALSE)) %>%
select(x, y, id)
bounds <- with(dat, range(x))
start <- rev(lm(y ~ splines::ns(x), data = dat)$coef)
fitnlme <- function(x, s1, a, c) {
a + drop(s1 * splines::ns(x * exp(c), B = bounds))
}
nlme::nlme(y ~ fitnlme(x, s1, a, c), fixed = s1 + a ~ 1, random = a + c ~ 1 | id, start = start, data = dat)
#> Nonlinear mixed-effects model fit by maximum likelihood
#> Model: y ~ fitnlme(x, s1, a, c)
#> Data: dat
#> Log-likelihood: -317.4776
#> Fixed: s1 + a ~ 1
#> s1 a
#> 48.53801 -18.49412
#>
#> Random effects:
#> Formula: list(a ~ 1, c ~ 1)
#> Level: id
#> Structure: General positive-definite, Log-Cholesky parametrization
#> StdDev Corr
#> a 5.5359464 a
#> c 0.2882972 -0.285
#> Residual 2.3225304
#>
#> Number of Observations: 124
#> Number of Groups: 12
```
<sup>Created on 2020-09-02 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: Tuesday, 11 August 2020 at 19:58
To: "Cole, Tim" <tim.cole using ucl.ac.uk>
Cc: "r-sig-mixed-models using r-project.org" <r-sig-mixed-models using r-project.org>
Subject: Re: [R-sig-ME] nlme bug ?
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<https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Freprex.tidyverse.org%2F&data=02%7C01%7C%7C4e3a3e85627446850a4308d83e288961%7C1faf88fea9984c5b93c9210a11d9a5c2%7C0%7C0%7C637327691137652975&sdata=RFEOtWcX5Zis44oS4zXQG%2Fgq9j9omUpwVZISuMaombI%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
From: Phillip Alday <phillip.alday using mpi.nl><mailto:phillip.alday using mpi.nl>
Date: Monday, 10 August 2020 at 16:44
To: "Cole, Tim" <tim.cole using ucl.ac.uk><mailto:tim.cole using ucl.ac.uk>, "r-sig-mixed-models using r-project.org"<mailto:r-sig-mixed-models using r-project.org> <r-sig-mixed-models using r-project.org><mailto: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]
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list