{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "4e45deeb",
            "cell_type": "markdown",
            "source": "<div class='alert alert-warning'>\n\nSciPy's interactive examples with Jupyterlite are experimental and may not always work as expected. Execution of cells containing imports may result in large downloads (up to 60MB of content for the first import from SciPy). Load times when importing from SciPy may take roughly 10-20 seconds. If you notice any problems, feel free to open an [issue](https://github.com/scipy/scipy/issues/new/choose).\n\n</div>",
            "metadata": {}
        },
        {
            "id": "5cb48cab",
            "cell_type": "markdown",
            "source": "We use the example from [3]: 10 students are asked to rate three\nteaching methods - tutorial, lecture, and seminar - on a scale of 1-5,\nwith 1 being the lowest and 5 being the highest. We have decided that\na confidence level of 99% is required to reject the null hypothesis in\nfavor of our alternative: that the seminar will have the highest ratings\nand the tutorial will have the lowest. Initially, the data have been\ntabulated with each row representing an individual student's ratings of\nthe three methods in the following order: tutorial, lecture, seminar.\n",
            "metadata": {}
        },
        {
            "id": "7652190c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "table = [[3, 4, 3],\n         [2, 2, 4],\n         [3, 3, 5],\n         [1, 3, 2],\n         [2, 3, 2],\n         [2, 4, 5],\n         [1, 2, 4],\n         [3, 4, 4],\n         [2, 4, 5],\n         [1, 3, 4]]",
            "outputs": []
        },
        {
            "id": "a2cf07bc",
            "cell_type": "markdown",
            "source": "Because the tutorial is hypothesized to have the lowest ratings, the\ncolumn corresponding with tutorial rankings should be first; the seminar\nis hypothesized to have the highest ratings, so its column should be last.\nSince the columns are already arranged in this order of increasing\npredicted mean, we can pass the table directly into `page_trend_test`.\n",
            "metadata": {}
        },
        {
            "id": "faac7bc4",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats import page_trend_test\nres = page_trend_test(table)\nres",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "PageTrendTestResult(statistic=133.5, pvalue=0.0018191161948127822,\n                    method='exact')"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "ac86de2b",
            "cell_type": "markdown",
            "source": "This *p*-value indicates that there is a 0.1819% chance that\nthe $L$ statistic would reach such an extreme value under the null\nhypothesis. Because 0.1819% is less than 1%, we have evidence to reject\nthe null hypothesis in favor of our alternative at a 99% confidence level.\n\nThe value of the $L$ statistic is 133.5. To check this manually,\nwe rank the data such that high scores correspond with high ranks, settling\nties with an average rank:\n",
            "metadata": {}
        },
        {
            "id": "f9ad6068",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats import rankdata\nranks = rankdata(table, axis=1)\nranks",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[1.5, 3. , 1.5],\n       [1.5, 1.5, 3. ],\n       [1.5, 1.5, 3. ],\n       [1. , 3. , 2. ],\n       [1.5, 3. , 1.5],\n       [1. , 2. , 3. ],\n       [1. , 2. , 3. ],\n       [1. , 2.5, 2.5],\n       [1. , 2. , 3. ],\n       [1. , 2. , 3. ]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "b157ee9f",
            "cell_type": "markdown",
            "source": "We add the ranks within each column, multiply the sums by the\npredicted ranks, and sum the products.\n",
            "metadata": {}
        },
        {
            "id": "782f24e8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nm, n = ranks.shape\npredicted_ranks = np.arange(1, n+1)\nL = (predicted_ranks * np.sum(ranks, axis=0)).sum()\nres.statistic == L",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "a2a40c6f",
            "cell_type": "markdown",
            "source": "As presented in [3], the asymptotic approximation of the *p*-value is the\nsurvival function of the normal distribution evaluated at the standardized\ntest statistic:\n",
            "metadata": {}
        },
        {
            "id": "330d3172",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.stats import norm\nE0 = (m*n*(n+1)**2)/4\nV0 = (m*n**2*(n+1)*(n**2-1))/144\nLambda = (L-E0)/np.sqrt(V0)\np = norm.sf(Lambda)\np",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.0012693433690751756"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "b08ee56c",
            "cell_type": "markdown",
            "source": "This does not precisely match the *p*-value reported by `page_trend_test`\nabove. The asymptotic distribution is not very accurate, nor conservative,\nfor $m \\leq 12$ and $n \\leq 8$, so `page_trend_test` chose to\nuse ``method='exact'`` based on the dimensions of the table and the\nrecommendations in Page's original paper [1]. To override\n`page_trend_test`'s choice, provide the `method` argument.\n",
            "metadata": {}
        },
        {
            "id": "119126e5",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res = page_trend_test(table, method=\"asymptotic\")\nres",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "PageTrendTestResult(statistic=133.5, pvalue=0.0012693433690751756,\n                    method='asymptotic')"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "d0c39aa6",
            "cell_type": "markdown",
            "source": "If the data are already ranked, we can pass in the ``ranks`` instead of\nthe ``table`` to save computation time.\n",
            "metadata": {}
        },
        {
            "id": "bf466b69",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res = page_trend_test(ranks,             # ranks of data\n                      ranked=True,       # data is already ranked\n                      )\nres",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "PageTrendTestResult(statistic=133.5, pvalue=0.0018191161948127822,\n                    method='exact')"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "a4452763",
            "cell_type": "markdown",
            "source": "Suppose the raw data had been tabulated in an order different from the\norder of predicted means, say lecture, seminar, tutorial.\n",
            "metadata": {}
        },
        {
            "id": "a24eee7f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "table = np.asarray(table)[:, [1, 2, 0]]",
            "outputs": []
        },
        {
            "id": "9cd7d8be",
            "cell_type": "markdown",
            "source": "Since the arrangement of this table is not consistent with the assumed\nordering, we can either rearrange the table or provide the\n`predicted_ranks`. Remembering that the lecture is predicted\nto have the middle rank, the seminar the highest, and tutorial the lowest,\nwe pass:\n",
            "metadata": {}
        },
        {
            "id": "8c86d756",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res = page_trend_test(table,             # data as originally tabulated\n                      predicted_ranks=[2, 3, 1],  # our predicted order\n                      )\nres",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "PageTrendTestResult(statistic=133.5, pvalue=0.0018191161948127822,\n                    method='exact')"
                    },
                    "execution_count": null
                }
            ]
        }
    ]
}