models can account for both additive and dominant effects of a genetic variant on expression levels. In the parametric ANOVA model, suppose the test samples are divided into $\ell+1$ groups by the SNP $j$. The mean expression level for gene $i$ in each group $m$ can be written as $$\begin{aligned}
  \overline{X_i^{(m,j)}} = \frac1{n^{(m,j)}} \sum_{\{l\colon G_{jl}=m\}} X_{il},\end{aligned}$$ where $n^{(m,j)}$ is the number of samples in genotype group $m$ for SNP $j$.

Again assuming that the expression data is standardised, the F-test statistic for testing gene $i$ against SNP $j$ can be written as $$\begin{aligned}
  F_i^{(j)} = \frac{n-\ell-1}{\ell} \frac{SS_i^{(j)}}{n-SS_i^{(j)}},\end{aligned}$$ where $SS_i^{(j)}$ is the sum of squares between groups, $$\begin{aligned}
  SS_i^{(j)} = \sum_{m=0}^\ell n^{(m,j)}\overline{X_i^{(m,j)}}^2.\end{aligned}$$

Let us define the $n\times s$ indicator matrix ${\mathbf{I}}^{(m)}$ for genotype group $m$, i.e. ${\mathbf{I}}_{lj}^{(m)} = 1$ if $G_{jl}=m$ and $0$ otherwise. Then $$\begin{aligned}
  \sum_{\{l\colon G_{jl}=m\}} X_{il} = \left({\mathbf{X}}{\mathbf{I}}^{(m)}\right)_{ij}.\end{aligned}$$ Hence, for each pair of expression level $X_i$ and SNP $G_j$, the sum of squares matrix $SS_i^{(j)}$ can be computed via $\ell-1$ matrix multiplications [^1].

In the non-parametric ANOVA model, the expression data matrix is converted to a matrix ${\mathbf{T}}$ of data ranks, independently over each row. In the absence of ties, the Kruskal-Wallis test statistic is given by $$\begin{aligned}
  