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

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

# Quantitative-trait study design in PAWEH

For a quantitative trait, genotype can shift the phenotype distribution rather
than simply changing a binary disease probability. `paweh` uses a
Falconer-style model to calculate prospective operating characteristics for
association-study design. Its core inputs are:

- `qtl_var`: the proportion of trait variance attributable to the modeled QTL;
- `tau`: dominance relative to additivity;
- `pd`: the frequency of the modeled allele; and
- residual variance or covariance: the remaining phenotype variation and, for
  multiple traits, how those traits co-vary.

These functions plan studies under a specified model. They do not estimate QTL
effects from individual-level phenotype and genotype observations.

# Genotype-specific trait distributions

The three genotype groups can have different expected trait distributions, but
residual variation makes those distributions overlap. Greater separation of
the genotype means generally creates a stronger potential association signal.

```{r genotype-distributions, fig.cap="Deterministic genotype-specific quantitative-trait distributions under the Falconer model."}
plot_qtl_genotype_distribution(
  qtl_var = 0.25,
  tau = 0.25,
  pd = 0.25,
  type = "binned",
  scale = "density"
)
```

Vertical lines identify the genotype-specific means. `type = "binned"` is a
deterministic histogram-style representation: each bar is calculated from the
difference between analytic Normal CDF values at its bin boundaries. No random
observations are simulated.

# Continuous-trait ANOVA design

When every recruited subject is measured on the continuous trait, a one-way
ANOVA can compare the three genotype groups. The following example calculates
power at a fixed sample size and the minimum total sample size (MSSN) for 80%
power.

```{r anova-design}
anova_power <- qtl_anova_power(
  N = 996, alpha = 0.0001,
  qtl_var = 0.025, tau = 0.5, pd = 0.15,
  count_method = "rounded", verbose = FALSE
)

anova_mssn <- qtl_anova_mssn(
  power = 0.80, alpha = 0.0001,
  qtl_var = 0.025, tau = 0.5, pd = 0.15,
  count_method = "rounded", multiple_of_three = TRUE,
  verbose = FALSE
)

data.frame(
  quantity = c("Power at N = 996", "Required total N", "Power at required N"),
  value = c(anova_power$power, anova_mssn$N, anova_mssn$achieved_power)
)
```

One power curve summarizes how the design changes as more participants are
measured.

```{r anova-plot, fig.cap="One-way genotype-group ANOVA power as total sample size increases."}
plot_qtl_anova_power(
  x_var = "N",
  x_values = seq(300, 1200, by = 150),
  alpha = 0.0001,
  qtl_var = 0.025,
  tau = 0.5,
  pd = 0.15
)
```

# Extreme-phenotype sampling

Instead of recruiting across the full phenotype distribution, investigators
may select subjects from its extremes. In `paweh`'s single-trait threshold
design, upper-tail individuals become selected cases, lower-tail individuals
become selected controls, and individuals in the middle are not selected.
Upper and lower events are distinct tails rather than complements.

Li et al. (2019) provide a **motivating published application** of extreme
quantitative-phenotype sampling. Their procedure is not the threshold-selected
genotype chi-square method implemented by `paweh`, so this is motivation rather
than a numerical reproduction.

```{r threshold-parameters}
threshold_model <- qtl_falconer_threshold_parameters(
  qtl_var = 0.025, tau = 0.5, pd = 0.15,
  x_upper = 5, x_lower = 5,
  verbose = FALSE
)

threshold_power <- qtl_threshold_chisq_power(
  N_case = 126, alpha = 0.0001,
  qtl_var = 0.025, tau = 0.5, pd = 0.15,
  x_upper = 5, x_lower = 5, k = 1,
  verbose = FALSE
)

threshold_mssn <- qtl_threshold_chisq_mssn(
  power = 0.80, alpha = 0.0001,
  qtl_var = 0.025, tau = 0.5, pd = 0.15,
  x_upper = 5, x_lower = 5, k = 1,
  verbose = FALSE
)

data.frame(
  quantity = c(
    "Upper threshold", "Lower threshold", "Power with 126 selected cases",
    "Required selected cases", "Required selected controls",
    "Expected source population screened for cases",
    "Expected source population screened for controls"
  ),
  value = c(
    threshold_model$upper_threshold,
    threshold_model$lower_threshold,
    threshold_power$power,
    threshold_mssn$N_case,
    threshold_mssn$N_control,
    threshold_mssn$expected_population_screened_cases,
    threshold_mssn$expected_population_screened_controls
  )
)
```

The selected study MSSN and the expected source population screened are
different quantities. Screening burden can be much larger because only a
fraction of the source population falls into either selected tail.

```{r threshold-plot, fig.cap="Required selected cases as upper-tail sampling becomes less extreme, with the lower tail fixed at 5%."}
plot_qtl_threshold_chisq_mssn(
  x_var = "x_upper",
  x_values = c(2.5, 5, 7.5, 10, 15),
  sample_size = "case",
  power = 0.80,
  alpha = 0.0001,
  qtl_var = 0.025,
  tau = 0.5,
  pd = 0.15,
  x_lower = 5,
  k = 1,
  title = "Required cases vs upper-tail selection"
)
```

Selection stringency changes genotype enrichment, selected-sample MSSN, and
the operational burden of finding eligible participants.

# Multiple quantitative traits

Multiple correlated or biologically related phenotypes may carry complementary
information about the same locus. Yang et al. (2010) provide a **motivating
published application** of joint phenotype analysis, but their method is not
the Pillai framework implemented by `paweh` and is not reproduced here.

For multivariate `paweh` designs, `qtl_var` and `tau` become
phenotype-specific vectors, `pd` remains shared, and `cor_matrix` describes
phenotype correlation. Use `test = "pillai"` for a continuous multivariate
design or `test = "threshold_chisq"` for jointly threshold-selected phenotypes.

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

The published pleiotropic example uses allele frequency 0.05, QTL variances
0.10 and 0.05, dominance ratios 0 and 0.50, independent phenotypes, target
power 0.80, and genome-wide `alpha = 5e-8`.

```{r gordon-multivariate-reproduction}
gordon_mv <- qtl_multivariate_mssn_full(
  power = 0.80,
  alpha = 5e-8,
  qtl_var = c(0.10, 0.05),
  tau = c(0, 0.50),
  pd = 0.05,
  cor_matrix = diag(2),
  test = "pillai",
  verbose = FALSE
)

stopifnot(gordon_mv$N == 326L)

data.frame(
  quantity = c(
    "Published integer MSSN",
    "paweh integer MSSN",
    "Historical fractional MSSN"
  ),
  value = c(326, gordon_mv$N, gordon_mv$historical_fractional_mssn)
)
```

The published result and `paweh` both require 326 participants. The historical
continuous root is approximately 325.5057; recruitment requires an integer, so
the result is rounded upward.

The inverse calculation evaluates power at the returned integer design using
the same model.

```{r gordon-multivariate-power}
gordon_mv_power <- qtl_multivariate_power_full(
  N = gordon_mv$N,
  alpha = 5e-8,
  qtl_var = c(0.10, 0.05),
  tau = c(0, 0.50),
  pd = 0.05,
  cor_matrix = diag(2),
  test = "pillai",
  verbose = FALSE
)

gordon_mv_power$power
```

# Joint multivariate threshold selection

For `test = "threshold_chisq"`, affected subjects must satisfy both upper-tail
conditions, and unaffected subjects must satisfy both lower-tail conditions:

\[
\text{affected}: Y_1 \ge T_{U1} \;\mathbf{AND}\; Y_2 \ge T_{U2},
\]

\[
\text{unaffected}: Y_1 \le T_{L1} \;\mathbf{AND}\; Y_2 \le T_{L2}.
\]

Subjects in mixed or middle regions are unselected. For example,
`x_upper = c(10, 10)` and `x_lower = c(10, 10)` define joint upper- and
lower-tail selection at 10% for each phenotype. These thresholds illustrate
the design and are not part of the published Gordon reproduction.

# Scope and limitations

Joint multivariate tail probabilities use deterministic Miwa integration
through 20 dimensions and locally seeded Genz--Bretz integration above 20.
The latter restores the caller's random-number state and retains integration
status information.

Important scope boundaries are:

- calculations are prospective study-design calculations, not model fitting;
- genotype-specific phenotypes follow the modeled normal assumptions;
- residual variance and covariance determine distributional overlap;
- single-trait threshold designs exclude the middle of the distribution;
- multivariate threshold selection uses joint AND semantics;
- subjects in mixed or middle multivariate regions are unselected;
- a multivariate trend test is not implemented; and
- plots visualize the design model rather than observed participant-level data.

# References

Genz A, Bretz F. *Computation of Multivariate Normal and t Probabilities*.
Springer; 2009. <https://doi.org/10.1007/978-3-642-01689-9>.

Gordon D, Londono D, Patel P, Kim W, Finch SJ, Heiman GA. An analytic solution
to computation of power and sample size for genetic association studies under
a pleiotropic mode of inheritance. *Human Heredity*. 2017;81(4):194--209.
<https://doi.org/10.1159/000457135>.

Li Y, Levran O, Kim J, Zhang T, Chen X, Suo C. Extreme sampling design in
genetic association mapping of quantitative trait loci using balanced and
unbalanced case-control samples. *Scientific Reports*. 2019;9:15504.
<https://doi.org/10.1038/s41598-019-51790-w>.

Pillai KCS. Some new test criteria in multivariate analysis. *Annals of
Mathematical Statistics*. 1955;26(1):117--121.
<https://doi.org/10.1214/aoms/1177728599>.

Yang Q, Wu H, Guo C. Analyze multivariate phenotypes in genetic association
studies by combining univariate association tests. *Genetic Epidemiology*.
2010;34(5):444--454. <https://doi.org/10.1002/gepi.20497>.

# Session information

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