search

MichaelChirico / r-bugs

A ⚠️read-only⚠️mirror of https://bugs.r-project.org/
20 stars 0 forks source link

[BUGZILLA #17712] nlme with varFixed multiplication error #6886

Open MichaelChirico opened 4 years ago

MichaelChirico commented 4 years ago

This was initially reported here: https://stat.ethz.ch/pipermail/r-sig-mixed-models/2020q1/028351.html

When trying to fit an NLME model using a vector of varFixed, I get an error that appears to be related to a multiplication issue at this line of code:

https://github.com/cran/nlme/blob/c006dfa23ad390948a74e67978a9828a6d60d89b/R/varFunc.R#L169

Is the below a bug or am I inaccurately applying varFixed (or something else)?

Here is a reproducible example:

library(nlme) d <- data.frame( obs=rnorm(n=97), # 97 chosen because it's prime and therefore can't be the size of a rectangular matrix groups=rep(c("A", "B"), each=50)[1:97], wt=abs(rnorm(n=97)) )

nlme( ob∼b, fixed=∼1, random=∼1|groups, weights=varFixed∼wt), start=c(b=0), data=d )

Which gets the error:

Error in recalc.varFunc(object[[i]], conLin) : dims [product 12] do not match the length of object [97] In addition: Warning message: In conLin$Xy * varWeights(object) : longer object length is not a multiple of shorter object length

METADATA

MichaelChirico commented 4 years ago

What is b in the formula?

"d" %in% names(d)

[1] FALSE

METADATA

MichaelChirico commented 4 years ago

Sorry I meant

"b" %in% names(d)

[1] FALSE

METADATA

MichaelChirico commented 4 years ago

Hi Benjamin,

b isn't in the data, it's the estimated mean for obs in the data. nlme handles formulae differently than lme (and many other parts of R) where the formulae are considered as normal math as opposed to statistical notation. In the example below my signature, note how the example with the weights argument commented out works correctly and how the mean of d$obs is similar to the estimated value for b (not the same because of the random effect).

Thanks,

Bill

library(nlme) d <- data.frame( obs=rnorm(n=97), # 97 chosen because it's prime and therefore can't be the size of a rectangular matrix groups=rep(c("A", "B"), each=50)[1:97], wt=abs(rnorm(n=97)) )

nlme( ob∼b, fixed=∼1, random=∼1|groups, weights=varFixed∼wt), start=c(b=0), data=d )

nlme( ob∼b, fixed=∼1, random=∼1|groups,

weights=varFixed∼wt),

start=c(b=0), data=d )

mean(d$obs)

METADATA

MichaelChirico commented 4 years ago

Thank you Bill, I hadn't seriously used nlme since grad school and forgotten that...

It appears you are not the first to experience this issue, for example see this report from 2012,

https://stat.ethz.ch/pipermail/r-help/2012-March/307356.html

Do you happen to have a comparable example that works as intended? I.e. one that follows the same codepath, still using weights=varFixed∼wt), and preferably still using 97 observations.

METADATA

MichaelChirico commented 4 years ago

I’ve never had success with nlme with varFixed, so I don’t have a working example.

METADATA

MichaelChirico commented 4 years ago

Looks like it was discussed even back in the S-PLUS days (2003)

https://web.archive.org/web/20041225075947/http://www.biostat.wustl.edu/archives/html/s-news/2003-09/msg00013.html

Douglas Bates replied to that with some information which I quote below, in case the wayback machine goes away.

conLin is an internal representation of a linear mixed-effects model
used in the lme step.  (The nlme function alternates between lme steps
and penalized nonlinear least squares steps when performing parameter
estimation.) The columns of conLin$Xy are the derivatives of the model
function with respect to the random effects and the fixed effects
followed by the working residuals.  In your case the fourth through
the sixth columns should be repetitions of the first through the
third.

There is a point in the lme optimization where the information in Xy
is condensed by means of a series of QR decompositions.  There will be
3 rows (the dimension of the random effects vector) retained for each
group, 3 rows for the fixed effects and 1 row for the response.
That's where the 46 rows come from.

METADATA

MichaelChirico commented 4 years ago

Thanks for that research! When the fix is put in place, perhaps Doug Bates' comments can be added as a code comment to help with future modifications.

Hopefully with the reports over all this time, we can fix it this time!

That said, if we are looking at a QR decomposition of the LME part, I'm surprised that it doesn't work because the equivalent LME model works.

library(nlme) d <- data.frame( obs=rnorm(n=97), # 97 chosen because it's prime and therefore can't be the size of a rectangular matrix groups=rep(c("A", "B"), each=50)[1:97], wt=abs(rnorm(n=97)) )

# nlme model fails nlme( ob∼b, fixed=∼1, random=∼1|groups, weights=varFixed∼wt), start=c(b=0), data=d )

# Equivalent lme model succeeds lme( ob∼1, random∼1|groups, weights=varFixed∼wt), data=d )

METADATA

MichaelChirico commented 4 years ago

In the code, decomp=TRUE signals that QR-decomposition is used, with the non-QR-decomposed version of conLin stored in oldConLin which gets used in various places throughout the code, in particular for this block:

## wrapping up nlmeFit <- if (decomp) MEestimate(nlmeSt, grpShrunk, oldConLin) else MEestimate(nlmeSt, grpShrunk)

...from which one might infer the author never intended for MEestimate to support a QR-decomposed conLin? If so, one might also suppose that when decomp is TRUE, oldConLin should need supplying here:

nlmeFit <- attr(nlmeSt, "lmeFit") <- MEestimate(nlmeSt, grpShrunk)

...but it seems doing so in isolation still does not suffice to remedy the issue, as conLin$Xy also gets used in lmeApVar via the full log-likelihood function (lmeApVar.fullLmeLogLik) passed to fdHess.

METADATA

MichaelChirico commented 4 years ago

I'm not positive that I'm following your logic correctly, but might it be as simple as setting decomp=FALSE in this situation? (One difficulty here is I'm not sure how to correctly define "this situation".)

I think that the change would occur at this location:

if (length(coef(nlmeSt)) == length(coef(nlmeSt$reStruct)) && !needUpdate(nlmeSt)) { # can do one decomposition ## need to save conLin for calculating updating between steps oldConLin <- attr(nlmeSt, "conLin") decomp <- TRUE } else decomp <- FALSE

METADATA

MichaelChirico commented 4 years ago

Right...hoping to avoid throwing out the decomp baby with the bathwater...hence my question about whether you had a comparable varFixed example that works as intended.

One idea is to also check length(nlmeSt): in your example it has length 2; maybe the decomp methodology can only work length(nlmeSt) == 1? That is the way lme determines decomp, which I think explains why your 'equivalent lme model succeeds'.

METADATA

MichaelChirico commented 4 years ago

We are rapidly moving past my knowledge of internals of (non)linear mixed modeling.

If testing length(nlmeSt) is the way that lme works, it seems rational that the same limit would apply to nlme in this case, but that is not based on specific knowledge of how it is supposed to work just the general rationale that something like that should work.

METADATA