---
title: "Case-control genetic association study design"
author: "paweh package authors"
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
vignette: >
  %\VignetteIndexEntry{2. Case-Control Study Design}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# From a scientific contrast to a design

A case-control association study asks whether genotype distributions differ
between affected and unaffected groups. `paweh` supports the two-degree-of-
freedom genotype chi-square test and a one-degree-of-freedom genotype trend
test. The control-to-case ratio `k` connects the two sample sizes:

\[
N_{control}=kN_{case}.
\]

In model-based designs, prevalence, disease-allele frequency, relative risk,
and mode of inheritance determine the conditional genotype distributions. In
model-free designs, the researcher supplies those distributions directly as
`g1` (case or affected) and `g0` (control or unaffected).

# Basic power and inverse design

The unified interfaces are the normal starting point. This model-based example
asks for power with 500 cases and an equal number of controls, then solves the
inverse problem for 80% power.

```{r basic-cc}
basic_cc_power <- cc_power(
  N_case = 500, alpha = 0.05,
  input_mode = "model_based",
  prev = 0.05, pd = 0.30, R2 = 1.8, MOI = "M",
  k = 1, verbose = FALSE
)

basic_cc_mssn <- cc_mssn(
  power = 0.80, alpha = 0.05,
  input_mode = "model_based",
  prev = 0.05, pd = 0.30, R2 = 1.8, MOI = "M",
  k = 1, verbose = FALSE
)

data.frame(
  quantity = c("Genotype-test power", "Required cases", "Required controls", "Required total"),
  value = c(
    basic_cc_power$tests$genotypes$power,
    basic_cc_mssn$tests$genotypes$MSSN_case,
    basic_cc_mssn$tests$genotypes$MSSN_ctrl,
    basic_cc_mssn$tests$genotypes$MSSN_total
  )
)
```

These are asymptotic operating characteristics under the stated disease
model. They do not account for uncertainty in the input parameters.

# Published-example reproduction: Gordon et al. (2002)

Gordon et al. studied case-control genotype tests under genotype
misclassification. Their Figure 1b DSB example assumes case allele frequency
0.40, control allele frequency 0.50, Hardy--Weinberg equilibrium (HWE), equal
numbers of cases and controls, 95% power, and `alpha = 0.05`.

Under HWE, allele frequency `p` gives genotype probabilities
`c((1-p)^2, 2*p*(1-p), p^2)`:

```{r gordon-setup}
g_gordon_case <- c(0.36, 0.48, 0.16)
g_gordon_control <- c(0.25, 0.50, 0.25)
```

In `paweh`'s two-parameter error model, Gordon's DSB parameters map as
`e1 = delta` and `e2 = epsilon / 2`. The division by two reflects equal
allocation of total heterozygote error to the two homozygotes.

```{r gordon-reproduction}
gordon_no_error <- cc_mssn(
  power = 0.95, alpha = 0.05,
  input_mode = "model_free",
  g1 = g_gordon_case, g0 = g_gordon_control, k = 1,
  geno_misclass = "2p", e1 = 0, e2 = 0,
  verbose = FALSE
)

gordon_dsb_10 <- cc_mssn(
  power = 0.95, alpha = 0.05,
  input_mode = "model_free",
  g1 = g_gordon_case, g0 = g_gordon_control, k = 1,
  geno_misclass = "2p", e1 = 0.10, e2 = 0.05,
  verbose = FALSE
)

data.frame(
  DSB_error = c("No error", "delta = epsilon = 0.10"),
  published_required_cases = c("approximately 386", "477"),
  paweh_required_cases = c(
    gordon_no_error$tests$genotypes$MSSN_case,
    gordon_dsb_10$tests$genotypes$MSSN_case
  )
)
```

This is a **Published-example reproduction**. The error result is exact at
477 cases. The baseline is near-exact: `paweh` returns 387 rather than the
paper's approximately 386 because the underlying continuous MSSN is about
386.08 and `paweh` ceilings the number of subjects.

## Extending the same design: how MSSN grows with DSB error

The following sweep keeps the published genotype distributions and design
fixed. A multiplier of one is `delta = epsilon = 0.10`; zero is no error.

```{r gordon-plot, fig.cap="Required case sample size as the Gordon DSB error pattern increases from no error to delta = epsilon = 0.10."}
plot_cc_mssn(
  x_var = "geno_error_multiplier",
  x_values = seq(0, 1, length.out = 6),
  test = "genotypes", input_mode = "model_free",
  sample_size = "case",
  power = 0.95, alpha = 0.05,
  g1 = g_gordon_case, g0 = g_gordon_control, k = 1,
  geno_misclass = "2p", e1_base = 0.10, e2_base = 0.05,
  x_label = "Fraction of the 10%/10% DSB error pattern",
  title = "Required cases under the Gordon DSB error pattern"
)
```

The curve is a design sensitivity analysis: it treats error rates as fixed
scenarios rather than estimating their uncertainty.

# Published-example reproduction: Edwards et al. (2005)

Edwards et al. examined phenotype error in a biallelic case-control design.
Their Figure 2 example uses 250 cases, 250 controls, `alpha = 0.01`, affected
minor-allele frequency 0.05, unaffected minor-allele frequency 0.15, and HWE.

```{r edwards-setup}
g_edwards_affected <- c(0.9025, 0.0950, 0.0025)
g_edwards_unaffected <- c(0.7225, 0.2550, 0.0225)

edwards_scenarios <- data.frame(
  prev = c(0.05, 0.01, 0.05, 0.01),
  theta = 0,
  phi = c(0.01, 0.01, 0.02, 0.02),
  published_power = c(0.91, 0.33, 0.76, 0.11)
)

edwards_scenarios$paweh_power <- vapply(seq_len(nrow(edwards_scenarios)), function(i) {
  cc_power(
    N_case = 250, alpha = 0.01, k = 1,
    input_mode = "model_free",
    g1 = g_edwards_affected, g0 = g_edwards_unaffected,
    prev = edwards_scenarios$prev[i],
    pheno_misclass = TRUE,
    theta = edwards_scenarios$theta[i],
    phi = edwards_scenarios$phi[i],
    verbose = FALSE
  )$tests$genotypes$power
}, numeric(1))

edwards_scenarios
```

This is an excellent **Published-example reproduction**; the paper reports
rounded powers. Here `theta = Pr(affected -> control)` and
`phi = Pr(unaffected -> case)`. For a rare disease, even a small `phi` can
matter greatly because the unaffected source population is much larger than
the affected source population. The no-error reference power is about
0.9896.

## Extending Edwards: prevalence and false-positive classification

Both plots use the same underlying genotype distributions. Only prevalence is
changed between panels.

```{r edwards-plot-common, echo=FALSE, results='hide'}
edwards_phi <- seq(0, 0.02, by = 0.0025)
```

```{r edwards-plot-prev05, fig.cap="Power as unaffected subjects are increasingly classified as cases, with prevalence 0.05."}
plot_cc_power(
  x_var = "phi", x_values = edwards_phi,
  test = "genotypes", input_mode = "model_free",
  N_case = 250, alpha = 0.01, k = 1,
  g1 = g_edwards_affected, g0 = g_edwards_unaffected,
  prev = 0.05, pheno_misclass = TRUE, theta = 0,
  title = "Power vs phenotype misclassification (prevalence 0.05)"
)
```

```{r edwards-plot-prev01, fig.cap="The same sensitivity analysis with prevalence 0.01."}
plot_cc_power(
  x_var = "phi", x_values = edwards_phi,
  test = "genotypes", input_mode = "model_free",
  N_case = 250, alpha = 0.01, k = 1,
  g1 = g_edwards_affected, g0 = g_edwards_unaffected,
  prev = 0.01, pheno_misclass = TRUE, theta = 0,
  title = "Power vs phenotype misclassification (prevalence 0.01)"
)
```

# Locus heterogeneity and the trend test

For case-control calculations, `pi` is the fraction of cases attributable to
the modeled locus; `1 - pi` is the heterogeneous fraction. The model-free
example below compares genotype and trend-test MSSN as heterogeneity grows.

```{r cc-heterogeneity-plot, fig.cap="Required cases for genotype and trend tests as the heterogeneous fraction increases."}
plot_cc_mssn(
  x_var = "locus_het_rate",
  x_values = seq(0, 0.50, by = 0.10),
  input_mode = "model_free", compare_tests = TRUE,
  sample_size = "case",
  power = 0.80, alpha = 0.05, k = 1,
  g1 = g_gordon_case, g0 = g_gordon_control,
  title = "Required cases vs locus heterogeneity"
)
```

Ahn et al. (2007) provide **Published methodological context** for how SNP
genotyping errors affect Cochran--Armitage trend-test designs. The plot above
does not claim to reproduce their specific Taylor-series error framework.

# Differential genotype error: a Type I error warning

Moskvina et al. (2006) showed that different error mechanisms in cases and
controls can substantially inflate Type I error, even when absolute error
rates are small. The following calculation is therefore context, not a
calibrated false-positive analysis.

```{r differential-error}
diff_error_design <- cc_power(
  N_case = 500, alpha = 0.05, k = 1,
  input_mode = "model_free",
  g1 = g_gordon_case, g0 = g_gordon_control,
  geno_misclass = "diff3p",
  case_e01 = 0.02, case_e02 = 0.01, case_e03 = 0.005,
  ctrl_e01 = 0.01, ctrl_e02 = 0.005, ctrl_e03 = 0.002,
  verbose = FALSE
)
diff_error_design$tests$genotypes$power
```

**Important:** `paweh` returns nominal asymptotic power or MSSN under
differential genotype error using the usual chi-square critical value. It does
not independently recalibrate an arbitrary differential-error null
distribution. These nominal operating characteristics must not be interpreted
as protection against the Type I error inflation described by Moskvina et al.

# Independent sensitivity scenarios

The unified interface accepts several modifiers in one call and reports each
as a separate sensitivity scenario from the same no-error baseline. It does
not combine the modifier transformations.

```{r independent-cc-scenarios}
scenario_cc <- cc_mssn(
  power = 0.80, alpha = 0.05, k = 1,
  input_mode = "model_free",
  g1 = g_gordon_case, g0 = g_gordon_control,
  locus_het = TRUE, pi = 0.80,
  prev = 0.05, pheno_misclass = TRUE, theta = 0.02, phi = 0.01,
  geno_misclass = "2p", e1 = 0.02, e2 = 0.01,
  verbose = FALSE
)

lapply(scenario_cc$scenarios, function(scenario) {
  c(
    genotype_cases = scenario$tests$genotypes$MSSN_case,
    trend_cases = scenario$tests$trend$MSSN_case
  )
})
```

# Choosing an interface and recognizing limits

Use `cc_power()` and `cc_mssn()` for most designs, particularly when switching
between model-based and model-free inputs or comparing modifier scenarios. Specialized
genotype-only functions remain available in the reference manual for focused
calculations.

Interpretation depends on assumptions close to each calculation:

- HWE was used to construct genotype probabilities in the published examples.
- Model-free probabilities are treated as supplied truth, not estimated with
  uncertainty.
- Model-based results inherit the selected prevalence, risk, and inheritance
  model.
- `theta` and `phi` describe opposite phenotype-error directions.
- `pi` is the homogeneous fraction; `1 - pi` is heterogeneous.
- Genotype and trend calculations use asymptotic chi-square approximations.
- Integer MSSNs are ceiling values; a paper using another rounding convention
  can differ by one subject.
- Differential genotype error can invalidate nominal Type I error calibration.

# References

Ahn K, Haynes C, Kim W, St Fleur R, Gordon D, Finch SJ. The effects of SNP
genotyping errors on the power of the Cochran-Armitage linear trend test for
case/control association studies. *Annals of Human Genetics*. 2007;71(Pt
2):249--261. <https://doi.org/10.1111/j.1469-1809.2006.00318.x>.

Edwards BJ, Haynes C, Levenstien MA, Finch SJ, Gordon D. Power and sample size
calculations in the presence of phenotype errors for case/control genetic
association studies. *BMC Genetics*. 2005;6:18.
<https://doi.org/10.1186/1471-2156-6-18>.

Gordon D, Finch SJ, Nothnagel M, Ott J. Power and sample size calculations for
case-control genetic association tests when errors are present: application to
single nucleotide polymorphisms. *Human Heredity*. 2002;54(1):22--33.
<https://doi.org/10.1159/000066696>.

Moskvina V, Craddock N, Holmans P, Owen MJ, O'Donovan MC. Effects of
differential genotyping error rate on the type I error probability of
case-control studies. *Human Heredity*. 2006;61(1):55--64.
<https://doi.org/10.1159/000092553>.

# Session information

```{r session-info}
sessionInfo()
```