A new statistical method to analyze Morris Water Maze data using Dirichlet distribution

2.1 The Dirichlet distribution

The Dirichlet distribution is a multivariate generalization of beta distributions. It describes the distribution of K-dimensional vectors p for which the sum of all the coordinates is fixed, i.e. ${\displaystyle {\sum }_{k=1}^K{p}_k=1}$. It is parametrized by a K-dimensional vector α of positive reals α_k > 0, 1 ≤ k ≤ K, such that its probability density function is given by

where $B(\mathit{\alpha })={\displaystyle {\prod }_{k=1}^K\text{Γ}({\alpha }_k)}/\text{Γ}(s)$ is the multivariate beta function and $s={\displaystyle {\sum }_{k=1}^K{\alpha }_k}$ is the precision. The marginal distributions are beta distributions with parameters (α_k, s – α_k) with expectation values m_k = α_k/s, variance m_k (1 – m_k)/(s + 1) and covariance between coordinates p_i and p_j given by – m_im_j/(s + 1). Therefore, the higher the precision s, the less diffuse coordinates are around their means. The Dirichlet distribution is the most general distribution for fixed-sum variables, motivating its use to describe compositional or fractional data such as MWM data, where K = 4. The likelihood of a sample of N independent observations D = {p₁, . . . , p_N} is given by

2.2 Likelihood-ratio test based on the Dirichlet distribution

Description of the test To reflect memory abilities, we would like to test whether the fraction of time spent in the four quadrants significantly differs from a uniform distribution, thus showing preference for one or several quadrants. To do so, we propose a likelihood-ratio test based on the Dirichlet distribution to distinguish between the null hypothesis of uniformity H₀ : {∃α > 0, ∀ k, α_k = α} (implying that all means m_k are equal to 1/K but the precision is not constrained), and the general hypothesis H₁ where the α_k’s are unconstrained. The likelihood-ratio statistic reads

In order to fit the distribution parameters to their maximum likelihood values, we refer to the numerical schemes developed in 5 and we used the open-source Python module dirichlet implemented by Eric Suh² and run with Python 3.6. In particular 5, proposes a technique to alternatively fit the means m_k and precision s, faster than fitting directly the α_k’s. The maximum likelihood parameters under the null hypothesis are thus estimated by setting the means to their uniform value, m_k = 1/K, and fitting the precision s. We provide a slightly modified version of the dirichlet package, forked from that of Eric Suh, which is publicly available³. Under the null hypothesis, the likelihood-ratio statistic Λ asymptotically follows a χ²-distribution with K – 1 degrees of freedom.

Bartlett correction Biological samples are usually limited and for small samples the statistic’s distribution deviates from a ${\chi }_{K-1}^2,$ as can be seen in Figure 2. We propose an approximate Bartlett-type correction 6 for small samples, which amounts to rescale the likelihood-ratio statistic to match its asymptotic mean, which is K – 1 in this test. Such a correction has been shown to correctly reproduce the first three moments of the asymptotic χ²-distribution 6. In order to derive the scaling factor, we needed to compute the expected value of the statistic as a function of the number of samples N and the precision s. To do so, we drew random samples of N observations from uniform Dirichlet distributions with precision s, varying N between 2 and 100 and s between 1 and 1000 (with logarithmically-spaced values), and measured the mean of the statistic Λ⁴. We found that the mean value of the likelihood-ratio statistic Λ depends very little on s in the probed range, and that the difference to the asymptotic value of K – 1 is well-fitted with a power law in N, i.e. 〈Λ〉 – (K – 1) ~ a_KN^b_K (data not shown). We found the approximate values a_K = 5.9 and b_K = –1.4 for K = 4. We therefore propose to use a corrected statistic

Figure 2. Histogram of the likelihood-ratio statistic Λ for s = 10 for various sample sizes N.

Histograms of the likelihood-ratio statistic Λ are represented in different colors according to the N value. Means are represented by vertical lines of the same color. The ${\chi }_3^2$-distribution is represented in blue. For small samples, the distribution of the likelihood ratio slightly deviates from a ${\chi }_3^2$ distribution and the mean is significantly greater than the theoretical value of 3.

Validation of the test To validate this correction, we compared the distribution of the uncorrected statistic Λ and the corrected statistic $\tilde{\Lambda }$ from our simulated samples to a ${\chi }_{K-1}^2$-distribution with probability-probability plots. As shown in Figure 3, the uncorrected statistic yields p-values significantly different from the theoretical ones while the corrected p-values are in perfect agreement with the ${\chi }_3^2$-distribution. Therefore, this correction significantly improves the reliability of the test for small samples.

Figure 3. Probability-probability plots for uncorrected and corrected statistics from our simulated data.

The uncorrected and corrected statistics are compared with a ${\chi }_{K-1}^2$ distribution for several sample sizes N. Grey lines represent equality between Λ percentiles and ${\chi }_{K-1}^2$ percentiles. Blue lines correspond to the uncorrected statistics and green lines to the corrected one. There is a difference between the p-values from the uncorrected statistics and the theoretical ones that disappears after correction, especially for small sample size.

We also computed the number of false negatives on the simulated data to evaluate the rate of type 1 error (Figure 4). We found that using the p-value from the corrected statistic leads to a consistent rate of type 1 error (i.e. α =5%), independent on the number of samples N or the precision s. On the contrary, using the non-corrected statistic leads to more false negatives, especially when the number of samples is small.

Figure 4. Rate of false negatives depending on the sample size N for different values of the precision s for a p-value of 0.05.

The rate of false-negative (i.e. type 1 error) using the corrected version of the statistic is represented by the blue line whereas the rate of false-negative using the uncorrected version of the statistic is represented by the blue dotted line. We found that the correction results in a consistent rate of false negatives, independent of the sample size N and precision s in the range tested. The red line represents the rate of false-negative using a single-sample t-test on the target quadrant.

Comparison with the one-sample Student t-test In order to compare the type 2 error obtained with the Dirichlet distribution with the results obtained using a one-sample t-test on the target quadrant as often done in the literature, we simulated data from a non-uniform Dirichlet distribution with the parameters α = (40; 20; 20; 20). We found that the p-value from the corrected statistic of the Dirichlet distribution is mainly lower than the one obtained with a single t-test on the target quadrant (75% of the p-values are lower for s = 30). This means that for some cases where the target quadrant is preferred, using a one-sample t-test on the target quadrant would not detect this preference whereas the test based on Dirichlet distribution would detect the divergence from uniformity. Beyond improving the description and the interpretation of data from the MWM, the test we propose extracts more information from the same experiment as it is based on a larger dataset and then decreases the number of false-positives.

Post-hoc analysis Using the test based on Dirichlet distribution, we can determine whether the fraction of time spent in the quadrants is uniformly distributed. In the case of a divergence from uniformity, we would like to evaluate what are the quadrants responsible for this divergence as a post-hoc analysis.

This can be performed by comparing the marginal distributions of each quadrant, that are Beta distributions, to a theoretical Beta distribution with parameters α = 0.25s and β = 0.75s. The only simple way to compare one distribution with a theoretical one is to apply a single sample t-test that compares a normal distribution with a theoretical normal distribution. However, we noticed that in the range of inter-individual variability we have in this kind of study (given by the parameter s of the Dirichlet distribution, usually found between 20 and 50), the marginal distribution are fairly close to a normal distribution. Seeking for simplification, we advise to apply single sample t-tests for a post-hoc characterization of the preference for a quadrant in groups showing a divergence from uniformity.

2.3 Bayesian inference of the parameters of the Dirichlet distribution

Bayesian analysis can be used to infer constraints on the parameters α_i’s of the Dirichlet distribution used to model the data (and subsequently the means m_i’s). Specifically, we performed nested sampling of the parameter space of the Dirichlet distribution using the PyMC3 package with Python 3.6. For simplicity, we used the Jeffreys prior⁵ π(α) which does not depend on the model parametrization (e.g., sampling over α or (m, s)) and leave the discussion about this choice for future work. The output is a sample of vectors α distributed as the posterior given the data, i.e. 𝒫 (α | D) ∝ 𝒫 (D | α)π(α), which enables us to compare confidence regions for different groups and visualize the consistency with uniformity.

3.1 Application of the likelihood-ratio test based on the Dirichlet distribution

We compared the distribution obtained for each group to a uniform distribution in the probe test of the standard reference memory task and we found that the Dirichlet distribution obtained for wild-type mice was significantly different from a uniform distribution (p = 0.0021), whereas the one obtained for 3xTg mice did not differ from a uniform distribution (p = 0.26). We also propose a module, included in the Dirichlet package, to draw charts showing at the same time mean values with uncertainties⁶ and inter-individual variability according to Dirichlet distribution (Figure 5). This result shows that 3xTg AD mice display long term memory deficits, which is in accordance with previous observations 9.

Figure 5. Time spent in the four quadrants by wild-type and 3xTg AD mice.

In this plot, each column represents a sample and each color represents a quadrant. Mean values for the fraction of time spent in each quadrant is represented by a dotted line and the error bars on the means are approximated with the inverse Fisher information. For 3xTg mice the fraction of time spent in each quadrant is approximately similar leading to a uniform distribution (p = 0.26) whereas for wild-type mice the time spent in the target quadrant is significantly higher leading to a non-uniform distribution (p = 0.0021).

To better characterize long term memory in wild-type mice, we applied single sample t-tests on the four quadrants as a post-hoc analysis. We performed a one-tailed single sample t-test to assess whether the fraction of time spent in the target quadrant by wild-type mice is greater than the theoretical value 25%. Conversely, we performed a one-tailed single sample t-test to assess whether the fraction of time spent in the opposite quadrant by wildtype mice is lower than the theoretical value 25%. For adjacent quadrants we performed two-tailed single sample t-tests. We observed that the fraction of time spent in the target and opposite quadrants were respectively significantly higher (p = 0.025) and lower (p = 0.013) than 25%. The fraction of time spent in the adjacent quadrants did not differ from 25%.

Using this dataset with usual sample sizes for behavioral studies (N=7), we confirmed that our test is able to discriminate efficient and deficient memory abilities on real data.

3.2 Perspectives using Bayesian inference

We inferred constraints on the parameters of the Dirichlet distribution for wild-type and 3xTg mice. Figure 6 indicates compatibility of the data with uniformity for 3xTg mice and shows a clear preference for the target quadrant for wild-type mice suggesting memory deficits in 3xTg mice compared with wild-type.

Figure 6. Fraction of time spent in the four quadrants for wild-type and 3xTg AD mice in the case of Bayesian inference.

Corner plot representing constraints on the mean fractions of time m_i’s for the two data sets wild-type (blue) and 3xTg (green). The diagonal plots show the marginal distributions of m_i’s (with shaded 68% confidence interval) and off-diagonal plots show the two-dimensional distributions of pairs of these variables (inner and outer contours represent the 68% and 95% confidence levels). The black dashed lines represent the case of uniformity (25%) and the red lines correspond to equal time spent in both considered quadrant. Constraints on m₁ (leftmost column) indicate that wild-type mice favor the target quadrant.

Underlying data

A Python notebook with the code to reproduce the simulations and figures is available at https://github.com/xuod/dirichlet 12. In vivo dataset available from: https://github.com/xuod/dirichlet/tree/master/example 12. (Experimental procedures and data acquisition are detailed in 8).