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.

> ## Example of scoping problem.
> ## Originally from a report by Markus Jantti:
> ## https://stat.ethz.ch/pipermail/r-help/2005-November/081382.html
> library(nlme)
> data(Ovary)
> ## stolen from example(anova.gls)
> # AR(1) errors within each Mare
> ## tolerance increased for flang (was 6e-6)
> fm1 <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), Ovary,
+            correlation = corAR1(form = ~ 1 | Mare))
> int1 <- intervals(fm1)
> ## no longer print attr(,"label"), PR#18196 :
> writeLines(intOut <- capture.output(int1))
Approximate 95% confidence intervals

 Coefficients:
                       lower       est.      upper
(Intercept)        10.908531 12.2163982 13.5242656
sin(2 * pi * Time) -4.044019 -2.7747122 -1.5054050
cos(2 * pi * Time) -2.272201 -0.8996047  0.4729919

 Correlation structure:
        lower      est.     upper
Phi 0.6684258 0.7532079 0.8186677

 Residual standard error:
   lower     est.    upper 
3.974697 4.616172 5.361173 
> stopifnot({
+     length(grep(',"label"', intOut, fixed=TRUE)) == 0
+     all.equal(int1$corStruct["Phi",],
+ 	      c(lower=0.66842829, est.=0.753207889, upper=0.81866619),
+ 	      tol = 1e-5)# 7e-6 needed for flan
+     all.equal(as.vector(int1$sigma),
+ 	      c(3.9747061, 4.61617157, 5.361161), tol = 1e-5)
+ })
> 
> # variance changes with a power of the absolute fitted values?
> fm2 <- update(fm1, weights = varPower())
> (a12 <- anova(fm1, fm2))
    Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fm1     1  5 1571.455 1590.056 -780.7273                        
fm2     2  6 1570.925 1593.247 -779.4626 1 vs 2 2.529261  0.1118
> stopifnot(identical(a12, anova(fm1, fm2, type = "seq")))# latter had failed
> 
> ## now define a little function
> dummy <- function(obj) anova(obj[[1]], obj[[2]])
> (d12 <- dummy(list(fm1, fm2)))
         Model df      AIC      BIC    logLik   Test  L.Ratio p-value
obj[[1]]     1  5 1571.455 1590.056 -780.7273                        
obj[[2]]     2  6 1570.925 1593.247 -779.4626 1 vs 2 2.529261  0.1118
> ## last failed < 3.1-66
> rownames(d12) <- rownames(a12)
> stopifnot(all.equal(a12, d12, tol = 1e-15),
+           all.equal(a12[2,"p-value"], 0.111752516, tol = 1e-5)
+ )
> 
> ## PR#13567
> fm1Orth.gls <- gls(distance ~ Sex * I(age - 11), Orthodont,
+                    correlation = corSymm(form = ~ 1 | Subject),
+                    weights = varIdent(form = ~ 1 | age))
> (aOr <- anova(fm1Orth.gls, Terms = "Sex"))
Denom. DF: 104 
 F-test for: Sex 
  numDF  F-value p-value
1     1 9.403075  0.0028
> stopifnot(all.equal(aOr[,"F-value"], 9.4030752449,
+ 		    aOr[,"p-value"], 0.0027608643857))
> 
> ## anova.gls(.) -- REML & ML
> (a1  <- anova(fm1))
Denom. DF: 305 
                   numDF  F-value p-value
(Intercept)            1 354.7332  <.0001
sin(2 * pi * Time)     1  18.5034  <.0001
cos(2 * pi * Time)     1   1.6633  0.1981
> (a1m <- anova(fm1, type="marginal"))
Denom. DF: 305 
                   numDF  F-value p-value
(Intercept)            1 337.8381  <.0001
sin(2 * pi * Time)     1  18.5034  <.0001
cos(2 * pi * Time)     1   1.6633  0.1981
> ##
> fm1M <- update(fm1, method = "ML")
> (a1M <- anova(fm1M))
Denom. DF: 305 
                   numDF  F-value p-value
(Intercept)            1 378.7745  <.0001
sin(2 * pi * Time)     1  19.1106  <.0001
cos(2 * pi * Time)     1   1.7133  0.1915
> (a1Mm <- anova(fm1M, type = "marginal"))
Denom. DF: 305 
                   numDF  F-value p-value
(Intercept)            1 360.2058  <.0001
sin(2 * pi * Time)     1  19.1106  <.0001
cos(2 * pi * Time)     1   1.7133  0.1915
> stopifnot(
+     all.equal(a1M[,"F-value"],
+               c(378.774471, 19.1105834, 1.71334909),
+               tolerance = 1e-7)
+    ,
+     all.equal(summary(fm1M)$tTable[,"t-value"] ^ 2,
+ 	      as.matrix(a1Mm)[,"F-value"], tolerance = 1e-14)
+    ,
+     all.equal(summary(fm1 )$tTable[,"t-value"] ^ 2,
+               as.matrix(a1m )[,"F-value"], tolerance = 1e-14)
+ )
> 
> proc.time()
   user  system elapsed 
   0.50    0.07    0.56