#  File src/library/stats/tests/nafns.R
#  Part of the R package, https://www.R-project.org
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  https://www.R-project.org/Licenses/

## Tests of functions handling NAs in fits
## These functions were introduced in 1.3.0.
## They are used by lm and glm in base R, and by
## packages MASS, rpart and survival.

## Comparison strictness is set at 256eps, increased from 100eps for v.2.1.0
## May look lenient, but notice that the Ozone levels that are being modeled
## can be more than 100.

dim(airquality)
nd <- airquality[c(6,25:27), ]

sm <- function(x) cat("length", length(x), "with", sum(is.na(x)), "NAs\n")

# default is to omit some rows
fit <- lm(Ozone ~ ., data=airquality, na.action=na.omit)
summary(fit)
sm(fitted(fit))
sm(resid(fit))
sm(predict(fit))
(pp <- predict(fit, nd))

fit2 <- lm(Ozone ~ ., data=airquality, na.action=na.exclude)
summary(fit2) # same as before
sm(fitted(fit2))
sm(resid(fit2))
sm(predict(fit2))
(pp2 <- predict(fit2, nd))

## same as before: napredict is only applied to predictions on the
## original data, following Therneau's original code (and S-PLUS).
## However, as from R 1.8.0 there is a separate na.action arg to predict.lm()
stopifnot(all.equal(pp, pp2))

## should fail
try(fit3 <- lm(Ozone ~ ., data=airquality, na.action=na.fail))

## more precise tests.
f1 <- fitted(fit)
f2 <- fitted(fit2)
common <- match(names(f1), names(f2))
stopifnot(max(abs(f1 - f2[common])) <= 256*.Machine$double.eps)
stopifnot(all(is.na(f2[-common])))

r1 <- resid(fit)
r2 <- resid(fit2)
common <- match(names(r1), names(r2))
stopifnot(max(abs(r1 - r2[common])) <= 256*.Machine$double.eps)
stopifnot(all(is.na(r2[-common])))

p1 <- predict(fit)
p2 <- predict(fit2)
common <- match(names(p1), names(p2))
stopifnot(max(abs(p1 - p2[common])) <= 256*.Machine$double.eps)
stopifnot(all(is.na(p2[-common])))


### now try out glm
gfit <- glm(Ozone ~ ., data=airquality, na.action=na.omit)
summary(gfit)
sm(fitted(gfit))
sm(resid(gfit))
sm(predict(gfit))
predict(gfit, nd)
(pp <- predict(gfit, nd))

gfit2 <- glm(Ozone ~ ., data=airquality, na.action=na.exclude)
summary(gfit2) # same as before
sm(fitted(gfit2))
sm(resid(gfit2))
sm(predict(gfit2))
(pp2 <- predict(gfit2, nd))
stopifnot(all.equal(pp, pp2))

## more precise tests.
f1 <- fitted(gfit)
f2 <- fitted(gfit2)
common <- match(names(f1), names(f2))
stopifnot(max(abs(f1 - f2[common])) <= 256*.Machine$double.eps)
stopifnot(all(is.na(f2[-common])))

r1 <- resid(gfit)
r2 <- resid(gfit2)
common <- match(names(r1), names(r2))
stopifnot(max(abs(r1 - r2[common])) <= 256*.Machine$double.eps)
stopifnot(all(is.na(r2[-common])))

p1 <- predict(gfit)
p2 <- predict(gfit2)
common <- match(names(p1), names(p2))
stopifnot(max(abs(p1 - p2[common])) <= 256*.Machine$double.eps)
stopifnot(all(is.na(p2[-common])))

## tests of diagnostic measures.
set.seed(11)
x <- 1:10
y <- c(rnorm(9),NA)
fit <- lm(y ~ x, na.action=na.exclude)
fit2 <- lm(y ~ x, subset=-10)

lm.influence(fit2); lm.influence(fit)

rstandard(fit2); rstandard(fit)
rstudent(fit2); rstudent(fit)

dffits(fit2); dffits(fit)

dfbetas(fit2); dfbetas(fit)

covratio(fit2); covratio(fit)

cooks.distance(fit2); cooks.distance(fit)

(inf <- influence.measures(fit))
(inf2 <- influence.measures(fit2))

summary(inf)
summary(inf2)

plot(fit)