Bug 17226 - predict with gnls fails if not all levels are given
Summary: predict with gnls fails if not all levels are given
Status: UNCONFIRMED
Alias: None
Product: R
Classification: Unclassified
Component: Add-ons (show other bugs)
Version: R-devel (trunk)
Hardware: All All
: P5 normal
Assignee: R-core
URL:
Depends on:
Blocks: 17228
  Show dependency treegraph
 
Reported: 2017-02-21 22:52 UTC by Bill Denney
Modified: 2017-02-22 16:22 UTC (History)
0 users

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Bill Denney 2017-02-21 22:52:22 UTC
In the following example, predict only works when all levels of a predictor are given, but it should work when all levels or any subset is given by applying the original contrasts rather than re-calculating the contrasts for the new object.  The issue appears to be similar to one posted on Stack Exchange which has a suggested patch: http://stats.stackexchange.com/questions/29513/error-in-getting-predictions-from-a-lme-object

library(nlme)
d.mod <- mtcars
d.mod$gear <- factor(d.mod$gear)

mymod <-
  gnls(mpg~e.gear*disp,
       data=d.mod,
       params=e.gear~gear-1,
       start=rep(0.1, nlevels(d.mod$gear)))
summary(mymod)

# Fails "contrasts can be applied only to factors with 2 or more levels"
predict(mymod, newdata=data.frame(disp=100, gear=factor("3")))
# Fails "Error in p %*% beta[pmap[[nm]]] : non-conformable arguments"
predict(mymod, newdata=data.frame(disp=100, gear=factor(c("3", "4"))))
# Succeeds
predict(mymod, newdata=data.frame(disp=100, gear=factor(c("3", "4", "5"))))
Comment 1 Bill Denney 2017-02-22 04:43:15 UTC
Changing the function getParsGnls to the following will fix the second error:

getParsGnls <- function (plist, pmap, beta, N) {
  pars <- array(0, c(N, length(plist)), list(NULL, names(plist)))
  for (nm in names(plist)) {
    pars[, nm] <- if (is.logical(p <- plist[[nm]])) {
      beta[pmap[[nm]]]
      } else {
        betatmp <- beta[pmap[[nm]]]
        p.cols <- paste(nm, colnames(p), sep=".")
        if (!all(mask.p <- p.cols %in% names(betatmp))) {
          stop("Invalid level(s) for covariate ", nm, ": ", colnames(p)[!mask.p])
        }
        p %*% betatmp[colnames(p)]
      }
  }
  pars
}

The above is not the right solution for anything but contr.treatment, but hopefully it provides an idea of how the solution could be generated.
Comment 2 Bill Denney 2017-02-22 16:22:31 UTC
A better version of getParsGnls is in the related bug 17228.

Also, related fixes should also be applied to nlme (and probably other predict functions within the nlme library should also be examined:

library(nlme)
d.mod <- mtcars
d.mod$gear <- factor(d.mod$gear)
d.mod$fcyl <- factor(d.mod$cyl)

mymod.nlme <-
  nlme(mpg~e.gear*disp,
       data=d.mod,
       fixed=e.gear~gear-1,
       random=e.gear~1|fcyl,
       start=rep(0.1, nlevels(d.mod$gear)))
summary(mymod.nlme)

# Fails "contrasts can be applied only to factors with 2 or more levels"
predict(mymod.nlme, newdata=data.frame(disp=100, gear=factor("3"), fcyl=factor("4")))
# Fails "Error in p %*% beta[pmap[[nm]]] : non-conformable arguments"
predict(mymod.nlme, newdata=data.frame(disp=100, gear=factor(c("3", "4"), fcyl=factor("4"))))
predict(mymod.nlme, newdata=data.frame(disp=100, gear=factor(c("3", "5"), fcyl=factor("4"))))
# Succeeds
predict(mymod.nlme,
        newdata=data.frame(disp=100,
                           gear=factor(c("3", "4", "5")),
                           fcyl=factor(c("4", "6", "8"))))
# Fails but should provide NA
predict(mymod.nlme,
        newdata=data.frame(disp=100,
                           gear=factor(c("6", "7", "8")),
                           fcyl=factor(c("4", "6", "8"))))
predict(mymod.nlme,
        newdata=data.frame(disp=100,
                           gear=factor(c("3", "4", "5", "6", "7", "8")),
                           fcyl=factor(c("4", "6", "8", "4", "6", "8"))))
predict(mymod.nlme,
        newdata=data.frame(disp=100,
                           gear=factor(c("3", "4", "5", "3", "4", "5")),
                           fcyl=factor(c("4", "6", "8", "4", "6", "8"))))