[Rd] predictions in nlme without fixed covariantes
Prof Brian Ripley
ripley at stats.ox.ac.uk
Thu Oct 3 09:17:36 CEST 2013
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
More information about the R-devel
mailing list