-
Notifications
You must be signed in to change notification settings - Fork 0
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
[BUGZILLA #17712] nlme with varFixed multiplication error #6886
Comments
What is b in the formula?
[1] FALSE METADATA
|
Sorry I meant
[1] FALSE METADATA
|
Hi Benjamin,
Thanks, Bill library(nlme) nlme( nlme( mean(d$obs) METADATA
|
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
|
I’ve never had success with nlme with varFixed, so I don’t have a working example. METADATA
|
Looks like it was discussed even back in the S-PLUS days (2003) Douglas Bates replied to that with some information which I quote below, in case the wayback machine goes away.
METADATA
|
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) # nlme model fails # Equivalent lme model succeeds METADATA
|
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 ...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
|
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)) && METADATA
|
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
|
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
|
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
The text was updated successfully, but these errors were encountered: