212 lines
7.1 KiB
Plaintext
212 lines
7.1 KiB
Plaintext
|
|
||
|
|
R version 3.5.2 Patched (2019-01-14 r75994) -- "Eggshell Igloo"
|
||
|
|
Copyright (C) 2019 The R Foundation for Statistical Computing
|
||
|
|
Platform: x86_64-pc-linux-gnu (64-bit)
|
||
|
|
|
||
|
|
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.
|
||
|
|
|
||
|
|
> library(nlme)
|
||
|
|
>
|
||
|
|
> options(digits = 6)# <==> less platform dependency in print() output
|
||
|
|
> if(!dev.interactive(orNone=TRUE)) pdf("test_lme.pdf")
|
||
|
|
>
|
||
|
|
> fm1 <- lmList(Oxboys)
|
||
|
|
> fm1
|
||
|
|
Call:
|
||
|
|
Model: height ~ age | Subject
|
||
|
|
Data: Oxboys
|
||
|
|
|
||
|
|
Coefficients:
|
||
|
|
(Intercept) age
|
||
|
|
10 130.262 3.72291
|
||
|
|
26 137.993 5.58878
|
||
|
|
25 139.210 4.02408
|
||
|
|
9 138.137 6.00906
|
||
|
|
2 142.858 5.44018
|
||
|
|
6 146.791 3.96317
|
||
|
|
7 146.128 4.99005
|
||
|
|
17 142.978 8.61178
|
||
|
|
16 147.545 4.54520
|
||
|
|
15 144.276 7.12426
|
||
|
|
8 148.293 6.46471
|
||
|
|
20 151.471 4.37447
|
||
|
|
1 148.120 7.17815
|
||
|
|
18 151.180 5.95779
|
||
|
|
5 151.429 6.24613
|
||
|
|
23 151.065 7.18512
|
||
|
|
11 150.047 8.50608
|
||
|
|
21 150.521 7.49779
|
||
|
|
3 155.651 4.77467
|
||
|
|
24 153.140 6.76470
|
||
|
|
22 154.567 8.08751
|
||
|
|
12 156.811 7.01547
|
||
|
|
13 156.071 8.49381
|
||
|
|
14 159.474 8.67089
|
||
|
|
19 164.576 9.06562
|
||
|
|
4 165.072 9.36056
|
||
|
|
|
||
|
|
Degrees of freedom: 234 total; 182 residual
|
||
|
|
Residual standard error: 0.659888
|
||
|
|
> fm2 <- lme(fm1)
|
||
|
|
> fm2
|
||
|
|
Linear mixed-effects model fit by REML
|
||
|
|
Data: Oxboys
|
||
|
|
Log-restricted-likelihood: -362.045
|
||
|
|
Fixed: height ~ age
|
||
|
|
(Intercept) age
|
||
|
|
149.37175 6.52547
|
||
|
|
|
||
|
|
Random effects:
|
||
|
|
Formula: ~age | Subject
|
||
|
|
Structure: General positive-definite, Log-Cholesky parametrization
|
||
|
|
StdDev Corr
|
||
|
|
(Intercept) 8.081077 (Intr)
|
||
|
|
age 1.680717 0.641
|
||
|
|
Residual 0.659889
|
||
|
|
|
||
|
|
Number of Observations: 234
|
||
|
|
Number of Groups: 26
|
||
|
|
> vc2 <- VarCorr(fm2)
|
||
|
|
> stopifnot(
|
||
|
|
+ all.equal(fixef(fm2), c("(Intercept)" = 149.371753,
|
||
|
|
+ age = 6.52546866), tol=1e-8),
|
||
|
|
+ all.equal(as.numeric(vc2[,"StdDev"]),
|
||
|
|
+ c(8.081077, 1.680717, 0.659889), tol=4e-7))
|
||
|
|
>
|
||
|
|
> # bug report from Arne.Mueller@sanofi-aventis.com
|
||
|
|
> mod <- distance ~ age + Sex
|
||
|
|
> fm3 <- lme(mod, Orthodont, random = ~ 1)
|
||
|
|
> pm3 <- predict(fm3, Orthodont)
|
||
|
|
> stopifnot(all.equal(mean(pm3), 24.023148148),
|
||
|
|
+ all.equal( sd(pm3), 2.4802195115),
|
||
|
|
+ all.equal(quantile(pm3), c("0%" = 17.0817792, "25%" = 22.3481813,
|
||
|
|
+ "50%" = 23.9271016, "75%" = 25.5740014,
|
||
|
|
+ "100%"= 30.8662241)))
|
||
|
|
>
|
||
|
|
>
|
||
|
|
> ## bug report and fix from Dimitris Rizopoulos and Spencer Graves:
|
||
|
|
> ## when 'returnObject = TRUE', do not stop() but give warning() on non-convergence:
|
||
|
|
> tools::assertWarning(
|
||
|
|
+ fm1 <- lme(distance ~ age, data = Orthodont,
|
||
|
|
+ control = lmeControl(msMaxIter = 1, returnObject = TRUE))
|
||
|
|
+ )
|
||
|
|
>
|
||
|
|
> ## "typo" in 'random=' -- giving 27-dim. vector random effect:
|
||
|
|
> ## PR#17524 B.Tyner: https://bugs.r-project.org/show_bug.cgi?id=17524
|
||
|
|
> try(lme(distance ~ 1, data=Orthodont, random = ~ Subject))
|
||
|
|
Error in lme.formula(distance ~ 1, data = Orthodont, random = ~Subject) :
|
||
|
|
fewer observations than random effects in all level 1 groups
|
||
|
|
> tools::assertError(lme(distance ~ age, data=Orthodont, random = ~ Subject))
|
||
|
|
> ## seg.faults in nlme <= 3.1-137 (2018) because of integer overflow
|
||
|
|
> ## The previous warning is now an *error* (unless new lmeControl(allow.n.lt.q=TRUE))
|
||
|
|
>
|
||
|
|
>
|
||
|
|
> ## based on bug report on R-help
|
||
|
|
> (p3.1 <- predict(fm3, Orthodont[1,]))
|
||
|
|
M01
|
||
|
|
25.3924
|
||
|
|
attr(,"label")
|
||
|
|
[1] "Predicted values (mm)"
|
||
|
|
> # failed in 3.1-88
|
||
|
|
> stopifnot(all.equal(pm3[1], p3.1,
|
||
|
|
+ check.attributes=FALSE, tolerance = 1e-14))
|
||
|
|
>
|
||
|
|
> ## Intervals failed in a patch proposal (Nov.2015):
|
||
|
|
> (fm4 <- lme(distance ~ age, Orthodont, random = ~ age | Subject))
|
||
|
|
Linear mixed-effects model fit by REML
|
||
|
|
Data: Orthodont
|
||
|
|
Log-restricted-likelihood: -221.318
|
||
|
|
Fixed: distance ~ age
|
||
|
|
(Intercept) age
|
||
|
|
16.761111 0.660185
|
||
|
|
|
||
|
|
Random effects:
|
||
|
|
Formula: ~age | Subject
|
||
|
|
Structure: General positive-definite, Log-Cholesky parametrization
|
||
|
|
StdDev Corr
|
||
|
|
(Intercept) 2.327034 (Intr)
|
||
|
|
age 0.226428 -0.609
|
||
|
|
Residual 1.310040
|
||
|
|
|
||
|
|
Number of Observations: 108
|
||
|
|
Number of Groups: 27
|
||
|
|
> i4 <- intervals(fm4)
|
||
|
|
> ## from dput(signif(i4$reStruct$Subject, 8))
|
||
|
|
> ## R-devel 2016-01-11; 64-bit :
|
||
|
|
> reSS <- data.frame(lower = c(0.9485605, 0.10250901, -0.93825047),
|
||
|
|
+ est. = c(2.3270341, 0.22642779, -0.60933286),
|
||
|
|
+ upper = c(5.7087424, 0.50014674, 0.29816857))
|
||
|
|
> ## R-devel 2016-01-11; 32-bit :
|
||
|
|
> ## reSS <- data.frame(lower = c(0.94962127,0.10262181, -0.93804767),
|
||
|
|
> ## est. = c(2.3270339, 0.22642779, -0.60933284),
|
||
|
|
> ## upper = c(5.7023648, 0.49959695, 0.29662651))
|
||
|
|
> rownames(reSS) <- rownames(i4$reStruct$Subject)
|
||
|
|
> sm4 <- summary(fm4)
|
||
|
|
> stopifnot(
|
||
|
|
+ all.equal(fixef(fm4),
|
||
|
|
+ c("(Intercept)" = 16.761111, age = 0.66018519)),
|
||
|
|
+ identical(fixef(fm4), sm4$tTable[,"Value"]),
|
||
|
|
+ all.equal(sm4$tTable[,"Std.Error"],
|
||
|
|
+ c("(Intercept)" = 0.77524603, age = 0.071253264), tol=6e-8),
|
||
|
|
+ all.equal(i4$reStruct$Subject[,"est."], reSS[,"est."], tol= 1e-7)
|
||
|
|
+ ## (lower, upper) cannot be very accurate for these : ==> tol = *e-4
|
||
|
|
+ ,## "interestingly" 32-bit values changed from 3.2.3 to R-devel(3.3.0):
|
||
|
|
+ all.equal(i4$reStruct$Subject[,c(1,3)], reSS[,c(1,3)], tol = .005)
|
||
|
|
+ ,
|
||
|
|
+ all.equal(as.vector(i4$sigma),
|
||
|
|
+ ## lower est. upper
|
||
|
|
+ c(1.0849772, 1.3100397, 1.5817881), tol=8e-4)
|
||
|
|
+ ,
|
||
|
|
+ all.equal(as.vector(i4$fixed),
|
||
|
|
+ as.vector(rbind(c(15.218322, 16.761111, 18.3039),
|
||
|
|
+ c(0.51838667, 0.66018519, 0.8019837))),
|
||
|
|
+ tol = 1e-6)
|
||
|
|
+ )
|
||
|
|
>
|
||
|
|
>
|
||
|
|
> ## wrong results from getData:
|
||
|
|
> ss2 <- readRDS("ss2.rds")
|
||
|
|
> m1 <- lme(PV1MATH ~ ESCS + Age +time ,
|
||
|
|
+ random = ~ time|SCHOOLID,
|
||
|
|
+ data = ss2,
|
||
|
|
+ weights = varIdent(form=~1|time),
|
||
|
|
+ corr = corCompSymm(form=~1|SCHOOLID/StIDStd),
|
||
|
|
+ na.action = na.omit)
|
||
|
|
> plot(m1, resid(.) ~ WEALTH)
|
||
|
|
>
|
||
|
|
> m2 <- lme(PV1MATH ~ ESCS + Age +time ,
|
||
|
|
+ random = ~ time|SCHOOLID,
|
||
|
|
+ data = ss2,
|
||
|
|
+ weights = varIdent(form=~1|time),
|
||
|
|
+ corr = corCompSymm(form=~1|SCHOOLID/StIDStd),
|
||
|
|
+ na.action = na.omit)
|
||
|
|
> plot(m2, resid(.) ~ WEALTH)
|
||
|
|
>
|
||
|
|
>
|
||
|
|
> ## Variogram() failing in the case of 1-observation groups (PR#17192):
|
||
|
|
> BW <- subset(BodyWeight, ! (Rat=="1" & Time > 1))
|
||
|
|
> if(interactive())
|
||
|
|
+ print( xtabs(~ Rat + Time, data = BW) )# Rat '1' only at Time == 1
|
||
|
|
> fm2 <- lme(fixed = weight ~ Time * Diet, random = ~ 1 | Rat, data = BW)
|
||
|
|
> Vfm2 <- Variogram(fm2, form = ~ Time | Rat)
|
||
|
|
> stopifnot(is.data.frame(Vfm2),
|
||
|
|
+ identical(dim(Vfm2), c(19L, 3L)),
|
||
|
|
+ all.equal(unlist(Vfm2[10,]), c(variog = 1.08575384191148,
|
||
|
|
+ dist = 22, n.pairs = 15))
|
||
|
|
+ )
|
||
|
|
> ## failed in nlme from 3.1-122 till 3.1-128
|
||
|
|
>
|
||
|
|
> proc.time()
|
||
|
|
user system elapsed
|
||
|
|
1.376 0.138 1.653
|