R devel - Sep 2013 - predictions in nlme without fixed covariantes

Dear all,

predict.lme() throws an error when the fixed part consists of only an intercept
and using newdata. See the reproducible example below. I've tracked the
error down to asOneFormula() which returns in this case NULL instead of a
formula. Changing NULL instead of ~1 in that function (see below) solves the
problem in the case of an intercept only model (m1). It does not solve the
problem in case of a model without intercept nor covariates (m2). The package
with altered asOneFormula() passes R CMD check on my machine.

Best regards,

Thierry Onkelinx

library(nlme)
data(Orthodont)
m0 <- lme(distance ~ Sex, random = ~1|Subject, data = Orthodont)
m1 <- lme(distance ~ 1, random = ~1|Subject, data = Orthodont)
m2 <- lme(distance ~ 0, random = ~1|Subject, data = Orthodont)

test.data <- Orthodont

test.data$Fitted <- predict(m0, level = 0)
test.data$Fitted.Newdata <- predict(m0, level = 0, newdata = test.data)
sum(abs(test.data$Fitted - test.data$Fitted.Newdata))

test.data$Fitted <- predict(m0, level = 1)
test.data$Fitted.Newdata <- predict(m0, level = 1, newdata = test.data)
sum(abs(test.data$Fitted - test.data$Fitted.Newdata))

test.data$Fitted <- predict(m1, level = 0)
test.data$Fitted.Newdata <- predict(m1, level = 0, newdata = test.data)
sum(abs(test.data$Fitted - test.data$Fitted.Newdata))

test.data$Fitted <- predict(m1, level = 1)
test.data$Fitted.Newdata <- predict(m1, level = 1, newdata = test.data)
sum(abs(test.data$Fitted - test.data$Fitted.Newdata))

test.data$Fitted <- predict(m2, level = 0)
test.data$Fitted.Newdata <- predict(m2, level = 0, newdata = test.data)
sum(abs(test.data$Fitted - test.data$Fitted.Newdata))

test.data$Fitted <- predict(m2, level = 1)
test.data$Fitted.Newdata <- predict(m2, level = 1, newdata = test.data)
sum(abs(test.data$Fitted - test.data$Fitted.Newdata))



#new version
asOneFormula <-
  ## Constructs a linear formula with all the variables used in a
  ## list of formulas, except for the names in omit
  function(..., omit = c(".", "pi"))
{
  names <- unique(allVarsRec((list(...))))
  names <- names[is.na(match(names, omit))]
  if (length(names)) {
    eval(parse(text = paste("~", paste(names, collapse =
"+")))[[1]])
  } else {
    ~ 1 #this was NULL
  }
}



sessionInfo()
R Under development (unstable) (2013-08-24 r63687)
Platform: i386-w64-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Dutch_Belgium.1252  LC_CTYPE=Dutch_Belgium.1252
[3] LC_MONETARY=Dutch_Belgium.1252 LC_NUMERIC=C
[5] LC_TIME=Dutch_Belgium.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

loaded via a namespace (and not attached):
[1] grid_3.1.0      lattice_0.20-15 tools_3.1.0
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en
binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer
and may not be regarded as stating an official position of INBO, as long as the
message is not confirmed by a duly signed document.

This should be fixed in nlme 3.1-112.

However, nlme has little support for formulae such as resp ~ 0, and does 
things like p:1 where p is the number of columns in the model matrix. 
3.1-112 does better but evidently the design did not consider this 
possibility.


On 30/09/2013 13:42, ONKELINX, Thierry wrote:
> Dear all,
>
> predict.lme() throws an error when the fixed part consists of only an
intercept and using newdata. See the reproducible example below. I've
tracked the error down to asOneFormula() which returns in this case NULL instead
of a formula. Changing NULL instead of ~1 in that function (see below) solves
the problem in the case of an intercept only model (m1). It does not solve the
problem in case of a model without intercept nor covariates (m2). The package
with altered asOneFormula() passes R CMD check on my machine.
>
> Best regards,
>
> Thierry Onkelinx
>
> library(nlme)
> data(Orthodont)
> m0 <- lme(distance ~ Sex, random = ~1|Subject, data = Orthodont)
> m1 <- lme(distance ~ 1, random = ~1|Subject, data = Orthodont)
> m2 <- lme(distance ~ 0, random = ~1|Subject, data = Orthodont)
>
> test.data <- Orthodont
>
> test.data$Fitted <- predict(m0, level = 0)
> test.data$Fitted.Newdata <- predict(m0, level = 0, newdata = test.data)
> sum(abs(test.data$Fitted - test.data$Fitted.Newdata))
>
> test.data$Fitted <- predict(m0, level = 1)
> test.data$Fitted.Newdata <- predict(m0, level = 1, newdata = test.data)
> sum(abs(test.data$Fitted - test.data$Fitted.Newdata))
>
> test.data$Fitted <- predict(m1, level = 0)
> test.data$Fitted.Newdata <- predict(m1, level = 0, newdata = test.data)
> sum(abs(test.data$Fitted - test.data$Fitted.Newdata))
>
> test.data$Fitted <- predict(m1, level = 1)
> test.data$Fitted.Newdata <- predict(m1, level = 1, newdata = test.data)
> sum(abs(test.data$Fitted - test.data$Fitted.Newdata))
>
> test.data$Fitted <- predict(m2, level = 0)
> test.data$Fitted.Newdata <- predict(m2, level = 0, newdata = test.data)
> sum(abs(test.data$Fitted - test.data$Fitted.Newdata))
>
> test.data$Fitted <- predict(m2, level = 1)
> test.data$Fitted.Newdata <- predict(m2, level = 1, newdata = test.data)
> sum(abs(test.data$Fitted - test.data$Fitted.Newdata))
>
>
>
> #new version
> asOneFormula <-
>    ## Constructs a linear formula with all the variables used in a
>    ## list of formulas, except for the names in omit
>    function(..., omit = c(".", "pi"))
> {
>    names <- unique(allVarsRec((list(...))))
>    names <- names[is.na(match(names, omit))]
>    if (length(names)) {
>      eval(parse(text = paste("~", paste(names, collapse =
"+")))[[1]])
>    } else {
>      ~ 1 #this was NULL
>    }
> }
>
>
>
> sessionInfo()
> R Under development (unstable) (2013-08-24 r63687)
> Platform: i386-w64-mingw32/i386 (32-bit)
>
> locale:
> [1] LC_COLLATE=Dutch_Belgium.1252  LC_CTYPE=Dutch_Belgium.1252
> [3] LC_MONETARY=Dutch_Belgium.1252 LC_NUMERIC=C
> [5] LC_TIME=Dutch_Belgium.1252
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> loaded via a namespace (and not attached):
> [1] grid_3.1.0      lattice_0.20-15 tools_3.1.0
> * * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
> Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver
weer en binden het INBO onder geen enkel beding, zolang dit bericht niet
bevestigd is door een geldig ondertekend document.
> The views expressed in this message and any annex are purely those of the
writer and may not be regarded as stating an official position of INBO, as long
as the message is not confirmed by a duly signed document.
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595