R Under development (unstable) (2024-06-05 r86695 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ## from PR#8905
> library(nlme)
> data(Orthodont)
> fm <- lme(distance ~ poly(age, 3) + Sex, data = Orthodont, random = ~ 1)
> # data for predictions
> Newdata <- head(Orthodont)
> Newdata$Sex <- factor(Newdata$Sex, levels = levels(Orthodont$Sex))
> (pr <- predict(fm, Newdata))
     M01      M01      M01      M01      M02      M02 
25.52963 26.51111 27.99259 29.43703 21.75617 22.73765 
attr(,"label")
[1] "Predicted values (mm)"
> stopifnot(all.equal(c(pr), fitted(fm)[1:6]))
> 
> ## https://stat.ethz.ch/pipermail/r-devel/2013-September/067600.html
> ## but with a different fix.
> 
> m0 <- lme(distance ~ Sex, random = ~1|Subject, data = Orthodont)
> Fitted <- predict(m0, level = 0)
> Fitted.Newdata <- predict(m0, level = 0, newdata = Orthodont)
> stopifnot(sum(abs(Fitted - Fitted.Newdata)) == 0)
> 
> Fitted <- predict(m0, level = 1)
> Fitted.Newdata <- predict(m0, level = 1, newdata = Orthodont)
> sum(abs(Fitted - Fitted.Newdata))
[1] 0
> stopifnot(sum(abs(Fitted - Fitted.Newdata)) == 0)
> 
> m1 <- lme(distance ~ 1, random = ~1|Subject, data = Orthodont)
> Fitted <- predict(m1, level = 0)
> Fitted.Newdata <- predict(m1, level = 0, newdata = Orthodont)
> stopifnot(sum(abs(Fitted - Fitted.Newdata)) == 0)
> 
> Fitted <- predict(m1, level = 1)
> Fitted.Newdata <- predict(m1, level = 1, newdata = Orthodont)
> stopifnot(sum(abs(Fitted - Fitted.Newdata)) == 0)
> 
> m2 <- lme(distance ~ 0, random = ~1|Subject, data = Orthodont)
> Fitted <- predict(m2, level = 0)
> Fitted.Newdata <- predict(m2, level = 0, newdata = Orthodont)
> stopifnot(sum(abs(Fitted - Fitted.Newdata)) == 0)
> 
> Fitted <- predict(m2, level = 1)
> Fitted.Newdata <- predict(m2, level = 1, newdata = Orthodont)
> stopifnot(sum(abs(Fitted - Fitted.Newdata)) == 0)
> 
> 
> m3 <- lme(fixed = distance ~ age, data = Orthodont,
+           random = ~ 1 | Subject)
> m4 <- update(m3, random = ~ age | Subject)
> m5 <- update(m4, fixed. = distance ~ age * Sex)
> 
> newD <- expand.grid(age = seq(7,15, by = .25),
+                     Sex = c("Male", "Female"),
+                     Subject = c("M01", "F01"))
> (n.age <- attr(newD, "out.attrs")$dim[["age"]]) # 33
[1] 33
> str(p5 <- predict(m5, newdata = newD, asList = TRUE, level=0:1))
'data.frame':	132 obs. of  3 variables:
 $ Subject        : Factor w/ 2 levels "M01","F01": 1 1 1 1 1 1 1 1 1 1 ...
 $ predict.fixed  : num  21.8 22 22.2 22.4 22.6 ...
 $ predict.Subject: num  24 24.2 24.4 24.6 24.8 ...
> pp5 <- cbind(newD, p5[,-1])
> stopifnot(identical(colnames(pp5),
+                     c("age", "Sex", "Subject", "predict.fixed", "predict.Subject")))
> fixef(m5) # (Intercept) age SexF age:SexF
  (Intercept)           age     SexFemale age:SexFemale 
   16.3406250     0.7843750     1.0321023    -0.3048295 
> p5Mf <- pp5[pp5$Sex == "Male", "predict.fixed"]
> p5MS <- subset(pp5, subset = Subject == "M01" & Sex == "Male",
+                select = "predict.Subject", drop=TRUE)
> X.1 <- cbind(1, newD[1:n.age,"age"])
> stopifnot(all.equal(p5Mf[  1:n.age],
+                     p5Mf[-(1:n.age)], tol = 1e-15)
+          ,
+           all.equal(p5Mf[  1:n.age],
+                     c(X.1 %*% fixef(m5)[1:2]), tol = 1e-15)
+          ,
+           all.equal(p5MS,
+                     c(X.1 %*% (fixef(m5)[1:2] + as.numeric(ranef(m5)["M01",]))), tole = 1e-15)
+           )
> 
> 
> ## PR#18312: predict with character vs. factor variables in newdata
> newOrth <- data.frame(Subject = "F03", Sex = "Female", age = 8,
+                       stringsAsFactors = FALSE)  # default in R >= 4.0.0
> stopifnot(all.equal(predict(m5, newdata = newOrth), # this failed
+                     fitted(m5)["F03"], # first obs of F03 is for age 8
+                     check.attributes = FALSE))
> ## failed in nlme <= 3.1-155 with
> ## Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
> ##   contrasts can be applied only to factors with 2 or more levels
> newOrth2 <- rbind(newOrth, list("M11", "Male", 16))
> stopifnot(all.equal(
+     predict(m5, newdata = type.convert(newOrth2, as.is = FALSE)), # factor
+     predict(m5, newdata = newOrth2)                               # character
+ ))
> ## predictions with character input were *wrong* in nlme <= 3.1-155
> 
> ## numeric newdata for a factor variable should at least warn
> tools::assertWarning(predict(m5, newdata = transform(newOrth, Sex = 2)), verbose = TRUE)
Asserted warning: variable 'Sex' is not a factor
> ## did not warn in nlme <= 3.1-155 and may return unexpected result
> ## (not the same as Sex=factor("Female", levels = c("Male", "Female")))
> 
> ## intercept-free model
> m0b <- lme(distance ~ Sex - 1, random = ~1|Subject, data = Orthodont)
> stopifnot(all.equal(predict(m0b, Orthodont[1,], level=0),
+                     fixef(m0b)[1], check.attributes = FALSE))
> ## predict wrongly returned c(0,0) in nlme <= 3.1-155
> 
> 
> ##--- simulate():---------
> 
> ## border cases
> ort.0 <- simulate(m3, method = character())# "nothing" stored
> ort.M <- simulate(m3, method = "ML",   seed=47)
> ort.R <- simulate(m3, method = "REML", seed=47)
> stopifnot(identical(names(ort.0), "null"),
+           identical(names(ort.M), "null"),
+           identical(names(ort.R), "null"),
+           identical(ort.0$null, list()),
+           identical(names(ort.M$null), "ML"),
+           identical(names(ort.R$null), "REML"),
+           all.equal(loM <- ort.M$null$ML  [,"logLik"], -215.437, tol = 2e-6)
+          ,
+           all.equal(loR <- ort.R$null$REML[,"logLik"], -217.325, tol = 2e-6)
+           )
> 
> system.time(
+  orthS3 <- simulate.lme(list(fixed = distance ~ age, data = Orthodont,
+                               random = ~ 1 | Subject), nsim = 3,
+                          m2 = list(random = ~ age | Subject), seed = 47)
+ )
   user  system elapsed 
   0.11    0.00    0.11 
> ## the same, starting from two fitted models :
> ort.S3 <- simulate(m3, m2 = m4, nsim = 3, seed = 47)
> attr(ort.S3, "call") <- attr(orthS3, "call")
> ## was 1e-15, larger tolerance needed with ATLAS
> stopifnot(all.equal(orthS3, ort.S3, tolerance = 1e-10))
> 
> logL <- sapply(orthS3, function(E) sapply(E,
+                        function(M) M[,"logLik"]), simplify="array")
> 
> stopifnot(is.array(logL), length(d <- dim(logL)) == 3, d == c(3,2,2),
+     sapply(orthS3, function(E) sapply(E, function(M) M[,"info"])) == 0
+    , # typically even identical(), but not with ATLAS
+     all.equal(logL[1,,"null"], c(ML = loM, REML = loR), tol = 1e-10)
+    ,
+     all.equal(c(logL) + 230,
+               c(14.563 , 2.86712, 1.00026, 12.6749, 1.1615,-0.602989,
+                 16.2301, 2.95877, 2.12854, 14.3586, 1.2534, 0.582263), tol=8e-6)
+ )
> 
> proc.time()
   user  system elapsed 
   0.53    0.06    0.59