*This program should be used only for inference on case-control designs with independent (i.e., unrelated) subjects and biallelic variants.* For all tests, any variant with only a single genotype in the sample or completely missing genotypes is eliminated from consideration. The number of polymorphic variants, defined as those having 2 or more different genotypes in the sample, is reported by the software with the output of each test for convenience.

The CA max/sum tests base inferences on the permutation null distributions of the maximum/sum over single-variant Cochran-Armitage chi-square statistics in the group of variants under consideration. Kinnamon et al. (2012) found that these tests had power comparable to or greater than the CMC and WSS tests in the presence of extensive neutral variation and missing genotypes. Therefore,*we recommend the use of these tests over the others implemented in RVASSOC.*

Letting $Q_t$ denote the value of the summary statistic in permutation $t$ and $Q_{obs}$ denote the observed value in the sample, the two-sided p-value for the CA max/sum test is estimated from $m$ permutations by:

$$\hat{p}=\frac{\#(Q_t \ge Q_{obs})+1}{m+1}$$ This p-value follows "`P:`" in the output. The program also estimates the permutation null distribution probability mass function at $Q_{obs}$ by: $$\hat{p}_{tied}=\frac{\#(Q_t=Q_{obs})+1}{m+1}$$ This probability mass follows "`P(tied):`" in the output. High values of $\hat{p}_{tied}$ are indicative of a highly discrete permutation null distribution.

Because two floating-point representations of the same number may differ due to rounding and truncation errors, all comparisons of test statistics in the above calculations are performed using Knuth's algorithm, as implemented in the GSL function `gsl_fcmp`. Briefly, Knuth's algorithm returns: $$\begin{array}{rcl}

Q_t \lt Q_{obs} & \mathrm{if} & Q_t-Q_{obs} \lt -2^{k} \epsilon \\

Q_t \gt Q_{obs} & \mathrm{if} & Q_t-Q_{obs} \gt 2^{k} \epsilon \\

Q_t=Q_{obs} & \mathrm{if} & |Q_t-Q_{obs}| \le 2^{k} \epsilon

\end{array}$$ In these formulas, $k$ is the larger of the base-2 exponents of $Q_t$ and $Q_{obs}$, and $\epsilon$ is the floating-point equality relative epsilon set by the flag `-fceps`.

The CA max/sum tests use all available genotype data at every variant. The average per-individual percentage of available genotypes is reported by the software. In cases in which there are as many or more individuals with missing genotypes at a variant in the entire sample than there are individuals in either the case or control groups, the CA statistic may not be calculable at that variant in the observed sample or in one or more permutations. If this occurs, RVASSOC issues a warning message to STDERR, and the analysis should be rerun with the offending variants excluded. *Make sure to keep the STDERR so that important error and warning messages are not lost!*

The CMC test (Li and Leal 2008) collapses variants having overall sample MAF (i.e., MAF estimate from both cases and controls) ≤ the value specified by `-cmccolmaf` into a single presence/absence indicator in the analysis. An approximate $F$ test based Hotelling's $T^2$ is used for inference on the vector comprising this indicator and minor allele counts at the common variants. This vector has a number of elements equaling the number of common variants plus 1.

One issue not considered in the original paper is that linkage disequilibrium (LD) among common variants can reduce the number of linearly independent elements in this random vector, which leads to a singular covariance matrix for which a standard inverse does not exist. However, calculating Hotelling’s $T^2$ statistic with any generalized inverse is equivalent to calculating the statistic with a standard inverse on a full-rank subset of linearly independent common variants (Kinnamon et al. 2012). Goodnight (1979) provides an algorithm for automatically calculating a $g2$ generalized inverse and the dimension of the full-rank subset without any prior knowledge of the full-rank subset. The algorithm involves applying the G2SWEEP operator once to each of the columns of the covariance matrix in succession. This operator zeros the rows and columns corresponding to common variants that are linearly dependent on the previous common variant minor allele counts and/or the rare variant indicator. Linear dependence is determined numerically using a relative tolerance specified by the flag `-cmcswptol`. The effective number of linearly independent vector elements, $v$, is then automatically obtained by subtracting the number of columns that are zeroed from the total number of columns in the covariance matrix. The p-value is then calculated using the approximate $F(v,N-v-1)$ distribution of the appropriately scaled Hotelling’s $T^2$ statistic calculated using the $g2$ generalized inverse of the covariance matrix, where $N$ is the total sample size.

Only individuals with complete genotype data at common variants can be used in calculating Hotelling’s $T^2$. Provided genotype data at all common variants are complete, individuals with missing genotype data at rare variants can be used if at least one minor allele is present for a variant with a non-missing genotype because the coding of the rare variant indicator would be 1 regardless of the other variant genotypes. However, if no minor allele is present at any variant with non-missing genotypes, the coding of the rare variant is ambiguous because it would depend on the values of the unobserved genotypes. Therefore, such individuals also have to be excluded from calculating Hotelling’s $T^2$. The software outputs the number of individuals with complete data who were used in the CMC. If there are not at least 2 cases, 2 controls, and 1 denominator degree of freedom (ddf), the CMC test cannot be performed and the results will read `F(-99,-99): -99 P: -99`. *With ddf≤4, results are reported, but they may not be reliable based on our experiments and the fact that the variance of the $F$ distribution exists only with ddf>4.*

In cases where there are as many common variants as polymorphic variants (i.e., no rare variants below MAF cutoff), the CMC yields the same result as Hotelling's $T^2$ test on all variants because the column of zeros for the rare variant indicator is eliminated as linearly dependent by the G2SWEEP operator.

The implementation of the WSS test follows the description in Madsen and Browning (2009) with four modifications. First, midranks are used to break ties in genetic scores when calculating the rank-sum statistic, $W$. Second, a two-sided p-value is used. Third, the two-sided p-value is estimated directly from the permutation distribution of $W$. Letting $W_t$ denote the value in permutation $t$, $W_{obs}$ denote the observed value in the sample, and $\overline{W}$ denote the mean of $W$ over all $m$ permutations, the two-sided p-value is estimated as: $$\hat{p}^{WSS}=\frac{\#(|W_t - \overline{W}|\ge|W_{obs} - \overline{W}|)+1}{m+1}$$ This p-value follows "`P:`" in the output. The program also estimates the permutation null distribution probability mass function at $|W_{obs} - \overline{W}|$ by: $$\hat{p}^{WSS}_{tied}=\frac{\#(|W_t - \overline{W}|=|W_{obs} - \overline{W}|)+1}{m+1}$$ This probability mass follows "`P(tied):`" in the output. High values of $\hat{p}^{WSS}_{tied}$ are indicative of a highly discrete permutation null distribution. Again, all comparisons in the above calculations are performed using the `gsl_fcmp` implementation of Knuth's algorithm. Finally, missing single-variant genotypes, which were not considered in the original paper, are not used in estimation of the MAF in controls and assigned values of 0 so as not to contribute to the WSS in an individual. This procedure is equivalent to calculating the genetic score over only nonmissing genotypes in each individual.