Data-Story/Task7_下_假设检验2_多元数值向量检验/假设检验2-多元数值向量的检验.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 加载必要的包\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. 多元数值向量的假设检验"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "完成了对一元数值假设检验的学习，接下来我们简单学习一波常用的多元数值向量的假设检验。\n",
    "\n",
    "在进入学习之前，我们先思考一个问题：单变量问题和多变量问题的区别是什么——主要在于微观与宏观的区别。单变量问题研究的是具体某一个指标的性质，而多变量问题则是将多个指标作为一个整体，用它们去描述一个对象的性质。\n",
    "\n",
    "为了让大家有一个感性的认识，我们比较一下均值检验与均值向量检验的区别。如果我们想要了解两个城市在某个指标上的差异，如企业的资产负债率，那么我们肯定使用的是单变量均值检验；然而，如果问题变成了了解两个城市的企业发展状况的差异，我们还能使用单变量均值检验吗？肯定是不行的，因为企业发展状况需要用多个指标去共同描述。\n",
    "\n",
    "接下来可能有人会问：那我能不能分别对这些指标做单变量检验呢？且不说如果指标有很多个，依次做单变量检验的工作量会很大；最主要的问题是，我们难以通过这些分散的单变量检验回答诸如“两个城市的企业发展状况是否存在显著差异”的问题，假设在十个指标中，有三个指标有显著差异，七个指标没有显著差异，我们可以认为两个城市的企业发展状况不存在显著差异吗？显然没有一个准确的答案。\n",
    "\n",
    "因此，我们需要将这些指标作为一个整体去研究，才能解决上述问题。而在多元数值向量的假设检验中，最常用的假设检验便是均值向量的假设检验。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 均值向量的检验"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "本章节覆盖了多元数值数据中的**组间均值向量**假设检验的内容，我们主要学习如何通过样本数据对其所在多元总体的**均值向量**的性质进行检验，并展示如何在python中实现每个检验。均值向量检验的内容架构与均值检验的架构是一样的，主要内容有：\n",
    "\n",
    "1. 比较一组数据的总体均值向量与一个固定向量是否相等的检验\n",
    "2. 比较两组数据的总体均值向量之间是否相等的检验\n",
    "3. 比较两组以上的多组数据的总体均值向量之间是否相等的检验\n",
    "\n",
    "只不过，多元部分就没有解决非正态性问题的非参数检验了（或许是有的，但是我不知道……）。在多元领域，有一个检验和t检验非常类似——Hotelling T2检验，该检验在变量个数等于1时就退化成了t检验。在下面的均值向量检验中，我们会反复使用Hotelling T2检验。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1.1 单组样本均值向量假定的检验"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "单组样本均值向量检验与单组样本的均值检验类似，只不过假定的对象从一个数值变成了一个向量。对于这个问题，我们可以使用单样本hotelling T2检验。\n",
    "\n",
    "设$\\boldsymbol{x}_1,\\cdots ,\\boldsymbol{x}_n\\,$来一个$p$元正态分布的一个样本。对于单样本hotelling T2检验，我们做如下两个假设：\n",
    "$$\n",
    "H_0:\\boldsymbol{\\mu }=\\boldsymbol{\\mu }_0\\leftrightarrow H_1:\\boldsymbol{\\mu }\\ne \\boldsymbol{\\mu }_0\\,\n",
    "$$\n",
    "注意，这里的$\\boldsymbol{\\mu }$是加粗的，代表向量的意思。\n",
    "\n",
    "检验统计量为：\n",
    "$$\n",
    "Test\\,\\,statistics=\\frac{n-p}{p\\left( n-1 \\right)}T^2\\,\\,, T^2=n\\left( \\boldsymbol{\\bar{x}}-\\boldsymbol{\\mu }_0 \\right) '\\boldsymbol{S}^{-1}\\left( \\boldsymbol{\\bar{x}}-\\boldsymbol{\\mu }_0 \\right) \n",
    "$$\n",
    "其中，$\\boldsymbol{S}$样本的无偏协方差矩阵\n",
    "\n",
    "检验统计量服从分布\n",
    "$$\n",
    "Test\\,\\,statistics\\sim F\\left( p,n-p \\right) \n",
    "$$\n",
    "p值计算公式为\n",
    "$$\n",
    "pvalue=P\\left( F>Test\\,\\,statistics \\right) \n",
    "$$\n",
    "注意，尽管备择假设的符号是$\\ne$，但这里p值并不是双边p值，多元检验中p值形式的判断并不能沿用一元均值检验的判断方法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example.5** 我们想研究某地区农村2岁男婴的发育状况是否达到国家参考标准，于是抽取了6名男婴并测量了他们的身高$x_1$、胸围$x_2$、上半臀围$x_3$。已知这三个指标的国家参考标准均值为$\\boldsymbol{\\mu }_0=\\left( 90,58,16 \\right) '$，6名男婴的数据如下所示，问能否认为该地区农村男婴的发育状况达到国家标准？\n",
    "\n",
    "由于scipy包中并没有Hotelling T2的api，因此我们需要使用上述的检验统计量及其分布自己编写假设检验。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>x3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>78</td>\n",
       "      <td>60.6</td>\n",
       "      <td>16.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>76</td>\n",
       "      <td>58.1</td>\n",
       "      <td>12.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>92</td>\n",
       "      <td>63.2</td>\n",
       "      <td>14.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>81</td>\n",
       "      <td>59.0</td>\n",
       "      <td>14.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>81</td>\n",
       "      <td>60.8</td>\n",
       "      <td>15.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>84</td>\n",
       "      <td>59.5</td>\n",
       "      <td>14.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   x1    x2    x3\n",
       "0  78  60.6  16.5\n",
       "1  76  58.1  12.5\n",
       "2  92  63.2  14.5\n",
       "3  81  59.0  14.0\n",
       "4  81  60.8  15.5\n",
       "5  84  59.5  14.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data=pd.read_excel('./data/ex5.xlsx')\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 构造一个函数，帮助我们进行单样本Hotelling T2检验\n",
    "def multi_checkmean(data:pd.DataFrame,mean:np.ndarray,significance=0.05):\n",
    "    '''\n",
    "    data:待比较的矩阵，格式最好为pd.DataFrame()，以列为变量，以行为样本。\n",
    "    mean:假定的均值向量，格式为numpy.ndarray。\n",
    "    '''\n",
    "    from scipy.stats import f\n",
    "\n",
    "    # 计算检验统计量\n",
    "    S=np.cov(data.T)\n",
    "    S_inv=np.linalg.inv(S)\n",
    "    x_bar=np.mean(data).values.T\n",
    "    T2=len(data)*np.dot((x_bar-mean).T,np.dot(S_inv,x_bar-mean))\n",
    "    Test_statistics=T2*(len(data)-len(mean))/(len(mean)*(len(data)-1))\n",
    "\n",
    "    # 计算p值\n",
    "    pvalue=f.sf(Test_statistics,len(mean),len(data)-len(mean))\n",
    "\n",
    "    # 输出样本的均值向量\n",
    "    print('样本均值向量为：{}'.format(x_bar))\n",
    "    print('目标均值向量为：{}'.format(mean))\n",
    "\n",
    "    # 比较p值与显著性水平\n",
    "    if pvalue<significance:\n",
    "        print('在显著性水平{0:}下，样本均值向量不等于假定的均值向量。(p={1:.4f})'.format(significance,pvalue))\n",
    "    else:\n",
    "        print('在显著性水平{0:}下，样本均值向量等于假定的均值向量。(p={1:.4f})'.format(significance,pvalue))\n",
    "    return pvalue\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "样本均值向量为：[82.  60.2 14.5]\n",
      "目标均值向量为：[90 58 16]\n",
      "在显著性水平0.05下，样本均值向量不等于假定的均值向量。(p=0.0022)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.002155272045932366"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean=np.array([90,58,16]).T #均值向量\n",
    "multi_checkmean(data,mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可见，该地区农村男婴的发育状况与国家标准存在显著差异，且通过观察两个均值向量可以发现，农村男婴的发育状况低于国家标准。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1.2 两组样本之间的均值向量相等性检验"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "与单变量的双样本均值检验一样，双样本均值向量的检验同样分为组别间独立与组别成对两种情况，且判别两者的方法也是相同的，此处不再赘述。我们直接开始构建假设检验。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 组别间独立"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example.6** 我们想知道3周岁的男女婴儿的身体指标是否存在显著差异，于是各抽取了男女各15名婴儿，并测量了他们的身高$X_1$与体重$X_2$，问：在显著性水平0.05下能否认为两者存在显著差异？\n",
    "\n",
    "由于两个样本性别不同，因此我们可以认为两个样本相互独立。\n",
    "\n",
    "设两个样本$\\boldsymbol{x}_1,\\cdots ,\\boldsymbol{x}_{n1}\\,$和$\\boldsymbol{y}_1,\\cdots ,\\boldsymbol{y}_{n2}\\,$分别来自两个$p$元正态分布，我们做如下两个假设：\n",
    "$$\n",
    "H_0:\\boldsymbol{\\mu_1 }=\\boldsymbol{\\mu }_2\\leftrightarrow H_1:\\boldsymbol{\\mu_1 }\\ne \\boldsymbol{\\mu }_2\\,\n",
    "$$\n",
    "检验统计量为：\n",
    "$$\n",
    "Test\\,\\,statistics=\\frac{n_1+n_2-p-1}{p\\left( n_1+n_2-2 \\right)}T^2\\,\\,, T^2=\\frac{n_1n_2}{n_1+n_2}(\\overline{\\boldsymbol{x}}-\\overline{\\boldsymbol{y}})^{'}\\boldsymbol{S}_{p}^{-1}(\\overline{\\boldsymbol{x}}-\\overline{\\boldsymbol{y}}) , \\boldsymbol{S}_p=\\frac{\\left( n_1-1 \\right) \\boldsymbol{S}_1+\\left( n_2-1 \\right) \\boldsymbol{S}_2}{n_1+n_2-2}\n",
    "$$\n",
    "检验统计量服从分布：\n",
    "$$\n",
    "Test\\,\\,statistics\\sim F\\left( p,n_1+n_2-p-1 \\right) \n",
    "$$\n",
    "p值形式为：\n",
    "$$\n",
    "pvalue=P\\left( F>Test\\,\\,statistics \\right) \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>X1</th>\n",
       "      <th>X2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>54.0</td>\n",
       "      <td>3.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>53.0</td>\n",
       "      <td>2.25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>51.5</td>\n",
       "      <td>2.50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>51.0</td>\n",
       "      <td>3.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>51.0</td>\n",
       "      <td>3.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>77.0</td>\n",
       "      <td>7.50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>77.0</td>\n",
       "      <td>10.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>77.0</td>\n",
       "      <td>9.50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>74.0</td>\n",
       "      <td>9.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>73.0</td>\n",
       "      <td>7.50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>91.0</td>\n",
       "      <td>12.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>91.0</td>\n",
       "      <td>13.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>94.0</td>\n",
       "      <td>15.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>92.0</td>\n",
       "      <td>12.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>91.0</td>\n",
       "      <td>12.50</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      X1     X2\n",
       "15  54.0   3.00\n",
       "16  53.0   2.25\n",
       "17  51.5   2.50\n",
       "18  51.0   3.00\n",
       "19  51.0   3.00\n",
       "20  77.0   7.50\n",
       "21  77.0  10.00\n",
       "22  77.0   9.50\n",
       "23  74.0   9.00\n",
       "24  73.0   7.50\n",
       "25  91.0  12.00\n",
       "26  91.0  13.00\n",
       "27  94.0  15.00\n",
       "28  92.0  12.00\n",
       "29  91.0  12.50"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data=pd.read_excel('./data/ex6.xlsx')\n",
    "# 将两个独立样本分开成两组\n",
    "group_boy=data[data['gender']=='b'].drop('gender',axis=1) # 去除性别标签列，只保留数值列\n",
    "group_girl=data[data['gender']=='g'].drop('gender',axis=1)\n",
    "group_girl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def multi_unparied_data(group1:pd.DataFrame,group2:pd.DataFrame,significance=0.05):\n",
    "    # 计算检验统计量\n",
    "    n1=len(group1)\n",
    "    n2=len(group2)\n",
    "    p=np.shape(group1)[1] # 变量维度\n",
    "    mean1=np.mean(group1).values.T\n",
    "    mean2=np.mean(group2).values.T\n",
    "    S1=np.cov(group1.T)\n",
    "    S2=np.cov(group2.T)\n",
    "    Sp=((n1-1)*S1+(n2-1)*S2)/(n1+n2-2)\n",
    "    T2=n1*n2*(mean1-mean2).T@np.linalg.inv(Sp)@(mean1-mean2)/(n1+n2)\n",
    "    Test_statistics=(n1+n2-p-1)*T2/(p*(n1+n2-2))\n",
    "\n",
    "    # 计算p值\n",
    "    from scipy.stats import f\n",
    "\n",
    "    pvalue=f.sf(Test_statistics,p,n1+n2-p-1)\n",
    "\n",
    "    # 比较p值与显著性水平\n",
    "    if pvalue<significance:\n",
    "        print('在显著性水平{0:}下，两组样本所在总体的均值向量不相等。(p={1:.4f})'.format(significance,pvalue))\n",
    "    else:\n",
    "        print('在显著性水平{0:}下，两组样本所在总体的均值向量相等。(p={1:.4f})'.format(significancee,pvalue))\n",
    "    return pvalue\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "在显著性水平0.05下，两组样本所在总体的均值向量相等。(p=0.9396)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9396308874768952"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "multi_unparied_data(group_boy,group_girl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 成对组别"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "多元均值向量检验中，成对检验的原理与一元均值检验的成对检验原理是相似的，本质上做两个均值向量之差与零向量的之间的单样本均值向量检验。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "设$\\left( \\boldsymbol{x}_i,\\boldsymbol{y}_i \\right) \\,\\,, i=1,\\cdots ,n $为服从$p$元正态分布的成对数据，令他们相减得：\n",
    "$$\n",
    "\\boldsymbol{d}_i=\\boldsymbol{x}_i-\\boldsymbol{y}_i\n",
    "$$\n",
    "检验统计量为：\n",
    "$$\n",
    "Test\\,\\,statistics=\\frac{n-p}{p\\left( n-1 \\right)}T^2\\,\\,, T^2=n\\boldsymbol{\\bar{d}}'\\boldsymbol{S}_{\\boldsymbol{d}}^{-1}\\boldsymbol{\\bar{d}}\n",
    "$$\n",
    "检验统计量服从分布：\n",
    "$$\n",
    "Test\\,\\,statistics\\sim F\\left( p,n-p \\right) \n",
    "$$\n",
    "p值计算公式为\n",
    "$$\n",
    "pvalue=P\\left( F>Test\\,\\,statistics \\right) \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def multi_paried_data(group1:pd.DataFrame,group2:pd.DataFrame,significance=0.05):\n",
    "    # 计算检验统计量值\n",
    "    n=len(group1)\n",
    "    p=np.shape(group1)[1]\n",
    "    mean1=np.mean(group1).values.T\n",
    "    mean2=np.mean(group2).values.T\n",
    "    d=group1-group2\n",
    "    d_mean=mean1-mean2\n",
    "    Sd=np.cov(d.T)\n",
    "    Sd_inv=np.linalg.inv(Sd)\n",
    "    T2=n*d_mean.T@Sd_inv@d_mean\n",
    "    Test_statistics=(n-p)*T2/(p*(n-1))\n",
    "    \n",
    "    # 计算p值\n",
    "    from scipy.stats import f\n",
    "\n",
    "    pvalue=f.sf(Test_statistics,p,n-p)\n",
    "\n",
    "    # 比较p值与显著性水平\n",
    "    if pvalue<significance:\n",
    "        print('在显著性水平{0:}下，两组样本所在总体的均值向量不相等。(p={1:.4f})'.format(significance,pvalue))\n",
    "    else:\n",
    "        print('在显著性水平{0:}下，两组样本所在总体的均值向量相等。(p={1:.4f})'.format(significance,pvalue))\n",
    "    return pvalue\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1.3 多元方差分析-多组样本间的均值向量相等性比较"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "本次课程我们只学习单因素多元方差分析的应用。\n",
    "\n",
    "多元方差分析与一元方差分析研究的内容是相似的，只是前者研究的是多个总体均值向量的相等性，而后者是均值的相等性。多元方差分析的检验统计量被称为威尔克斯(Wilks)$\\varLambda $统计量，该统计量服从的分布为$\\varLambda $分布。这个分布的分位点值可通过查F分布表得到，不过两者的转换非常复杂，在这里就不做介绍了。\n",
    "\n",
    "万幸的是，statsmodels.multivariate.manova模块有多元方差分析(MANOVA)的api，既然有现成的api可以使用，那我们就直接拿来使用就好了！"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example.8** 某医生测量了16个正常人的早晨中3个小时各小时的低频心电频值谱（LF）与高频心电频值谱（HF），数据在下面给出，问：这16个人的两个指标表现在3次测量中有无显著的差别？\n",
    "\n",
    "我们简单对这个例子进行分析：由于问题考察的变量有两个：LF与HF，因此这是一个多元问题；又由于样本有三个：早晨的3次测量，因此这是一个多组均值比较问题。结合以上分析，该问题是一个（单因素）多元方差分析问题。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LF</th>\n",
       "      <th>HF</th>\n",
       "      <th>Time</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4.66</td>\n",
       "      <td>2.89</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.54</td>\n",
       "      <td>4.65</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5.91</td>\n",
       "      <td>4.53</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.95</td>\n",
       "      <td>3.31</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.51</td>\n",
       "      <td>3.78</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>4.22</td>\n",
       "      <td>2.61</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>4.61</td>\n",
       "      <td>3.10</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>5.08</td>\n",
       "      <td>4.38</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>3.71</td>\n",
       "      <td>1.76</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>3.63</td>\n",
       "      <td>3.17</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>4.49</td>\n",
       "      <td>4.08</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>5.70</td>\n",
       "      <td>4.78</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>4.96</td>\n",
       "      <td>3.39</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>5.83</td>\n",
       "      <td>4.02</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>5.22</td>\n",
       "      <td>5.08</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>4.15</td>\n",
       "      <td>2.39</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>4.29</td>\n",
       "      <td>3.03</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>4.69</td>\n",
       "      <td>4.77</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>5.28</td>\n",
       "      <td>4.41</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>5.05</td>\n",
       "      <td>3.28</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>4.94</td>\n",
       "      <td>3.56</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>4.54</td>\n",
       "      <td>3.28</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>4.26</td>\n",
       "      <td>3.11</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>5.56</td>\n",
       "      <td>5.36</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>3.96</td>\n",
       "      <td>2.47</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>3.64</td>\n",
       "      <td>3.19</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>4.86</td>\n",
       "      <td>4.12</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>5.72</td>\n",
       "      <td>5.44</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>5.14</td>\n",
       "      <td>3.88</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>5.64</td>\n",
       "      <td>4.06</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>5.03</td>\n",
       "      <td>4.99</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>4.15</td>\n",
       "      <td>2.08</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>4.77</td>\n",
       "      <td>3.57</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>4.58</td>\n",
       "      <td>3.04</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>5.37</td>\n",
       "      <td>4.79</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35</th>\n",
       "      <td>4.65</td>\n",
       "      <td>2.86</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>4.68</td>\n",
       "      <td>3.97</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>37</th>\n",
       "      <td>4.61</td>\n",
       "      <td>4.40</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>38</th>\n",
       "      <td>5.27</td>\n",
       "      <td>3.88</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39</th>\n",
       "      <td>5.55</td>\n",
       "      <td>5.00</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>4.16</td>\n",
       "      <td>2.70</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>3.30</td>\n",
       "      <td>3.10</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>4.64</td>\n",
       "      <td>3.87</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>5.54</td>\n",
       "      <td>4.89</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>5.21</td>\n",
       "      <td>3.88</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>5.26</td>\n",
       "      <td>3.84</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46</th>\n",
       "      <td>5.43</td>\n",
       "      <td>4.50</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>47</th>\n",
       "      <td>4.57</td>\n",
       "      <td>2.32</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      LF    HF  Time\n",
       "0   4.66  2.89     1\n",
       "1   4.54  4.65     1\n",
       "2   5.91  4.53     1\n",
       "3   4.95  3.31     1\n",
       "4   5.51  3.78     1\n",
       "5   4.22  2.61     1\n",
       "6   4.61  3.10     1\n",
       "7   5.08  4.38     1\n",
       "8   3.71  1.76     1\n",
       "9   3.63  3.17     1\n",
       "10  4.49  4.08     1\n",
       "11  5.70  4.78     1\n",
       "12  4.96  3.39     1\n",
       "13  5.83  4.02     1\n",
       "14  5.22  5.08     1\n",
       "15  4.15  2.39     1\n",
       "16  4.29  3.03     2\n",
       "17  4.69  4.77     2\n",
       "18  5.28  4.41     2\n",
       "19  5.05  3.28     2\n",
       "20  4.94  3.56     2\n",
       "21  4.54  3.28     2\n",
       "22  4.26  3.11     2\n",
       "23  5.56  5.36     2\n",
       "24  3.96  2.47     2\n",
       "25  3.64  3.19     2\n",
       "26  4.86  4.12     2\n",
       "27  5.72  5.44     2\n",
       "28  5.14  3.88     2\n",
       "29  5.64  4.06     2\n",
       "30  5.03  4.99     2\n",
       "31  4.15  2.08     2\n",
       "32  4.77  3.57     3\n",
       "33  4.58  3.04     3\n",
       "34  5.37  4.79     3\n",
       "35  4.65  2.86     3\n",
       "36  4.68  3.97     3\n",
       "37  4.61  4.40     3\n",
       "38  5.27  3.88     3\n",
       "39  5.55  5.00     3\n",
       "40  4.16  2.70     3\n",
       "41  3.30  3.10     3\n",
       "42  4.64  3.87     3\n",
       "43  5.54  4.89     3\n",
       "44  5.21  3.88     3\n",
       "45  5.26  3.84     3\n",
       "46  5.43  4.50     3\n",
       "47  4.57  2.32     3"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data=pd.read_excel('./data/ex8.xlsx')\n",
    "group_1=data[data['Time']==1].drop('Time',axis=1).reset_index(drop=True)\n",
    "group_2=data[data['Time']==2].drop('Time',axis=1).reset_index(drop=True)\n",
    "group_3=data[data['Time']==3].drop('Time',axis=1).reset_index(drop=True)\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                           Value Num DF Den DF   F Value   Pr > F\n",
      "Wilks' lambda            0.99114      2   45.0  0.201139  0.81853\n",
      "Pillai's trace           0.00886    2.0   45.0  0.201139  0.81853\n",
      "Hotelling-Lawley trace  0.008939      2   45.0  0.201139  0.81853\n",
      "Roy's greatest root     0.008939      2     45  0.201139  0.81853\n"
     ]
    }
   ],
   "source": [
    "from statsmodels.multivariate.manova import MANOVA\n",
    "\n",
    "model=MANOVA.from_formula('LF + HF ~ Time', data=data).mv_test()\n",
    "# 在''中填入公式，其中~左侧填入自变量名称，~右侧填入因素名称\n",
    "print(model.results['Time']['stat'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在上述报告表中，我们只需要看Wilks' lambda的p值即可。此处p值为0.81853，显然我们不能拒绝原假设，3次测量没有显著区别。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 从多元到一元——寻找显著的变量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "至此，一元数值与多元数值向量的常用假设检验已经学习完毕。虽然从研究对象来看，两者存在一定的区别，但是这并不意味着它们是相互独立、没有联系的。相反，在进行实际分析中，两者的功能是互补的。\n",
    "\n",
    "多元均值向量的检验可以告诉我们不同分析对象的某个宏观性质（如城市企业的发展状况）是否存在显著区别，如果我们拒绝了原假设，便可以说明它们确实存在显著区别。但是在实际分析中，仅仅得出“它们存在显著区别”的结论是远远不够的，我们会想知道，究竟是什么变量导致了这种显著区别呢？这个时候，我们就要是用一元的均值向量检验来回答这一问题。\n",
    "\n",
    "也就是说，在实际分析当中，往往都是始于“多元”，终于“一元”。为了让大家更好地理解，我们对以下例子进行一次实操演练。\n",
    "\n",
    "**Example.9** 为了研究3种销售方式对商品销售额的影响，我们选择了4种商品按照这3种销售方式进行销售。这四种商品的销售额分别为x1,x2,x3,x4。数据如下所示，请用多元/一元假设检验方法对该问题进行研究分析。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x1</th>\n",
       "      <th>x2</th>\n",
       "      <th>x3</th>\n",
       "      <th>x4</th>\n",
       "      <th>sale</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>125</td>\n",
       "      <td>60</td>\n",
       "      <td>338</td>\n",
       "      <td>210</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>119</td>\n",
       "      <td>80</td>\n",
       "      <td>233</td>\n",
       "      <td>330</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>63</td>\n",
       "      <td>51</td>\n",
       "      <td>260</td>\n",
       "      <td>203</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>65</td>\n",
       "      <td>51</td>\n",
       "      <td>429</td>\n",
       "      <td>150</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>130</td>\n",
       "      <td>65</td>\n",
       "      <td>403</td>\n",
       "      <td>205</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>69</td>\n",
       "      <td>45</td>\n",
       "      <td>350</td>\n",
       "      <td>190</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>46</td>\n",
       "      <td>60</td>\n",
       "      <td>585</td>\n",
       "      <td>200</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>146</td>\n",
       "      <td>66</td>\n",
       "      <td>273</td>\n",
       "      <td>250</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>87</td>\n",
       "      <td>54</td>\n",
       "      <td>585</td>\n",
       "      <td>240</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>110</td>\n",
       "      <td>77</td>\n",
       "      <td>507</td>\n",
       "      <td>270</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>107</td>\n",
       "      <td>60</td>\n",
       "      <td>364</td>\n",
       "      <td>200</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>130</td>\n",
       "      <td>61</td>\n",
       "      <td>391</td>\n",
       "      <td>200</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>80</td>\n",
       "      <td>45</td>\n",
       "      <td>429</td>\n",
       "      <td>270</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>60</td>\n",
       "      <td>50</td>\n",
       "      <td>442</td>\n",
       "      <td>190</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>81</td>\n",
       "      <td>54</td>\n",
       "      <td>260</td>\n",
       "      <td>280</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>135</td>\n",
       "      <td>87</td>\n",
       "      <td>507</td>\n",
       "      <td>260</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>57</td>\n",
       "      <td>48</td>\n",
       "      <td>400</td>\n",
       "      <td>285</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>75</td>\n",
       "      <td>52</td>\n",
       "      <td>520</td>\n",
       "      <td>260</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>76</td>\n",
       "      <td>65</td>\n",
       "      <td>403</td>\n",
       "      <td>250</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>55</td>\n",
       "      <td>42</td>\n",
       "      <td>411</td>\n",
       "      <td>170</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>66</td>\n",
       "      <td>54</td>\n",
       "      <td>455</td>\n",
       "      <td>310</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>82</td>\n",
       "      <td>45</td>\n",
       "      <td>403</td>\n",
       "      <td>210</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>65</td>\n",
       "      <td>65</td>\n",
       "      <td>312</td>\n",
       "      <td>280</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>40</td>\n",
       "      <td>51</td>\n",
       "      <td>477</td>\n",
       "      <td>280</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>67</td>\n",
       "      <td>54</td>\n",
       "      <td>481</td>\n",
       "      <td>293</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>38</td>\n",
       "      <td>50</td>\n",
       "      <td>468</td>\n",
       "      <td>210</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>42</td>\n",
       "      <td>45</td>\n",
       "      <td>351</td>\n",
       "      <td>190</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>113</td>\n",
       "      <td>40</td>\n",
       "      <td>390</td>\n",
       "      <td>310</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>80</td>\n",
       "      <td>55</td>\n",
       "      <td>520</td>\n",
       "      <td>200</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>76</td>\n",
       "      <td>60</td>\n",
       "      <td>507</td>\n",
       "      <td>189</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>94</td>\n",
       "      <td>33</td>\n",
       "      <td>260</td>\n",
       "      <td>280</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>60</td>\n",
       "      <td>51</td>\n",
       "      <td>429</td>\n",
       "      <td>190</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>55</td>\n",
       "      <td>40</td>\n",
       "      <td>390</td>\n",
       "      <td>295</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>65</td>\n",
       "      <td>48</td>\n",
       "      <td>481</td>\n",
       "      <td>177</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>69</td>\n",
       "      <td>48</td>\n",
       "      <td>442</td>\n",
       "      <td>225</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35</th>\n",
       "      <td>125</td>\n",
       "      <td>63</td>\n",
       "      <td>312</td>\n",
       "      <td>270</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>120</td>\n",
       "      <td>56</td>\n",
       "      <td>416</td>\n",
       "      <td>280</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>37</th>\n",
       "      <td>70</td>\n",
       "      <td>45</td>\n",
       "      <td>468</td>\n",
       "      <td>370</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>38</th>\n",
       "      <td>62</td>\n",
       "      <td>66</td>\n",
       "      <td>416</td>\n",
       "      <td>224</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39</th>\n",
       "      <td>69</td>\n",
       "      <td>60</td>\n",
       "      <td>377</td>\n",
       "      <td>280</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>65</td>\n",
       "      <td>33</td>\n",
       "      <td>480</td>\n",
       "      <td>260</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>100</td>\n",
       "      <td>34</td>\n",
       "      <td>468</td>\n",
       "      <td>295</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>65</td>\n",
       "      <td>63</td>\n",
       "      <td>416</td>\n",
       "      <td>265</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>117</td>\n",
       "      <td>48</td>\n",
       "      <td>468</td>\n",
       "      <td>250</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>114</td>\n",
       "      <td>63</td>\n",
       "      <td>395</td>\n",
       "      <td>380</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>55</td>\n",
       "      <td>30</td>\n",
       "      <td>546</td>\n",
       "      <td>235</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46</th>\n",
       "      <td>64</td>\n",
       "      <td>51</td>\n",
       "      <td>507</td>\n",
       "      <td>320</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>47</th>\n",
       "      <td>110</td>\n",
       "      <td>90</td>\n",
       "      <td>442</td>\n",
       "      <td>225</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>48</th>\n",
       "      <td>60</td>\n",
       "      <td>62</td>\n",
       "      <td>440</td>\n",
       "      <td>248</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>49</th>\n",
       "      <td>110</td>\n",
       "      <td>69</td>\n",
       "      <td>377</td>\n",
       "      <td>260</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50</th>\n",
       "      <td>88</td>\n",
       "      <td>78</td>\n",
       "      <td>299</td>\n",
       "      <td>360</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>51</th>\n",
       "      <td>73</td>\n",
       "      <td>63</td>\n",
       "      <td>390</td>\n",
       "      <td>320</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52</th>\n",
       "      <td>114</td>\n",
       "      <td>55</td>\n",
       "      <td>494</td>\n",
       "      <td>240</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>53</th>\n",
       "      <td>103</td>\n",
       "      <td>54</td>\n",
       "      <td>416</td>\n",
       "      <td>310</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54</th>\n",
       "      <td>100</td>\n",
       "      <td>33</td>\n",
       "      <td>273</td>\n",
       "      <td>312</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55</th>\n",
       "      <td>140</td>\n",
       "      <td>61</td>\n",
       "      <td>312</td>\n",
       "      <td>345</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>56</th>\n",
       "      <td>80</td>\n",
       "      <td>36</td>\n",
       "      <td>286</td>\n",
       "      <td>250</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>57</th>\n",
       "      <td>135</td>\n",
       "      <td>54</td>\n",
       "      <td>468</td>\n",
       "      <td>345</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>58</th>\n",
       "      <td>130</td>\n",
       "      <td>69</td>\n",
       "      <td>325</td>\n",
       "      <td>360</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59</th>\n",
       "      <td>60</td>\n",
       "      <td>57</td>\n",
       "      <td>273</td>\n",
       "      <td>260</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     x1  x2   x3   x4  sale\n",
       "0   125  60  338  210     1\n",
       "1   119  80  233  330     1\n",
       "2    63  51  260  203     1\n",
       "3    65  51  429  150     1\n",
       "4   130  65  403  205     1\n",
       "5    69  45  350  190     1\n",
       "6    46  60  585  200     1\n",
       "7   146  66  273  250     1\n",
       "8    87  54  585  240     1\n",
       "9   110  77  507  270     1\n",
       "10  107  60  364  200     1\n",
       "11  130  61  391  200     1\n",
       "12   80  45  429  270     1\n",
       "13   60  50  442  190     1\n",
       "14   81  54  260  280     1\n",
       "15  135  87  507  260     1\n",
       "16   57  48  400  285     1\n",
       "17   75  52  520  260     1\n",
       "18   76  65  403  250     1\n",
       "19   55  42  411  170     1\n",
       "20   66  54  455  310     2\n",
       "21   82  45  403  210     2\n",
       "22   65  65  312  280     2\n",
       "23   40  51  477  280     2\n",
       "24   67  54  481  293     2\n",
       "25   38  50  468  210     2\n",
       "26   42  45  351  190     2\n",
       "27  113  40  390  310     2\n",
       "28   80  55  520  200     2\n",
       "29   76  60  507  189     2\n",
       "30   94  33  260  280     2\n",
       "31   60  51  429  190     2\n",
       "32   55  40  390  295     2\n",
       "33   65  48  481  177     2\n",
       "34   69  48  442  225     2\n",
       "35  125  63  312  270     2\n",
       "36  120  56  416  280     2\n",
       "37   70  45  468  370     2\n",
       "38   62  66  416  224     2\n",
       "39   69  60  377  280     2\n",
       "40   65  33  480  260     3\n",
       "41  100  34  468  295     3\n",
       "42   65  63  416  265     3\n",
       "43  117  48  468  250     3\n",
       "44  114  63  395  380     3\n",
       "45   55  30  546  235     3\n",
       "46   64  51  507  320     3\n",
       "47  110  90  442  225     3\n",
       "48   60  62  440  248     3\n",
       "49  110  69  377  260     3\n",
       "50   88  78  299  360     3\n",
       "51   73  63  390  320     3\n",
       "52  114  55  494  240     3\n",
       "53  103  54  416  310     3\n",
       "54  100  33  273  312     3\n",
       "55  140  61  312  345     3\n",
       "56   80  36  286  250     3\n",
       "57  135  54  468  345     3\n",
       "58  130  69  325  360     3\n",
       "59   60  57  273  260     3"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data=pd.read_excel('./data/ex9.xlsx')\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们先简单分析一波问题：该例中有三个样本，因此我们会使用方差分析进行研究；销售额方面，本例用了四个商品的销售额去衡量，也就是说变量有4个，因此在分析的第一步我们可以使用多元方差分析，分析三种销售方式在“宏观”上是否会给商品销售额带来影响。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                           Value Num DF Den DF   F Value    Pr > F\n",
      "Wilks' lambda           0.749149      4   55.0  4.604152  0.002807\n",
      "Pillai's trace          0.250851    4.0   55.0  4.604152  0.002807\n",
      "Hotelling-Lawley trace  0.334847      4   55.0  4.604152  0.002807\n",
      "Roy's greatest root     0.334847      4     55  4.604152  0.002807\n"
     ]
    }
   ],
   "source": [
    "# 先做多元方差分析\n",
    "model=MANOVA.from_formula('x1 + x2 + x3 + x4 ~ sale', data=data).mv_test()\n",
    "print(model.results['sale']['stat'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "p值为0.002807，显然三种销售方式的均值向量不全部相等。分析到这一步肯定是不够的，我们可以进一步思考：这三种销售方式的显著差异究竟是由哪些商品引起的呢？这个问题可以转化成——这三种销售方式会对哪种商品的销售产生显著的差异，为了研究这个问题，我们可以分别对变量x1,x2,x3,x4做一元方差分析。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F_onewayResult(statistic=3.3766410140335, pvalue=0.04112546209315577)\n",
      "F_onewayResult(statistic=1.6154481009614068, pvalue=0.2077698930751813)\n",
      "F_onewayResult(statistic=0.16557150166792145, pvalue=0.8478153535976287)\n",
      "F_onewayResult(statistic=8.008491041468082, pvalue=0.0008607020541776114)\n"
     ]
    }
   ],
   "source": [
    "# 分别检验\n",
    "\n",
    "## 先分组\n",
    "group_1=data[data['sale']==1].drop('sale',axis=1).reset_index(drop=True)\n",
    "group_2=data[data['sale']==2].drop('sale',axis=1).reset_index(drop=True)\n",
    "group_3=data[data['sale']==3].drop('sale',axis=1).reset_index(drop=True)\n",
    "\n",
    "## 分别做单变量方差检验\n",
    "print(stats.f_oneway(group_1.x1.values,group_2.x1.values,group_3.x1.values))\n",
    "print(stats.f_oneway(group_1.x2.values,group_2.x2.values,group_3.x2.values))\n",
    "print(stats.f_oneway(group_1.x3.values,group_2.x3.values,group_3.x3.values))\n",
    "print(stats.f_oneway(group_1.x4.values,group_2.x4.values,group_3.x4.values))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "x2~x3的p值很大，因此我们可以认为这三种销售方式对于商品2/3而言不会产生显著的影响；x1、x4的p值小于0.05，且x4的p值尤其小，我们猜测，可能是商品4销售额均值的差异引起了均值向量的差异。对于该问题，我们可以剔除x4，对x1~x3再做一次多元方差分析。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                           Value Num DF Den DF   F Value    Pr > F\n",
      "Wilks' lambda           0.976477      3   56.0  0.449666  0.718524\n",
      "Pillai's trace          0.023523    3.0   56.0  0.449666  0.718524\n",
      "Hotelling-Lawley trace  0.024089      3   56.0  0.449666  0.718524\n",
      "Roy's greatest root     0.024089      3     56  0.449666  0.718524\n"
     ]
    }
   ],
   "source": [
    "# 去掉变量x4，再做一次多元方差分析\n",
    "model=MANOVA.from_formula('x1 + x2 + x3 ~ sale', data=data).mv_test()\n",
    "print(model.results['sale']['stat'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "去除了x4后，多元方差分析检验不显著，因此说明对于商品1/2/3而言，三种销售方式的总体均值向量之间没有显著差异；此外，尽管对x1单独进行方差分析是显著的，可是x1~x3一起做多元方差分析时却不显著，这说明商品1对三种销售方式的差异无明显影响，即“你的差异影响不了大伙”~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 作业"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为研究东、中、西部各省市规模以上的企业发展状况，我们收集了各城市企业的主要经济指标，包括：总资产贡献率、资产负债率、流动资产周转次数、工业成本费用利润率、产品销售率。我们用变量“类别”定义了各类城市，其中1为东部城市；2为中部城市；3为西部城市。数据文件为homework2.xlsx。假设显著性水平为$\\alpha$，问：\n",
    "\n",
    "1. 对三个类别的城市进行均值向量间的两两比较，查看结果\n",
    "\n",
    "2. 对三个类别的城市同时进行均值向量间的比较，查看结果\n",
    "\n",
    "3. 承接问题2，你认为哪些变量导致了三个类别城市均值向量的差异？说出你的理由。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.7 ('base')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "7758e92e9a61d7a3490898707f7eeb937c85e9d1e8d4e877cc6c187218f226d5"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}