[R-sig-ME] Nonlinear mixed models with nlme::nlme() and lme4::nlmer()

[David Sirl]

Hi all,

Following up my email below about getting different results when fitting 
the same model to the same data using nlme::nlme() and lme4::nlmer() 
[they agree fairly well on reported deviance, quartiles of residuals, 
estimated random effect variances and estimated fixed effects; but 
differ wildly on the standard errors and correlations associated with 
fixed effects]; I wonder if anyone has any suggestion/s for where or how 
I might be able to tempt someone to speculate what's going on please?

(I am very willing to follow instructions / try things out to help 
figure out what's happening; though due to other commitments I may be a 
few days between opportunities to do so.)

Best wishes,
Dave


On 21/11/2024 21:12, David Sirl wrote:
>
> Hi R-sig-mixed-models-ers,
>
> I'm fairly new to working with mixed models (and only have a couple of 
> years of working with R behind me) and have reason to want to fit some 
> nonlinear mixed models, so I'm working through the analysis of the 
> Soybean data that's in Sec 6.3 of Pinheiro & Bates' book 
> "Mixed-Effects Models in S and S-plus" (2000, Springer series on 
> Statistics & Computing).
>
> To help learn about these models I've been reproducing this analysis 
> using both nlme::nlme() (as P&B use, though in S rather than R) and 
> lme4::nlmer() (which I understand to have broadly similar 
> functionality but with better numerics in some places) and the results 
> from the two different packages agreed fairly well for the first 
> analyses I did. But now I have tried to fit a particular model (to the 
> same data) using these two methods and I get results that (i) suggest 
> I really am fitting the same model to the same data, (ii) in some 
> respects agree closely, but (iii) in other respects differ wildly.
>
> A couple of colleagues have checked my code and don't see any issues, 
> so I wonder if anyone here might be able to suggest what's going wrong 
> with either one of these functions (which I presume is not very 
> likely) or my use of these functions (presumably much more likely!) 
> please?
>
> I've set out what I've done below; and I'd be very grateful for any 
> suggestions that anyone can offer re what has gone wrong and where.
>
> Best wishes,
> Dave
>
>
> *A first sanity check*
>
> First I'll just write out the key parts of the calls I make to these 
> functions, in case I've been really daft and not noticed that my calls 
> should not be expected to fit the same model:
>
> nlme::nlme(weight ~ SSlogis(Time, Asym, xmid, scal),
>            data = ...,
>            fixed = Asym+xmid+scal~1,
>            random = Asym~1,
>            groups = ~Plot,
>            start = ...)
>
> lme4::nlmer(weight ~ SSlogis(Time, Asym, xmid, scal) ~ (Asym + xmid + 
> scal) + (Asym|Plot),
>             data = ...,
>             start = ...)
>
>
> *Some more detail*
>
> Assuming/hoping that the code above suggests that I am fitting the 
> same model with both methods, here's a bit more detail of the context, 
> what I'm trying to do and what's happening that suggests to me 
> something is going wrong.
>
> The dataset is Soybean (has "weight" recorded at 8-10 values of "Time" 
> for each of 48 values of the grouping factor "Plot") and the model I'm 
> trying to fit is of logistic growth using SSlogis(Time, Asym, xmid, 
> scal) with a random effect for Asym only (and fixed effects for all 
> three parameters).
>
> The starting parameter values I use are the (marginal) medians of the 
> parameter estimates for group-wise OLS fits of a logistic function to 
> the data for each Plot (fitted using nlme::olsList()) - I've omitted 
> this and entered the values explicitly here:
>
> startParams = c(Asym=19.813031, xmid=56.243394, scal=8.718537)
>
> fm3Soy.nlme =
> nlme::nlme(weight ~ SSlogis(Time, Asym, xmid, scal),
>              data = nlme::Soybean
>              fixed = Asym+xmid+scal~1,
>              random = Asym~1,
>              groups = ~Plot,
>              start = startParams)
>
>
> fm3Soy.nlmer =
> lme4::nlmer(weight ~ SSlogis(Time, Asym, xmid, scal) ~ (Asym + xmid + 
> scal) + (Asym|Plot),
>               data = nlme::Soybean
>               start = startParams)
>
> From the summary() of each fit I get what appear to be exactly the 
> same parameters being estimated / reported on, and they're very much 
> in agreement regarding likelihood/AIC/BIC values, variances of random 
> effects, point estimates of fixed effects. But they are qualitatively 
> very different in terms of standard errors and correlations for fixed 
> effects. The standard errors differ by a factor of ~100 and the 
> correlations are in the range 0.3 - 0.7 vs -0.3 - 0.0 and, based on 
> what I know of the data (and the output in P&B's book!)  I'm fairly 
> confident that the results from lme4::nlmer() which are incorrect.
>
> A further suggestion of things going wrong arises if I try to find CIs 
> for model parameters [using intervals(fm3Soy.nlme) and 
> confint(fm3Soy.nlmer)]. The latter gives an error message "Error in 
> eval(expr, envir, enclos): step factor reduced below 0.001 without 
> reducing pwrss" which I don't understand and for which my Google-ing 
> doesn't yield any help - I only see help for when this error arises in 
> fitting the model.
>
> Another possibly relevant point is that this is not a model that's 
> actually fitted in P&B (the first thing they do is deal with the 
> non-constant error variance). But even if some of the model 
> assumptions are violated I am still fitting the same model to the same 
> data, so I feel quite un-nerved by the big discrepancy in the results 
> from the two methods.
>
> The fits both seem reasonably good in the sense that the augmented 
> predictions [e.g. with nlme::augPred(fm3Soy.nlme, level=0:1) |> 
> plot()] for both methods agree with each other and, visually, agree 
> fairly well to the data. (It's only standard errors and correlations 
> of fixed effects that differ.)
>
> (I have also fitted the model in a Bayesian way using brms::brm() and 
> get results that are in keeping with those from nlme::nlme() above.)
>
> Thanks very much if you've read this far. I hope I've included all 
> relevant information to help track down what might be happening, but 
> if someone is interested in helping and would like some more 
> information then please do let me know.
>
	[[alternative HTML version deleted]]