Coverage for /usr/lib/python3/dist-packages/scipy/stats/_page_trend_test.py: 24%
100 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from itertools import permutations
2import numpy as np
3import math
4from ._continuous_distns import norm
5import scipy.stats
6from dataclasses import dataclass
9@dataclass
10class PageTrendTestResult:
11 statistic: float
12 pvalue: float
13 method: str
16def page_trend_test(data, ranked=False, predicted_ranks=None, method='auto'):
17 r"""
18 Perform Page's Test, a measure of trend in observations between treatments.
20 Page's Test (also known as Page's :math:`L` test) is useful when:
22 * there are :math:`n \geq 3` treatments,
23 * :math:`m \geq 2` subjects are observed for each treatment, and
24 * the observations are hypothesized to have a particular order.
26 Specifically, the test considers the null hypothesis that
28 .. math::
30 m_1 = m_2 = m_3 \cdots = m_n,
32 where :math:`m_j` is the mean of the observed quantity under treatment
33 :math:`j`, against the alternative hypothesis that
35 .. math::
37 m_1 \leq m_2 \leq m_3 \leq \cdots \leq m_n,
39 where at least one inequality is strict.
41 As noted by [4]_, Page's :math:`L` test has greater statistical power than
42 the Friedman test against the alternative that there is a difference in
43 trend, as Friedman's test only considers a difference in the means of the
44 observations without considering their order. Whereas Spearman :math:`\rho`
45 considers the correlation between the ranked observations of two variables
46 (e.g. the airspeed velocity of a swallow vs. the weight of the coconut it
47 carries), Page's :math:`L` is concerned with a trend in an observation
48 (e.g. the airspeed velocity of a swallow) across several distinct
49 treatments (e.g. carrying each of five coconuts of different weight) even
50 as the observation is repeated with multiple subjects (e.g. one European
51 swallow and one African swallow).
53 Parameters
54 ----------
55 data : array-like
56 A :math:`m \times n` array; the element in row :math:`i` and
57 column :math:`j` is the observation corresponding with subject
58 :math:`i` and treatment :math:`j`. By default, the columns are
59 assumed to be arranged in order of increasing predicted mean.
61 ranked : boolean, optional
62 By default, `data` is assumed to be observations rather than ranks;
63 it will be ranked with `scipy.stats.rankdata` along ``axis=1``. If
64 `data` is provided in the form of ranks, pass argument ``True``.
66 predicted_ranks : array-like, optional
67 The predicted ranks of the column means. If not specified,
68 the columns are assumed to be arranged in order of increasing
69 predicted mean, so the default `predicted_ranks` are
70 :math:`[1, 2, \dots, n-1, n]`.
72 method : {'auto', 'asymptotic', 'exact'}, optional
73 Selects the method used to calculate the *p*-value. The following
74 options are available.
76 * 'auto': selects between 'exact' and 'asymptotic' to
77 achieve reasonably accurate results in reasonable time (default)
78 * 'asymptotic': compares the standardized test statistic against
79 the normal distribution
80 * 'exact': computes the exact *p*-value by comparing the observed
81 :math:`L` statistic against those realized by all possible
82 permutations of ranks (under the null hypothesis that each
83 permutation is equally likely)
85 Returns
86 -------
87 res : PageTrendTestResult
88 An object containing attributes:
90 statistic : float
91 Page's :math:`L` test statistic.
92 pvalue : float
93 The associated *p*-value
94 method : {'asymptotic', 'exact'}
95 The method used to compute the *p*-value
97 See Also
98 --------
99 rankdata, friedmanchisquare, spearmanr
101 Notes
102 -----
103 As noted in [1]_, "the :math:`n` 'treatments' could just as well represent
104 :math:`n` objects or events or performances or persons or trials ranked."
105 Similarly, the :math:`m` 'subjects' could equally stand for :math:`m`
106 "groupings by ability or some other control variable, or judges doing
107 the ranking, or random replications of some other sort."
109 The procedure for calculating the :math:`L` statistic, adapted from
110 [1]_, is:
112 1. "Predetermine with careful logic the appropriate hypotheses
113 concerning the predicted ording of the experimental results.
114 If no reasonable basis for ordering any treatments is known, the
115 :math:`L` test is not appropriate."
116 2. "As in other experiments, determine at what level of confidence
117 you will reject the null hypothesis that there is no agreement of
118 experimental results with the monotonic hypothesis."
119 3. "Cast the experimental material into a two-way table of :math:`n`
120 columns (treatments, objects ranked, conditions) and :math:`m`
121 rows (subjects, replication groups, levels of control variables)."
122 4. "When experimental observations are recorded, rank them across each
123 row", e.g. ``ranks = scipy.stats.rankdata(data, axis=1)``.
124 5. "Add the ranks in each column", e.g.
125 ``colsums = np.sum(ranks, axis=0)``.
126 6. "Multiply each sum of ranks by the predicted rank for that same
127 column", e.g. ``products = predicted_ranks * colsums``.
128 7. "Sum all such products", e.g. ``L = products.sum()``.
130 [1]_ continues by suggesting use of the standardized statistic
132 .. math::
134 \chi_L^2 = \frac{\left[12L-3mn(n+1)^2\right]^2}{mn^2(n^2-1)(n+1)}
136 "which is distributed approximately as chi-square with 1 degree of
137 freedom. The ordinary use of :math:`\chi^2` tables would be
138 equivalent to a two-sided test of agreement. If a one-sided test
139 is desired, *as will almost always be the case*, the probability
140 discovered in the chi-square table should be *halved*."
142 However, this standardized statistic does not distinguish between the
143 observed values being well correlated with the predicted ranks and being
144 _anti_-correlated with the predicted ranks. Instead, we follow [2]_
145 and calculate the standardized statistic
147 .. math::
149 \Lambda = \frac{L - E_0}{\sqrt{V_0}},
151 where :math:`E_0 = \frac{1}{4} mn(n+1)^2` and
152 :math:`V_0 = \frac{1}{144} mn^2(n+1)(n^2-1)`, "which is asymptotically
153 normal under the null hypothesis".
155 The *p*-value for ``method='exact'`` is generated by comparing the observed
156 value of :math:`L` against the :math:`L` values generated for all
157 :math:`(n!)^m` possible permutations of ranks. The calculation is performed
158 using the recursive method of [5].
160 The *p*-values are not adjusted for the possibility of ties. When
161 ties are present, the reported ``'exact'`` *p*-values may be somewhat
162 larger (i.e. more conservative) than the true *p*-value [2]_. The
163 ``'asymptotic'``` *p*-values, however, tend to be smaller (i.e. less
164 conservative) than the ``'exact'`` *p*-values.
166 References
167 ----------
168 .. [1] Ellis Batten Page, "Ordered hypotheses for multiple treatments:
169 a significant test for linear ranks", *Journal of the American
170 Statistical Association* 58(301), p. 216--230, 1963.
172 .. [2] Markus Neuhauser, *Nonparametric Statistical Test: A computational
173 approach*, CRC Press, p. 150--152, 2012.
175 .. [3] Statext LLC, "Page's L Trend Test - Easy Statistics", *Statext -
176 Statistics Study*, https://www.statext.com/practice/PageTrendTest03.php,
177 Accessed July 12, 2020.
179 .. [4] "Page's Trend Test", *Wikipedia*, WikimediaFoundation,
180 https://en.wikipedia.org/wiki/Page%27s_trend_test,
181 Accessed July 12, 2020.
183 .. [5] Robert E. Odeh, "The exact distribution of Page's L-statistic in
184 the two-way layout", *Communications in Statistics - Simulation and
185 Computation*, 6(1), p. 49--61, 1977.
187 Examples
188 --------
189 We use the example from [3]_: 10 students are asked to rate three
190 teaching methods - tutorial, lecture, and seminar - on a scale of 1-5,
191 with 1 being the lowest and 5 being the highest. We have decided that
192 a confidence level of 99% is required to reject the null hypothesis in
193 favor of our alternative: that the seminar will have the highest ratings
194 and the tutorial will have the lowest. Initially, the data have been
195 tabulated with each row representing an individual student's ratings of
196 the three methods in the following order: tutorial, lecture, seminar.
198 >>> table = [[3, 4, 3],
199 ... [2, 2, 4],
200 ... [3, 3, 5],
201 ... [1, 3, 2],
202 ... [2, 3, 2],
203 ... [2, 4, 5],
204 ... [1, 2, 4],
205 ... [3, 4, 4],
206 ... [2, 4, 5],
207 ... [1, 3, 4]]
209 Because the tutorial is hypothesized to have the lowest ratings, the
210 column corresponding with tutorial rankings should be first; the seminar
211 is hypothesized to have the highest ratings, so its column should be last.
212 Since the columns are already arranged in this order of increasing
213 predicted mean, we can pass the table directly into `page_trend_test`.
215 >>> from scipy.stats import page_trend_test
216 >>> res = page_trend_test(table)
217 >>> res
218 PageTrendTestResult(statistic=133.5, pvalue=0.0018191161948127822,
219 method='exact')
221 This *p*-value indicates that there is a 0.1819% chance that
222 the :math:`L` statistic would reach such an extreme value under the null
223 hypothesis. Because 0.1819% is less than 1%, we have evidence to reject
224 the null hypothesis in favor of our alternative at a 99% confidence level.
226 The value of the :math:`L` statistic is 133.5. To check this manually,
227 we rank the data such that high scores correspond with high ranks, settling
228 ties with an average rank:
230 >>> from scipy.stats import rankdata
231 >>> ranks = rankdata(table, axis=1)
232 >>> ranks
233 array([[1.5, 3. , 1.5],
234 [1.5, 1.5, 3. ],
235 [1.5, 1.5, 3. ],
236 [1. , 3. , 2. ],
237 [1.5, 3. , 1.5],
238 [1. , 2. , 3. ],
239 [1. , 2. , 3. ],
240 [1. , 2.5, 2.5],
241 [1. , 2. , 3. ],
242 [1. , 2. , 3. ]])
244 We add the ranks within each column, multiply the sums by the
245 predicted ranks, and sum the products.
247 >>> import numpy as np
248 >>> m, n = ranks.shape
249 >>> predicted_ranks = np.arange(1, n+1)
250 >>> L = (predicted_ranks * np.sum(ranks, axis=0)).sum()
251 >>> res.statistic == L
252 True
254 As presented in [3]_, the asymptotic approximation of the *p*-value is the
255 survival function of the normal distribution evaluated at the standardized
256 test statistic:
258 >>> from scipy.stats import norm
259 >>> E0 = (m*n*(n+1)**2)/4
260 >>> V0 = (m*n**2*(n+1)*(n**2-1))/144
261 >>> Lambda = (L-E0)/np.sqrt(V0)
262 >>> p = norm.sf(Lambda)
263 >>> p
264 0.0012693433690751756
266 This does not precisely match the *p*-value reported by `page_trend_test`
267 above. The asymptotic distribution is not very accurate, nor conservative,
268 for :math:`m \leq 12` and :math:`n \leq 8`, so `page_trend_test` chose to
269 use ``method='exact'`` based on the dimensions of the table and the
270 recommendations in Page's original paper [1]_. To override
271 `page_trend_test`'s choice, provide the `method` argument.
273 >>> res = page_trend_test(table, method="asymptotic")
274 >>> res
275 PageTrendTestResult(statistic=133.5, pvalue=0.0012693433690751756,
276 method='asymptotic')
278 If the data are already ranked, we can pass in the ``ranks`` instead of
279 the ``table`` to save computation time.
281 >>> res = page_trend_test(ranks, # ranks of data
282 ... ranked=True, # data is already ranked
283 ... )
284 >>> res
285 PageTrendTestResult(statistic=133.5, pvalue=0.0018191161948127822,
286 method='exact')
288 Suppose the raw data had been tabulated in an order different from the
289 order of predicted means, say lecture, seminar, tutorial.
291 >>> table = np.asarray(table)[:, [1, 2, 0]]
293 Since the arrangement of this table is not consistent with the assumed
294 ordering, we can either rearrange the table or provide the
295 `predicted_ranks`. Remembering that the lecture is predicted
296 to have the middle rank, the seminar the highest, and tutorial the lowest,
297 we pass:
299 >>> res = page_trend_test(table, # data as originally tabulated
300 ... predicted_ranks=[2, 3, 1], # our predicted order
301 ... )
302 >>> res
303 PageTrendTestResult(statistic=133.5, pvalue=0.0018191161948127822,
304 method='exact')
306 """
308 # Possible values of the method parameter and the corresponding function
309 # used to evaluate the p value
310 methods = {"asymptotic": _l_p_asymptotic,
311 "exact": _l_p_exact,
312 "auto": None}
313 if method not in methods:
314 raise ValueError(f"`method` must be in {set(methods)}")
316 ranks = np.array(data, copy=False)
317 if ranks.ndim != 2: # TODO: relax this to accept 3d arrays?
318 raise ValueError("`data` must be a 2d array.")
320 m, n = ranks.shape
321 if m < 2 or n < 3:
322 raise ValueError("Page's L is only appropriate for data with two "
323 "or more rows and three or more columns.")
325 if np.any(np.isnan(data)):
326 raise ValueError("`data` contains NaNs, which cannot be ranked "
327 "meaningfully")
329 # ensure NumPy array and rank the data if it's not already ranked
330 if ranked:
331 # Only a basic check on whether data is ranked. Checking that the data
332 # is properly ranked could take as much time as ranking it.
333 if not (ranks.min() >= 1 and ranks.max() <= ranks.shape[1]):
334 raise ValueError("`data` is not properly ranked. Rank the data or "
335 "pass `ranked=False`.")
336 else:
337 ranks = scipy.stats.rankdata(data, axis=-1)
339 # generate predicted ranks if not provided, ensure valid NumPy array
340 if predicted_ranks is None:
341 predicted_ranks = np.arange(1, n+1)
342 else:
343 predicted_ranks = np.array(predicted_ranks, copy=False)
344 if (predicted_ranks.ndim < 1 or
345 (set(predicted_ranks) != set(range(1, n+1)) or
346 len(predicted_ranks) != n)):
347 raise ValueError(f"`predicted_ranks` must include each integer "
348 f"from 1 to {n} (the number of columns in "
349 f"`data`) exactly once.")
351 if type(ranked) is not bool:
352 raise TypeError("`ranked` must be boolean.")
354 # Calculate the L statistic
355 L = _l_vectorized(ranks, predicted_ranks)
357 # Calculate the p-value
358 if method == "auto":
359 method = _choose_method(ranks)
360 p_fun = methods[method] # get the function corresponding with the method
361 p = p_fun(L, m, n)
363 page_result = PageTrendTestResult(statistic=L, pvalue=p, method=method)
364 return page_result
367def _choose_method(ranks):
368 '''Choose method for computing p-value automatically'''
369 m, n = ranks.shape
370 if n > 8 or (m > 12 and n > 3) or m > 20: # as in [1], [4]
371 method = "asymptotic"
372 else:
373 method = "exact"
374 return method
377def _l_vectorized(ranks, predicted_ranks):
378 '''Calculate's Page's L statistic for each page of a 3d array'''
379 colsums = ranks.sum(axis=-2, keepdims=True)
380 products = predicted_ranks * colsums
381 Ls = products.sum(axis=-1)
382 Ls = Ls[0] if Ls.size == 1 else Ls.ravel()
383 return Ls
386def _l_p_asymptotic(L, m, n):
387 '''Calculate the p-value of Page's L from the asymptotic distribution'''
388 # Using [1] as a reference, the asymptotic p-value would be calculated as:
389 # chi_L = (12*L - 3*m*n*(n+1)**2)**2/(m*n**2*(n**2-1)*(n+1))
390 # p = chi2.sf(chi_L, df=1, loc=0, scale=1)/2
391 # but this is insentive to the direction of the hypothesized ranking
393 # See [2] page 151
394 E0 = (m*n*(n+1)**2)/4
395 V0 = (m*n**2*(n+1)*(n**2-1))/144
396 Lambda = (L-E0)/np.sqrt(V0)
397 # This is a one-sided "greater" test - calculate the probability that the
398 # L statistic under H0 would be greater than the observed L statistic
399 p = norm.sf(Lambda)
400 return p
403def _l_p_exact(L, m, n):
404 '''Calculate the p-value of Page's L exactly'''
405 # [1] uses m, n; [5] uses n, k.
406 # Switch convention here because exact calculation code references [5].
407 L, n, k = int(L), int(m), int(n)
408 _pagel_state.set_k(k)
409 return _pagel_state.sf(L, n)
412class _PageL:
413 '''Maintains state between `page_trend_test` executions'''
415 def __init__(self):
416 '''Lightweight initialization'''
417 self.all_pmfs = {}
419 def set_k(self, k):
420 '''Calculate lower and upper limits of L for single row'''
421 self.k = k
422 # See [5] top of page 52
423 self.a, self.b = (k*(k+1)*(k+2))//6, (k*(k+1)*(2*k+1))//6
425 def sf(self, l, n):
426 '''Survival function of Page's L statistic'''
427 ps = [self.pmf(l, n) for l in range(l, n*self.b + 1)]
428 return np.sum(ps)
430 def p_l_k_1(self):
431 '''Relative frequency of each L value over all possible single rows'''
433 # See [5] Equation (6)
434 ranks = range(1, self.k+1)
435 # generate all possible rows of length k
436 rank_perms = np.array(list(permutations(ranks)))
437 # compute Page's L for all possible rows
438 Ls = (ranks*rank_perms).sum(axis=1)
439 # count occurrences of each L value
440 counts = np.histogram(Ls, np.arange(self.a-0.5, self.b+1.5))[0]
441 # factorial(k) is number of possible permutations
442 return counts/math.factorial(self.k)
444 def pmf(self, l, n):
445 '''Recursive function to evaluate p(l, k, n); see [5] Equation 1'''
447 if n not in self.all_pmfs:
448 self.all_pmfs[n] = {}
449 if self.k not in self.all_pmfs[n]:
450 self.all_pmfs[n][self.k] = {}
452 # Cache results to avoid repeating calculation. Initially this was
453 # written with lru_cache, but this seems faster? Also, we could add
454 # an option to save this for future lookup.
455 if l in self.all_pmfs[n][self.k]:
456 return self.all_pmfs[n][self.k][l]
458 if n == 1:
459 ps = self.p_l_k_1() # [5] Equation 6
460 ls = range(self.a, self.b+1)
461 # not fast, but we'll only be here once
462 self.all_pmfs[n][self.k] = {l: p for l, p in zip(ls, ps)}
463 return self.all_pmfs[n][self.k][l]
465 p = 0
466 low = max(l-(n-1)*self.b, self.a) # [5] Equation 2
467 high = min(l-(n-1)*self.a, self.b)
469 # [5] Equation 1
470 for t in range(low, high+1):
471 p1 = self.pmf(l-t, n-1)
472 p2 = self.pmf(t, 1)
473 p += p1*p2
474 self.all_pmfs[n][self.k][l] = p
475 return p
478# Maintain state for faster repeat calls to page_trend_test w/ method='exact'
479_pagel_state = _PageL()