---
title: "Statistical test workflows"
output: rmarkdown::pdf_document
vignette: >
  %\VignetteIndexEntry{Statistical test workflows}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, echo = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = TRUE)
```

## Introduction

`testflow` provides statistical testing workflows organized by study design.

## One numerical variable

```{r}
library(testflow)
cardio <- make_cardio_data()
test_one_sample(cardio, sbp_3m, mu = 140)
```

## Two independent groups

```{r}
test_two_groups(sbp_3m ~ sex, data = cardio)
```

## Paired measurements

```{r}
test_paired(sbp_3m ~ sbp_baseline, data = cardio)
```

## More than two groups

```{r}
test_groups(sbp_3m ~ treatment, data = cardio)
```

## Factorial designs

```{r}
test_factorial(sbp_3m ~ sex * treatment, data = cardio)
```

## Repeated measurements

```{r}
test_repeated(cardio, c(sbp_baseline, sbp_3m, sbp_6m), id = id)
```

The repeated numeric workflow chooses repeated-measures ANOVA when the
within-time normality checks are acceptable and Friedman otherwise. Post-hoc
comparisons are paired t-tests for the parametric branch and paired Wilcoxon
tests for the non-parametric branch.

## Categorical outcomes

```{r}
test_categorical(treatment ~ controlled_3m, data = cardio)
```

## Repeated categorical outcomes

```{r}
test_repeated_categorical(cardio, c(controlled_baseline, controlled_3m, controlled_6m))
```

The repeated categorical workflow uses Cochran Q for binary repeated outcomes
and pairwise McNemar tests for follow-up comparisons.

## References

- Fisher, R. A. (1925). \emph{Statistical Methods for Research Workers}.
- Gosset, W. S. (1908). The probable error of a mean.
- Welch, B. L. (1947). Generalization of Student's problem with unequal variances.
- Wilcoxon, F. (1945). Individual comparisons by ranking methods.
- Mann, H. B., & Whitney, D. R. (1947). On a test of whether one of two random variables is stochastically larger than the other.
- Levene, H. (1960). Robust tests for equality of variances.
- Kruskal, W. H., & Wallis, W. A. (1952). Use of ranks in one-criterion variance analysis.
- Tukey, J. W. (1949). Comparing individual means in the analysis of variance.
- Dunn, O. J. (1964). Multiple comparisons using rank sums.
- Friedman, M. (1937). The use of ranks to avoid the assumption of normality implicit in the analysis of variance.
- Cochran, W. G. (1950). The comparison of percentages in matched samples.
- McNemar, Q. (1947). Note on the sampling error of the difference between correlated proportions or percentages.
- Pearson, K. (1895, 1900).
- Spearman, C. (1904). The proof and measurement of association between two things.
- Kendall, M. G. (1938). A new measure of rank correlation.
- Cramer, H. (1946). \emph{Mathematical Methods of Statistics}.
- Clopper, C. J., & Pearson, E. S. (1934). The use of confidence or fiducial limits illustrated in the case of the binomial.
- Cohen, J. (1988). \emph{Statistical Power Analysis for the Behavioral Sciences}.
- Draper, N. R., & Smith, H. (1998). \emph{Applied Regression Analysis} (3rd ed.).
- Hosmer, D. W., Lemeshow, S., & Sturdivant, R. X. (2013). \emph{Applied Logistic Regression} (3rd ed.).
- Kaplan, E. L., & Meier, P. (1958). Nonparametric estimation from incomplete observations.
- Mantel, N. (1966). Evaluation of survival data and two new rank order statistics arising in its consideration.
- Cox, D. R. (1972). Regression models and life-tables.
- Grambsch, P. M., & Therneau, T. M. (1994). Proportional hazards tests and diagnostics based on weighted residuals.
- Simel, D. L., Samsa, G. P., & Matchar, D. B. (1991). Likelihood ratios with confidence: sample size estimation for diagnostic test studies.
- Altman, D. G., & Bland, J. M. (1994). Diagnostic tests 1: sensitivity and specificity.
- Hanley, J. A., & McNeil, B. J. (1982). The meaning and use of the area under a receiver operating characteristic (ROC) curve.
- Youden, W. J. (1950). Index for rating diagnostic tests.
- Fleiss, J. L., Cohen, J., & Everitt, B. S. (1969). Large sample standard errors of kappa and weighted kappa.
- Landis, J. R., & Koch, G. G. (1977). The measurement of observer agreement for categorical data.
- Shrout, P. E., & Fleiss, J. L. (1979). Intraclass correlations: uses in assessing rater reliability.
- McGraw, K. O., & Wong, S. P. (1996). Forming inferences about some intraclass correlation coefficients.
- Koo, T. K., & Li, M. Y. (2016). A guideline of selecting and reporting intraclass correlation coefficients for reliability research.

## Correlation

```{r}
test_correlation(sbp_3m ~ age, data = cardio)
```

## Regression

```{r}
test_linear_regression(sbp_3m ~ age + ldl, data = cardio)
```

Linear regression reports the overall model F-test as the primary result,
the per-term coefficient table (with confidence intervals) in
`alternative_tests$coefficients`, and R squared / adjusted R squared as the
effect size. Residual normality (Shapiro-Wilk), homoscedasticity
(Breusch-Pagan), and multicollinearity (variance inflation factor, when
there is more than one predictor) are checked as assumptions.

```{r}
test_logistic_regression(controlled_3m ~ age + ldl, data = cardio)
```

Logistic regression reports the likelihood-ratio test against the
intercept-only model as the primary result, coefficients on both the
log-odds and odds-ratio scale in `alternative_tests`, and McFadden's pseudo R
squared as the effect size. Multicollinearity and influential observations
(Cook's distance) are checked as assumptions.

## Survival analysis

`make_cardio_data()` does not include time-to-event data, so this section
uses a small synthetic example instead.

```{r}
set.seed(1)
n <- 100
survival_dat <- tibble::tibble(
  time = rexp(n, 0.1),
  status = rbinom(n, 1, 0.7),
  arm = rep(c("control", "treatment"), each = n / 2),
  age = rnorm(n, 60, 10)
)
test_survival(Surv(time, status) ~ arm, data = survival_dat)
```

`test_survival()` reports the log-rank test as the primary result and a
companion univariate Cox hazard ratio as the effect size; the hazard ratio
is not itself part of the log-rank test and can disagree with it under
strongly non-proportional hazards. `Surv()` is re-exported from the
`survival` package, so it is available after `library(testflow)` alone.

```{r}
test_cox(Surv(time, status) ~ age + arm, data = survival_dat)
```

`test_cox()` reports the overall likelihood-ratio test as the primary
result, hazard ratios per term, the concordance index as the effect size,
and checks the proportional-hazards assumption via the Schoenfeld residual
test (`survival::cox.zph()`).

## Diagnostic test accuracy and agreement

`make_cardio_data()` has no diagnostic-test or multi-rater columns, so this
section uses small synthetic examples instead.

```{r}
set.seed(1)
diag_dat <- tibble::tibble(
  test = c(rep("positive", 55), rep("negative", 98)),
  reference = c(rep("positive", 45), rep("negative", 10), rep("positive", 8), rep("negative", 90))
)
test_diagnostic(diag_dat, test, reference)
```

`test_diagnostic()` reports sensitivity, specificity, predictive values, and
likelihood ratios (each with an exact confidence interval) in
`alternative_tests$diagnostic_table`, and tests overall accuracy against the
no-information rate as the primary result, following the same convention as
`caret::confusionMatrix()`.

```{r}
roc_dat <- tibble::tibble(
  marker = c(rnorm(60, 2, 1), rnorm(50, 0, 1)),
  disease = c(rep("yes", 60), rep("no", 50))
)
test_roc(roc_dat, marker, disease)
```

`test_roc()` computes the AUC from the Mann-Whitney relationship, a
closed-form Hanley-McNeil confidence interval, and the Youden's-J-optimal
threshold.

```{r}
agree_dat <- tibble::tibble(
  rater1 = sample(c("mild", "moderate", "severe"), 100, replace = TRUE),
  rater2 = sample(c("mild", "moderate", "severe"), 100, replace = TRUE)
)
test_agreement(agree_dat, rater1, rater2)
```

`test_agreement()` reports Cohen's kappa with the Fleiss-Cohen-Everitt
(1969) large-sample confidence interval.

```{r}
icc_dat <- tibble::tibble(
  rater1 = rnorm(30, 50, 10),
  rater2 = rnorm(30, 50, 10),
  rater3 = rnorm(30, 50, 10)
)
test_icc(icc_dat, c(rater1, rater2, rater3))
```

`test_icc()` reports ICC(2,1) (two-way random effects, absolute agreement,
single measurement) as the primary result, following the reliability-study
recommendation of Koo & Li (2016), alongside ICC(1,1) and ICC(3,1) for
comparison in `alternative_tests$icc_table`.

## Outliers

```{r}
test_outliers(c(sbp_3m, ldl, crp), data = cardio)
```

## Reporting and plotting

Every workflow returns a `testflow` object. Use `report(x)`, `plot(x)`, and `as_tibble(x)`.
See `effect-size-formulas.Rmd` for the exact formulas used by the reported
effect-size estimates.
