--- title: "Transmission disequilibrium test study design" author: "paweh package authors" output: BiocStyle::html_document: toc: true toc_float: true vignette: > %\VignetteIndexEntry{3. TDT 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) ``` # Why use an affected-child trio design? The transmission disequilibrium test (TDT) compares an allele transmitted to an affected child with the alternative parental allele that was not transmitted. This within-family comparison can avoid confounding from population stratification that affects unrelated case-control comparisons. Let `ET` and `ENT` denote expected transmitted and non-transmitted counts. The one-degree-of-freedom TDT non-centrality parameter is based on their contrast, \[ \lambda = \frac{(ET-ENT)^2}{ET+ENT}. \] In this vignette, `N` and every reported MSSN are numbers of affected-child trios, not numbers of individual people. # Power and MSSN for trios The canonical model-based interfaces answer inverse design questions under a common genetic model. ```{r basic-tdt} basic_tdt_power <- tdt_power( N = 600, input_mode = "model_based", pd = 0.30, prev = 0.05, R1 = 1.5, R2 = 2.25, alpha = 0.05, delta_prime = 1, misclass_rate = 0.01, heter_rate = 0.10, verbose = FALSE ) basic_tdt_mssn <- tdt_mssn( target_power = 0.80, input_mode = "model_based", pd = 0.30, prev = 0.05, R1 = 1.5, R2 = 2.25, alpha = 0.05, delta_prime = 1, misclass_rate = 0.01, heter_rate = 0.10, verbose = FALSE ) data.frame( scenario = c("No error", "Phenotype misclassification", "Locus heterogeneity"), power_at_600_trios = unlist(basic_tdt_power$power), required_trios = unlist(basic_tdt_mssn$N), row.names = NULL ) ``` The canonical backend reports separate modifier scenarios. It does not claim that the misclassification and heterogeneity components are a combined-error analysis. # Model-based and model-free TDT inputs Model-based input uses `pd`, `prev`, `R1`, `R2`, and `delta_prime`. The lower-level probability and count functions make the connection to `ET` and `ENT` visible without requiring a long derivation. ```{r expected-counts} gT <- tdt_expected_transmission_probability( pd = 0.30, prev = 0.05, R1 = 1.5, R2 = 2.25, delta_prime = 1, verbose = FALSE ) gNT <- tdt_expected_nontransmission_probability( pd = 0.30, prev = 0.05, R1 = 1.5, R2 = 2.25, delta_prime = 1, verbose = FALSE ) expected_counts <- tdt_expected_transmission_counts( N_star = 500, pd = 0.30, prev = 0.05, R1 = 1.5, R2 = 2.25, delta_prime = 1, pi = 1, verbose = FALSE ) data.frame( quantity = c("Transmission probability", "Non-transmission probability", "ET", "ENT"), value = c(gT$gT_star, gNT$gNT_star, expected_counts$ET_star, expected_counts$ENT_star) ) ``` If a pilot study or external calculation already supplies expected counts, they can be used directly: ```{r model-free-power} count_power <- tdt_power_from_expected_counts( ET = expected_counts$ET_star, ENT = expected_counts$ENT_star, alpha = 0.05 ) canonical_model_free <- tdt_power( N = 500, input_mode = "model_free", ET = expected_counts$ET_star, ENT = expected_counts$ENT_star, misclass_rate = 0, heter_rate = 0, verbose = FALSE ) c(lower_level_power = count_power$power, canonical_power = canonical_model_free$power$no_error) ``` For model-free MSSN, `n_trios` records the trio count represented by the supplied `ET` and `ENT`. Applying heterogeneity additionally requires `pd`; applying phenotype misclassification requires `pd` and `prev`. Model-free inputs therefore relocate assumptions rather than eliminate them. # Chapter 5 genotype misclassification The public `effect` argument can select `"none"`, `"phenotype_misclassification"`, `"locus_heterogeneity"`, or `"genotype_misclassification"`. The last option implements the Chapter 5 model-based TDT genotype-error model through the scalar argument `genotype_misclassification_rate`. For adjacent-genotype error rate (e), with (0 \le e \le 0.5), the transition matrix from true to observed genotype is \[ M(e) = \begin{bmatrix} 1-e & e & 0 \\ e & 1-2e & e \\ 0 & e & 1-e \end{bmatrix}. \] Thus a true genotype 0 can be observed as 0 or 1, a true genotype 1 can be observed as 0, 1, or 2, and a true genotype 2 can be observed as 1 or 2. Direct errors between genotypes 0 and 2 are not allowed. The calculation begins with the true Mendelian trio distribution, applies the transition model to the observed trio genotypes, and removes observed Mendelian-inconsistent states. It renormalizes the retained consistent states, calculates their expected transmitted and non-transmitted counts, and uses the resulting TDT non-centrality parameter for prospective power or MSSN. ```{r genotype-misclassification} genotype_error_power <- tdt_power( N = 600, input_mode = "model_based", pd = 0.30, prev = 0.05, R1 = 1.5, R2 = 2.25, alpha = 0.05, delta_prime = 1, effect = "genotype_misclassification", genotype_misclassification_rate = 0.01, verbose = FALSE ) genotype_error_mssn <- tdt_mssn( target_power = 0.80, input_mode = "model_based", pd = 0.30, prev = 0.05, R1 = 1.5, R2 = 2.25, alpha = 0.05, delta_prime = 1, effect = "genotype_misclassification", genotype_misclassification_rate = 0.01, verbose = FALSE ) genotype_error_mssn_values <- c( genotype_error_mssn$N$no_error, genotype_error_mssn$N$genotype_misclassification ) data.frame( scenario = c("No genotype error", "1% adjacent-genotype error"), power_at_600 = c( genotype_error_power$power$no_error, genotype_error_power$power$genotype_misclassification ), fractional_mssn = genotype_error_mssn_values, recruitment_trios = ceiling(genotype_error_mssn_values) ) ``` In this example, genotype error reduces power at 600 trios and increases the recruitment target for 80% power. This branch is available only with `input_mode = "model_based"` and uses one scalar adjacent-genotype error rate; it does not accept a custom transition matrix. It cannot be combined with phenotype misclassification or locus heterogeneity. It is neither TDTae nor the sequencing/NGS genotype-error implementation. # Published-example reproduction: Buyske et al. (2009) Buyske et al. examined phenotype misclassification in affected-child trios. The validated multiplicative example uses target power 0.80, `alpha = 1e-5`, `pd = 0.50`, prevalence 0.01, `R2 = 2.5`, `R1 = sqrt(2.5)`, maximum positive LD, and phenotype misclassification rate 0.05. ```{r buyske-reproduction} buyske <- tdt_mssn( target_power = 0.80, alpha = 1e-5, pd = 0.50, prev = 0.01, R1 = sqrt(2.5), R2 = 2.5, delta_prime = 1, misclass_rate = 0.05, heter_rate = 0, verbose = FALSE ) buyske_ratio <- buyske$N$misclassification / buyske$N$no_error data.frame( quantity = c("No-error MSSN", "Misclassification MSSN", "Final sample-size ratio"), published = c("approximately 546", "approximately 21,400", "approximately 39-fold"), paweh = c(buyske$N$no_error, buyske$N$misclassification, buyske_ratio) ) ``` This is an excellent **Published-example reproduction**: approximately 546 versus 21,401 affected-child trios, a 39.23-fold final sample-size ratio. A 39.23-fold ratio is a 3,822.7% increase, not a 39.23% increase. The `percent_increase` component stores this conventional percentage as `100 * ((new / baseline) - 1)`; the validated final-to-baseline ratio remains 39.23-fold. ## Extending Buyske: a misclassification curve ```{r buyske-curve, fig.cap="Required affected-child trios as phenotype misclassification increases under the Buyske model."} plot_tdt_mssn( x_var = "misclass_rate", x_values = seq(0, 0.05, by = 0.005), scenario = "misclassification", input_mode = "model_based", target_power = 0.80, alpha = 1e-5, pd = 0.50, prev = 0.01, R1 = sqrt(2.5), R2 = 2.5, delta_prime = 1, heter_rate = 0, title = "Required trios vs phenotype misclassification" ) ``` The steep curve identifies a fragile design region. It does not represent a confidence interval or a probability distribution over the error rate. ## Extending Buyske: prevalence by misclassification The optional interactive surface can evaluate the same Buyske parameters over a modest prevalence-by-error grid. The example is left unevaluated in the installed vignette to avoid embedding a multi-megabyte Plotly widget; run the copyable code in an interactive R session to explore the surface. ```{r buyske-surface, eval=FALSE} plot_tdt_surface3d( metric = "mssn", scenario = "misclassification", x = "prev", y = "misclass_rate", x_values = c(0.005, 0.01, 0.02, 0.04), y_values = seq(0, 0.05, length.out = 5), target_power = 0.80, alpha = 1e-5, pd = 0.50, R1 = sqrt(2.5), R2 = 2.5, delta_prime = 1, heter_rate = 0 ) ``` When run interactively, this surface extends discrete calculations into a sensitivity analysis; it is not an uncertainty distribution. # Published-example reproduction: Chen et al. (2009) Chen et al. studied TDT design under locus heterogeneity. Here `heter_rate` is the heterogeneous or unlinked trio fraction, equal to `1 - pi`. ```{r chen-reproduction} chen_no_heterogeneity <- tdt_mssn( target_power = 0.80, alpha = 1e-5, pd = 0.25, prev = 0.10, R1 = sqrt(1.5), R2 = 1.5, delta_prime = 1, misclass_rate = 0, heter_rate = 0, verbose = FALSE ) chen_half_heterogeneous <- tdt_mssn( target_power = 0.80, alpha = 1e-5, pd = 0.25, prev = 0.10, R1 = sqrt(1.5), R2 = 1.5, delta_prime = 1, misclass_rate = 0, heter_rate = 0.50, verbose = FALSE ) data.frame( heterogeneous_fraction = c(0, 0.50), published_required_trios = c("approximately 3430", "approximately 13721"), paweh_required_trios = c( chen_no_heterogeneity$N$heterogeneity, chen_half_heterogeneous$N$heterogeneity ) ) ``` This is a near-exact **Published-example reproduction**. Small differences reflect the paper's reported integer values and rounding. ## Extending Chen: dilution by heterogeneous trios ```{r chen-curve, fig.cap="Required affected-child trios as the heterogeneous or unlinked fraction increases under the Chen model."} plot_tdt_mssn( x_var = "heter_rate", x_values = seq(0, 0.50, by = 0.05), scenario = "heterogeneity", input_mode = "model_based", target_power = 0.80, alpha = 1e-5, pd = 0.25, prev = 0.10, R1 = sqrt(1.5), R2 = 1.5, delta_prime = 1, misclass_rate = 0, title = "Required trios vs locus heterogeneity" ) ``` Heterogeneous trios dilute locus-specific signal. A flatter curve would indicate a more robust design; a steep curve signals that modest departures from homogeneity require substantially more families. # Assumptions and limitations - The design is an affected-child trio TDT using an asymptotic chi-square approximation. - In the current parameterization, `delta_prime` lies in `[0, 1]` and scales maximum **positive** LD. Negative LD requires another allele-frequency- dependent normalization and is not implemented. - `misclass_rate` and `heter_rate` describe distinct reported scenarios. - Chapter 5 genotype misclassification is model-based, uses one scalar adjacent-genotype error rate, and cannot be combined with the other error mechanisms. - TDTae and sequencing/NGS genotype error are outside this implementation. - Model-free `ET` and `ENT` are treated as supplied expected counts. - MSSN is measured in trios. - Curves and surfaces vary fixed assumptions; they are not uncertainty bands. # References Buyske S, Yang G, Matise TC, Gordon D. When a case is not a case: effects of phenotype misclassification on power and sample size requirements for the transmission disequilibrium test with affected child trios. *Human Heredity*. 2009;67(4):287--292. . Chen C, Yang G, Buyske S, Matise T, Finch SJ, Gordon D. Transmission disequilibrium test power and sample size in the presence of locus heterogeneity. *Statistical Applications in Genetics and Molecular Biology*. 2009;8:Article 44. . Gordon D, Finch SJ, Kim W. *Heterogeneity in Statistical Genetics: How to Assess, Address, and Account for Mixtures in Association Studies*. Springer; 2020. . # Session information ```{r session-info} sessionInfo() ```