S_{ij} = \frac{12}{n(n+1)} \sum_{m=0}^\ell n^{(m,j)}\,\overline{T_i^{(m,j)}}^2 - 3(n+1),\end{aligned}$$ where $\overline{T_i^{(m,j)}}$ is the average expression rank of gene $i$ in genotype group $m$ of SNP $j$, defined as $$\begin{aligned}
  \overline{T_i^{(m,j)}}&=\frac1{n^{(m,j)}} \sum_{\{l\colon G_{jl}=m\}}T_{il},\end{aligned}$$ which can be similarly obtained from the $\ell-1$ matrix multiplications.

There is as yet no consensus about which statistical model is most appropriate for eQTL detection. Non-parametric methods were introduced in the earliest eQTL studies [@brem2002; @schadt2008] and have remained popular, as they are robust against variations in the underlying genetic model and trait distribution. More recently, the linear model implemented in matrix-eQTL has been used in a number of large-scale studies [@lappalainen2013transcriptome; @ardlie2015genotype]. A comparison on a dataset of 102 human whole blood samples showed that the parametric ANOVA method was highly sensitive to the presence of outlying gene expression values and SNPs with singleton genotype group. Linear models reported the highest number of eQTL associations after empirical False Discovery Rate (FDR) correction, with an expected bias towards additive linear associations. The Kruskal-Wallis test was most robust against data outliers and heterogeneous genotype group sizes and detected a higher proportion of non-linear associations, but was more conservative for calling additive linear associations than linear models [@qi2014].

In summary, when 