---
title: "Diagnostics"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Diagnostics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

# car, lmtest and moments back individual assumption tests; gridExtra
# arranges the dashboard. All are in Suggests.
has_assumption_pkgs <- requireNamespace("car", quietly = TRUE) &&
  requireNamespace("lmtest", quietly = TRUE)
has_gridextra <- requireNamespace("gridExtra", quietly = TRUE)
```

```{r setup}
library(tidylearn)
library(dplyr)
```

## Overview

A model that fits is not the same as a model you should use. These functions
answer four questions:

- **Do the assumptions hold?** `tl_check_assumptions()`
- **Which observations drove the fit?** `tl_influence_measures()`,
  `tl_detect_outliers()`
- **Is the difference between two models real?** `tl_compare_cv()`,
  `tl_test_model_difference()`
- **Are there interactions I have not modelled?** `tl_test_interactions()`,
  `tl_interaction_effects()`

`tl_explore()` runs an unsupervised sweep over a dataset before you model it
at all.

We use a linear model throughout, because that is where assumption checking
has teeth.

```{r}
model <- tl_model(mtcars, mpg ~ wt + hp + disp, method = "linear")
```

## Checking Assumptions

`tl_check_assumptions()` runs six checks and returns a verdict on each.

```{r, eval = has_assumption_pkgs}
assumptions <- tl_check_assumptions(model, verbose = FALSE)
names(assumptions)
```

```{r, eval = has_assumption_pkgs}
assumptions$overall
```

Each entry carries the test that was run, the verdict, and what to do about
it. Nothing here is a pass/fail gate — the recommendation is a prompt, not an
instruction.

```{r, eval = has_assumption_pkgs}
assumptions$normality
```

```{r, eval = has_assumption_pkgs}
assumptions$multicollinearity
```

A compact table of every check:

```{r, eval = has_assumption_pkgs}
checks <- c("linearity", "independence", "homoscedasticity",
            "normality", "multicollinearity", "outliers")

data.frame(
  assumption = vapply(checks, function(x) assumptions[[x]]$assumption,
                      character(1)),
  holds = vapply(checks, function(x) isTRUE(assumptions[[x]]$check),
                 logical(1)),
  detail = vapply(checks, function(x) assumptions[[x]]$details, character(1)),
  row.names = NULL
)
```

`disp` correlating with both `wt` and `hp` is what drives the VIF here, and it
is the kind of thing that is invisible in a coefficient table.

### The dashboard

`tl_diagnostic_dashboard()` draws the standard panels in one grid.

```{r, eval = has_assumption_pkgs && has_gridextra, fig.height = 7}
tl_diagnostic_dashboard(model)
```

Switch off any section you do not want with `include_influence`,
`include_assumptions` or `include_performance`.

## Influence

`tl_influence_measures()` returns one row per observation with Cook's
distance, leverage, DFFITS, standardised and studentised residuals, DFBETAS
per coefficient, and a flag for each.

```{r}
influence <- tl_influence_measures(model)
dim(influence)
```

```{r}
influence %>%
  filter(is_influential) %>%
  select(observation, cooks_distance, leverage, dffits, std_residual)
```

The flags use conventional cutoffs, which you can override with
`threshold_cook`, `threshold_leverage` and `threshold_dffits`.

The DFBETAS columns say *which coefficient* an observation moved, which is
usually the more useful question:

```{r}
influence %>%
  select(observation, starts_with("dfbetas_")) %>%
  arrange(desc(abs(dfbetas_wt))) %>%
  head(4)
```

### Refit without the influential rows

The point of the exercise is to see whether the conclusion survives.

```{r}
keep <- !influence$is_influential
refit <- tl_model(mtcars[keep, ], mpg ~ wt + hp + disp, method = "linear")

data.frame(
  term = names(coef(model$fit)),
  all_rows = round(unname(coef(model$fit)), 4),
  without_influential = round(unname(coef(refit$fit)), 4)
)
```

```{r}
sum(!keep)
```

If dropping a handful of rows moves a coefficient materially, that
coefficient describes those rows rather than the population you sampled.

## Outliers in the Data

`tl_influence_measures()` is about a fitted model. `tl_detect_outliers()`
works on the data itself, before or independently of any fit.

```{r}
outliers <- tl_detect_outliers(
  mtcars,
  variables = c("mpg", "hp", "wt"),
  method = "iqr",
  plot = FALSE
)

outliers$outlier_counts$total
outliers$outlier_counts$by_variable
```

```{r}
mtcars[outliers$outlier_indices, c("mpg", "hp", "wt")]
```

`method` also takes `"zscore"` and `"mahalanobis"`. The first two treat each
variable separately; Mahalanobis distance accounts for the correlation
between them, so it finds points that are unremarkable on every single axis
and unusual in combination.

```{r}
mahal <- tl_detect_outliers(
  mtcars,
  variables = c("mpg", "hp", "wt"),
  method = "mahalanobis",
  plot = FALSE
)

mahal$outlier_indices
```

Set `plot = TRUE` to get a ggplot2 object back in `$plot`.

## Comparing Models

A difference in a single held-out score is not evidence. `tl_compare_cv()`
scores several fitted models over the same folds.

```{r}
simple <- tl_model(mtcars, mpg ~ wt, method = "linear")
full <- tl_model(mtcars, mpg ~ wt + hp + disp, method = "linear")
tree <- tl_model(mtcars, mpg ~ wt + hp + disp, method = "tree")

cv <- tl_compare_cv(
  mtcars,
  models = list(simple = simple, full = full, tree = tree),
  folds = 5,
  metrics = c("rmse", "rsq")
)

names(cv)
```

```{r}
cv$summary
```

Per-fold scores are kept as well, which is what makes a test possible:

```{r}
head(cv$fold_metrics)
```

### Is the difference real?

`tl_test_model_difference()` compares each model against a baseline using the
per-fold scores.

```{r}
tl_test_model_difference(
  cv,
  baseline_model = "simple",
  metric = "rmse",
  test = "t.test"
)
```

With five folds this has very little power, so treat a non-significant result
as "these folds do not separate the models" rather than as evidence they are
equivalent. `test = "wilcox.test"` drops the normality assumption, which
matters more at small fold counts than the loss of power costs you.

## Interactions

`tl_test_interactions()` fits each candidate interaction and reports whether
it earns its degrees of freedom.

```{r}
interactions <- tl_test_interactions(
  mtcars, mpg ~ wt + hp + disp,
  all_pairs = TRUE
)

interactions
```

`delta_r2` is the more useful column: a p-value tells you the term is
detectable, `delta_r2` tells you whether it is worth carrying.

Restrict the search with `numeric_only`, `categorical_only` or `mixed_only`,
or name a single pair with `var1` and `var2`.

### Reading an interaction

Once a term is in the model, `tl_interaction_effects()` says what it does at
different levels of the moderator.

```{r}
model_int <- tl_model(mtcars, mpg ~ wt * hp, method = "linear")
effects <- tl_interaction_effects(model_int, var = "wt", by_var = "hp")

effects$slopes
```

The slope of `mpg` on `wt` weakens as `hp` rises — extra weight costs less
fuel economy in a high-powered car, which already had little to lose.

`slope_se` describes the straight line fitted to the prediction grid rather
than the sampling uncertainty of the marginal effect — for a linear model the
grid is exactly linear, so it is near zero by construction. Use
`summary(model_int$fit)` for inference on the interaction coefficient.

```{r}
summary(model_int$fit)$coefficients
```

`tl_auto_interactions()` does the search and the refit in one step, returning
a model with the surviving interactions already in the formula:

```{r}
auto <- tl_auto_interactions(mtcars, mpg ~ wt + hp + disp)
auto$spec$formula
```

## Exploring Before Modelling

`tl_explore()` runs PCA, picks a cluster count, clusters, and computes a
distance summary in one call. It is a first look at a dataset, not a
diagnostic of a fit.

```{r}
eda <- tl_explore(iris, response = "Species", max_components = 4, k_range = 2:5)
names(eda)
```

```{r}
eda$optimal_k
```

```{r}
get_pca_variance(eda$pca)
```

```{r, fig.height = 6}
plot(eda)
```

## A Checklist

For a linear model, in order:

1. `tl_check_assumptions()` — six checks, with the reason each one failed.
2. `tl_influence_measures()` — refit without the flagged rows and see whether
   the coefficients hold.
3. `tl_test_interactions()` — the effect you assumed was additive may not be.
4. `tl_compare_cv()` then `tl_test_model_difference()` — before preferring
   one model over another.

For tree-based and other non-parametric methods, steps 1 and 3 do not apply;
step 2 is available through `tl_detect_outliers()` on the data, and step 4
works unchanged.
