Statistics_Based_Data_Analysis9.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7ebfbf0d-60d2-45fa-83fd-eb19fc75a1d6",
   "metadata": {},
   "source": [
    "# 통계적 검정 방법\n",
    "## 분산분석(ANOVA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e1e36a4f-bfe7-4065-ae2b-43c4f4d7d8cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)   # FutureWarning 제거\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import math\n",
    "from scipy import stats\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline\n",
    "\n",
    "font_family = 'Malgun Gothic'\n",
    "\n",
    "sns.set(font=font_family, rc={'axes.unicode_minus' : False})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0911fe11-c87f-4bf0-b6b7-5a02db4166cc",
   "metadata": {},
   "source": [
    "### 일원분산분석 ( One-way ANOVA )\n",
    "- 3개 이상의 범주를 갖는 하나의 실험 요인(독립변수)에 대한 하나의 종속변수의 평균 차이를 검정하는 분석\n",
    "- stats.f_oneway() : scipy f 통계량\n",
    "- statsmodels.stats.anova() : statsmodels anova\n",
    "---\n",
    "- 증명하고자 하는 가설\n",
    "  - 4개의 각기 다른 신입사원 교육훈련 기법의 효과성을 평가하고자 한다. 새로 입사한 32명의 시입사원에게 4개 교육기관에서 다른 학습 방법을 적용시켰다. 한 달 간의 훈련기간이 끝난 후 시험을 봤는데 그 점수는 아래와 같다. 4개의 교육훈련 기법 간 차이가 있는가? 만약 있다면 어떻게 다른가?\n",
    "  - -> 4개 이상의 대응표본을 비교해야 하므로 일원배치 분산분석 기법을 사용해야 한다.\n",
    "  \n",
    "  - 귀무가설 : 4개의 교육훈련 기법 간의 차이가 없다.\n",
    "  - 대립가설 : 4개의 교육훈련 기법 간의 차이가 있다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "349b3b63-0913-4e52-bbd2-a68a4add74a0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a 평균 : 79.0\n",
      "b 평균 : 65.375\n",
      "c 평균 : 67.125\n",
      "d 평균 : 69.5\n"
     ]
    }
   ],
   "source": [
    "# 데이터셋\n",
    "a = [66, 74, 82, 75, 73, 97, 87, 78]\n",
    "b = [72, 51, 59, 62, 74, 64, 78, 63]\n",
    "c = [61, 60, 57, 60, 81, 77, 70, 71]\n",
    "d = [63, 61, 76, 84, 58, 65, 69, 80]\n",
    "print(f'a 평균 : {np.mean(a)}')\n",
    "print(f'b 평균 : {np.mean(b)}')\n",
    "print(f'c 평균 : {np.mean(c)}')\n",
    "print(f'd 평균 : {np.mean(d)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "85becef6-e867-429a-b1bb-c2e334424fea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD9CAYAAAChtfywAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQVklEQVR4nO3dTWxU9f7H8U/PDIXSKVTGjhlqbiIJZoAEJIRUBUYTWkWxArGAxFYe04XKQhbiwhsXikajG25NE3RBjNUEgyjVkJT4NDwUEzcqT0VTIiWxLam32AKlpXP+C/JvLkKdmfZMz3zL+7UrM5x++UHeHX7nzJk813VdAQDMcvweAAAwOoQcAIwj5ABgHCEHAOMIOQAYR8gBwDhCDgDGBf36xv/97yUlk7l9CXs4HFJXV6/fY4wLrKW3WE9vWVhPx8nTHXcU3vIx30KeTLo5H3JJJma0grX0FuvpLcvrydYKABhHyAHAOEIOAMYRcgAwjpADgHGEHACMI+QAYJxv15H7KR4v0+nTpzw5Viw2S4nED54cCwBG4rYMebrhjUSmqLPzryxPAwCjw9YKABhHyAHAOEIOAMYRcgAwjpADgHGEHACMI+QAYBwhBwDjCDkAGEfIAcA4Qg4AxhFyADCOkAOAcYQcAIwj5ABgHCEHAOMIOQAYR8gBwDhCDgDGEXIAMI6QA4BxhBwAjCPkAGAcIQcA4wg5ABgXTOdJu3bt0nfffadr165pzZo1qqqq0sGDB1VfX69gMKjy8nLV1tZme1YAwC2kDHlzc7OOHz+uhoYGDQwMaP369Vq8eLHq6urU0NCggoIC1dTUKB6PKxaLjcXMAID/kXJr5cSJE1q0aJHy8vKUn5+veDyu5uZmxeNxhUIhBQIBVVZWKpFIjMW8AIC/SfmKfM6cOfrggw+0atUqJZNJHTt2TEuWLFFpaenQc6LRqFpbWzP6xuFwKPNpfVBSUuT3COMGa+kt1tNbltczZcgfeOAB/frrr9qwYYPuvPNOzZw5U8lkUoFAYOg5juPIcTI7b9rV1atk0s184jF24UKP3yOMCyUlRaylh1hPb1lYT8fJG/YFcFr1ffbZZ/Xxxx9r586dunTpkgYGBtTR0TH0eHt7u6LRqDfTAgAykjLkyWRSyWRSkvTzzz+rra1NTz/9tJqamtTX16fBwUE1NjaqoqIi68MCAG6Wcmvl8uXL2rhxoxzHUSgU0ptvvqlwOKxNmzapurparuuqqqrqhj1zAMDYyXNd15eNagt75JHIFHV2/uX3GOOChT1IS1hPb1lYz1HvkQMAchchBwDjCDkAGEfIAcA4Qg4AxhFyADCOkAOAcYQcAIwj5ABgHCEHAOMIOQAYR8gBwDhCDgDGEXIAMI6QA4BxhBwAjCPkAGAcIQcA4wg5ABhHyAHAuKDfAwBAtsTjZTp9+pQnx4rFZimR+MGTY3mNkAMYt9INbyQyRZ2df2V5muxhawUAjBt3r8jvvfdf6u7u9ux4kciUUR+juLhYZ86c82AaALjZuAt5d3e3Z/9FKikp0oULPaM+jhc/DABgOGytAIBxhBwAjCPkAGAcIQcA4wg5ABhHyAHAOEIOAMaNu+vIAeu8uj9ILt8bBN4i5ECOSSe+1u8NAm+lFfJ9+/Zp79696u/v18MPP6znnntOBw8eVH19vYLBoMrLy1VbW5vtWQEAt5Ay5D09Pdq1a5f279+vYDCompoaLV26VHV1dWpoaFBBQYFqamoUj8cVi8XGYmYAwP9IebIzPz9f+fn56u3t1ZUrV9Tf36+WlhbF43GFQiEFAgFVVlYqkUiMxbwAgL9J+Yp84sSJ2rBhg5YtW6YJEyZo3bp1+vPPP1VaWjr0nGg0qtbW1oy+cTgcynzaNJWUFOXcsbycySrWwFusp7csr2fKkLe0tOjLL7/Ut99+q2AwqJdeekmzZ89WYWHh0HMcx5HjZHYlY1dXr5JJN/OJ0+DFHQsl7+5+KHk3k1VeriWuYz29levr6Th5w74ATlnfQ4cO6ZFHHtHkyZOVn5+vlStXas+ePero6Bh6Tnt7u6LRqHcTAwDSljLkM2fO1NGjR+W61189JxIJrVmzRk1NTerr69Pg4KAaGxtVUVGR9WEBADdLubXy0EMP6eTJk1q3bp0kac6cOdq8ebMikYiqq6vluq6qqqpu2DMHAIydPPf/X2qPsWztkXv5RgkvPyHodn/zBnvk3uLflLcsrOeo9sgBALmNkAOAcdxrBaPm1U2eJG70BIwEIceocZMnwF9srQCAcYQcAIwj5ABgHHvkAEy6995/qbu727PjRSJTRn2M4uJinTlzzoNpMkPIAZjU3d2dk2/+88O4C/mBrUvVs2uDJ8fy6n2IB7Yu9ehIAHCzcRfyx/7zdc79lH4sMkWd//ZgIAC4BU52AoBxhBwAjCPkAGAcIQcA4wg5ABhHyAHAOEIOAMYRcgAwjpADgHGEHACMG3dv0Zf8u3HNcIqLi/0eAcA4Nu5C7uXHifHxZAAsYGsFAIwj5ABgHCEHAOPG3R45vJOLH6Ul+fdxWqPFeiJbCDmGlYsfpSXl3lVJ6WI9kS1srQCAcYQcAIwj5ABgHCEHAOMIOQAYR8gBwLiUlx9++umn2r9//9DXv//+uzZv3qzp06ervr5ewWBQ5eXlqq2tzeqgAIBbSxny1atXa/Xq1ZKk/v5+rVu3Tk899ZSeeeYZNTQ0qKCgQDU1NYrH44rFYlkfGABwo4y2Vj777DMtW7ZMR44cUTweVygUUiAQUGVlpRKJRLZmBAD8g7Tf2em6rj755BN99NFH2rt3r0pLS4cei0ajam1tzegbh8OhjJ7vl5KSIr9H8JWXf/5cPdZYytU1YD29O5Yfa5l2yJubmzV//nwVFRWpv79fhYWFQ485jiPHyey8aVdXr5JJN6Pf4wev3gZtlVd/fi/fUi7Z/XthPb2Vi+uZrbV0nLxhXwCnXd+9e/dq+fLlkqRIJKKOjo6hx9rb2xWNRkc5JgBgJNIKueu6+umnn7RgwQJJ0pIlS9TU1KS+vj4NDg6qsbFRFRUVWR0UAHBraW2t/Pbbb7r77ruHtk/C4bA2bdqk6upqua6rqqqqG/bMAQBjJ62Qz5w5U7t3777h11auXKmVK1dmYSQAQCZ4ZycAGEfIAcA4Qg4AxhFyADCOkAOAcYQcAIwj5ABgXNr3WgGAXHJg61L17NrgybG8ujvKga1LPTpSZgg5AJMe+8/X6uz8y5NjeXXTrMciU9T5bw8GyhBbKwBgHCEHAONuy62VeLxMp0+fSuu5kciUf3w8FpulROIHL8YCgBG5LUOebni9vnk/AGTDbRlywA+5eJWF5N+VFvAOIQfGSC5eZSH5d6UFvMPJTgAwjpADgHGEHACMI+QAYBwhBwDjCDkAGEfIAcA4riPHsHgDC2ADIceweAMLYANbKwBgHCEHAOPYWgFgVqrbTI+14uJiX74vIQdgklfnb6TrPxC8PN5YY2sFAIwj5ABgHCEHAOPYI8c/yrWTSZJ/J5SAXEXIMSxOJgE2sLUCAMal9Yq8t7dXr776qs6ePatgMKidO3fql19+UX19vYLBoMrLy1VbW5vtWQEAt5BWyHfs2KF4PK53331X0vWw19XVqaGhQQUFBaqpqVE8HlcsFsvqsACAm6XcWunq6tK5c+e0YsWKoV87cuSI4vG4QqGQAoGAKisrlUgksjooAODWUr4iP3PmjKZPn65t27apvb1d8+fPVzgcVmlp6dBzotGoWltbM/rG4XAo82l9UFJS5PcI4wZrmZtXAd1xxx383cj2v8+UIe/q6tKpU6f04Ycfqri4WK+88oqmTp2qQCAw9BzHceQ4mZ037erqVTLpZj7xGPLy1qvQbb+WuXwV0O3+dyPl/ho4Tt6wL4BT1jccDmv+/PmaNm2aHMfRo48+qj179qijo2PoOe3t7YpGo95NDABIW8qQ33fffTp58qR6e3slSUePHtXzzz+vpqYm9fX1aXBwUI2NjaqoqMj6sACAm6XcWikoKNCLL76oLVu2KBAIaPbs2VqxYoXy8vJUXV0t13VVVVV1w545AGDs5Lmu68tGNXvktxfe2ekt1tNbFtZzVHvkAIDcRsgBwDhumgVg3IrHy3T69Km0npvqGv9YbJYSiR+8GMtzhBzAuJVueK2fD2NrBQCMI+QAYBwhBwDj2CMHcky6J+gsn5yDtwg5kGPSia/1k3PwFlsrAGAcIQcA4wg5ABhHyAHAOEIOAMYRcgAwjpADgHGEHACMI+QAYBwhBwDjCDkAGEfIAcA4Qg4AxhFyADCOkAOAcYQcAIwj5ABgHCEHAOMIOQAYR8gBwDhCDgDGEXIAMI6QA4BxQb8HgH3xeJlOnz6V8nmRyJSUz4nFZimR+MGLsYDbBiHHqKUT3pKSIl240DMG0wC3n7RCvnXrVnV0dGjixImSpN27d+ubb75RfX29gsGgysvLVVtbm9VBAQC3llbIL168qPr6eoXDYUnSpUuXVFdXp4aGBhUUFKimpkbxeFyxWCyrwwIAbpbWyc6enh5NnTp16OvDhw8rHo8rFAopEAiosrJSiUQia0MCAIaX9h75+vXrFQgEtHHjRv3xxx8qLS0deiwajaq1tTWjbxwOhzJ6vl9KSor8HmHcYC29xXp6y/J6phXyffv2SZI6Ojq0ZcsWVVZWqrCwcOhxx3HkOJldydjV1atk0s3o94w1TtB5h7X0FuvpLQvr6Th5w74Azqi+d911lxYvXqxIJKKOjo6hX29vb1c0Gh3dlACAEUkZ8mvXrqmn5/pPqt7eXh07dkwzZsxQU1OT+vr6NDg4qMbGRlVUVGR9WADAzVJurQwMDGj9+vWaNGmSBgcHtWXLFs2dO1ebNm1SdXW1XNdVVVXVDXvm6XCcvBEPPZaszGkBa+kt1tNbub6e/zRfnuu6ub1RDQD4R9xrBQCMI+QAYBwhBwDjCDkAGEfIAcA4Qg4AxhFyADCOkAOAcYQcAIwj5Miq8+fPa8OGDX6PAaT0448/6uWXX/Z7jBEh5ABgHCEHAOPS/oSg28n27dvV1tamK1euaNu2bVqyZInfI5mWTCb1+uuv69SpU0omk3rjjTd0zz33+D2WWd9//7127dol13U1d+5cs9sBueDzzz9XQ0ODioqKNGPGDL/HGTkXN2lpaXFd13XPnTvnrl271udpbGtra3Pnzp3rnjhxwnVd1z127Ji7ceNGn6eyq62tzV21apXb3d3tuq7rXr582eeJ7Dp79qy7Zs0a98qVK67ruu5rr73mbt++3eepRoatlVu4ePGi3nnnHb311lvq7Oz0exzzZs2apdmzZ0uSysrKdP78eZ8nsuvQoUN68sknhz4MvaCgwOeJ7Dpy5IiWL1+uSZMmSZLpD8ch5H9z+PBhvf/++3riiSf09ttvZ/xZpLjZhAkTbvg6Ly+3b+Cfy65evapgkB1RL1y9elWBQGDo64GBAR+nGR0q9TctLS168MEHFYvF1NzcrGQy6fdI5h0/fnzoVfhXX32lefPm+TyRXffff7+++OILXb58WdL1j1/EyMybN08HDhzQwMCAXNdVY2Oj3yONGD/a/+bxxx/X1q1b9fXXX2vRokU3/MTGyMRiMb333ns6f/68CgsLtWPHDr9HMisWi2nt2rWqqanR5MmTVVZWphdeeMHvsUxasGCBFi5cqKqqKk2bNk0LFy7UuXPn/B5rRPioNwAwjq0VADCOkAOAcYQcAIwj5ABgHCEHAOMIOQAYR8gBwDhCDgDGEXIAMO7/AFSLFCRslAunAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 이상값 확인\n",
    "plot_data = [a, b, c, d]\n",
    "plt.boxplot(plot_data)\n",
    "plt.xticks([1,2,3,4],['a', 'b', 'c', 'd'])\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8766f8f1-6011-43a0-8539-d473fe18c021",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 데이터프레임 변환\n",
    "check = ['a', 'b', 'c', 'd']\n",
    "check = [check[j] for j in range(4) for i in range(8)]\n",
    "data = pd.DataFrame({'점수':a+b+c+d, '교육기관':check})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b7fc5680-6a69-4e1f-87de-0c1a7202f6b1",
   "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>점수</th>\n",
       "      <th>교육기관</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>66</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>74</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>82</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>75</td>\n",
       "      <td>a</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   점수 교육기관\n",
       "0  66    a\n",
       "1  74    a\n",
       "2  82    a\n",
       "3  75    a"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "b259316e-73ec-4c2c-9bde-cf0268644b28",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "정규성 검정 결과 ( 샤피로 )\n",
      "ShapiroResult(statistic=0.9516036510467529, pvalue=0.7273865938186646)\n",
      "ShapiroResult(statistic=0.9632648825645447, pvalue=0.8405646681785583)\n",
      "ShapiroResult(statistic=0.8988886475563049, pvalue=0.2823662757873535)\n",
      "ShapiroResult(statistic=0.93294757604599, pvalue=0.5432827472686768)\n",
      "--------------------------------------------------------------------------------\n",
      "정규성 검정 결과 ( 콜모고로프-스미르노프 )\n",
      "KstestResult(statistic=0.625, pvalue=0.08702408702408701)\n",
      "KstestResult(statistic=0.625, pvalue=0.08702408702408701)\n",
      "KstestResult(statistic=0.5, pvalue=0.2826728826728826)\n",
      "KstestResult(statistic=0.25, pvalue=0.98010878010878)\n",
      "KstestResult(statistic=0.25, pvalue=0.98010878010878)\n",
      "KstestResult(statistic=0.25, pvalue=0.98010878010878)\n",
      "--------------------------------------------------------------------------------\n",
      "등분산성 검정 결과 ( levene )\n",
      "LeveneResult(statistic=0.07560756075607561, pvalue=0.972627101462453)\n",
      "--------------------------------------------------------------------------------\n",
      "등분산성 검정 결과 ( bartlett )\n",
      "BartlettResult(statistic=0.07373983712569482, pvalue=0.9947906337756077)\n"
     ]
    }
   ],
   "source": [
    "## 정규성 검정 : 자료의 모집단 분포는 정규분포를 따름 \n",
    "from scipy.stats import shapiro\n",
    "from scipy.stats import levene\n",
    "from scipy.stats import bartlett\n",
    "\n",
    "normal_a = shapiro(a)\n",
    "normal_b = shapiro(b)\n",
    "normal_c = shapiro(c)\n",
    "normal_d = shapiro(d)\n",
    "\n",
    "print('정규성 검정 결과 ( 샤피로 )')\n",
    "print(normal_a, normal_b, normal_c, normal_d, sep='\\n')\n",
    "\n",
    "print('-'*80)\n",
    "\n",
    "print('정규성 검정 결과 ( 콜모고로프-스미르노프 )')\n",
    "print(stats.ks_2samp(a, b), stats.ks_2samp(a, c), stats.ks_2samp(a, d),\n",
    "      stats.ks_2samp(b, c), stats.ks_2samp(b, d),\n",
    "      stats.ks_2samp(c, d), sep='\\n')\n",
    "\n",
    "print('-'*80)\n",
    "\n",
    "## 등분산성 검정 : 모든 집단의 모분산은 동일함\n",
    "print('등분산성 검정 결과 ( levene )')\n",
    "print(levene(a, b, c, d))\n",
    "\n",
    "print('-'*80)\n",
    "\n",
    "print('등분산성 검정 결과 ( bartlett )')\n",
    "print(bartlett(a,b,c,d))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "798757b4-2e1f-4572-80ca-9e7a1aff17c6",
   "metadata": {
    "tags": []
   },
   "source": [
    "*[검정 결과 ]*\n",
    "\n",
    "a, b, c, d 모든 데이터가 p-value > 0.05 보다 크므로 정규성과 등분산성 만족"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "95bee825-da8a-46c7-8873-da19f623fb5c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[F = 3.5], [p = 0.028896983449837425]\n",
      "귀무가설 기각, 4개의 교육훈련 기법 간 차이가 있다.\n"
     ]
    }
   ],
   "source": [
    "# F통계량을 이용한 가설 검정\n",
    "f_statistic, pval = stats.f_oneway(a,b,c,d)\n",
    "\n",
    "print(f'[F = {f_statistic:.1f}], [p = {pval}]')\n",
    "\n",
    "if pval < 0.05:\n",
    "    print('귀무가설 기각, 4개의 교육훈련 기법 간 차이가 있다.')\n",
    "else:\n",
    "    print('귀무가설 채택, 4개의 교육훈련 기법 간 차이가 없다.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "257644f8-c3d5-4d81-b9d1-22c654d08cb8",
   "metadata": {},
   "source": [
    "---\n",
    "Statsmodel을 이용한 가설 검정\n",
    "- model = ols('종속변수 ~ 독립변수', data).fit() : model\n",
    "- sm.stats.anova_lm(model, type=2) : Type 2 Anova DataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "ca92d416-6682-4d98-a9b3-8678871836ed",
   "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>df</th>\n",
       "      <th>sum_sq</th>\n",
       "      <th>mean_sq</th>\n",
       "      <th>F</th>\n",
       "      <th>PR(&gt;F)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>C(교육기관)</th>\n",
       "      <td>3.0</td>\n",
       "      <td>885.25</td>\n",
       "      <td>295.083333</td>\n",
       "      <td>3.482176</td>\n",
       "      <td>0.028897</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Residual</th>\n",
       "      <td>28.0</td>\n",
       "      <td>2372.75</td>\n",
       "      <td>84.741071</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            df   sum_sq     mean_sq         F    PR(>F)\n",
       "C(교육기관)    3.0   885.25  295.083333  3.482176  0.028897\n",
       "Residual  28.0  2372.75   84.741071       NaN       NaN"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Statsmodel을 사용한 일원분산분석\n",
    "from statsmodels.formula.api import ols\n",
    "from statsmodels.stats.anova import anova_lm\n",
    "\n",
    "model_1A = ols('점수 ~ C(교육기관)', data).fit()\n",
    "anova_lm_result = anova_lm(model_1A)\n",
    "anova_lm_result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "d40119c3-d8d0-41f7-af7b-380afa64f641",
   "metadata": {},
   "outputs": [],
   "source": [
    "# data.shape(32, 2)\n",
    "# SSA k-1 = 4-1 = 3\n",
    "# SSE N-k = 32-4 = 28\n",
    "# SST N-1 = 32-1 = 31\n",
    "# MSA / MSW = 295.08333 / 84.731071 = 3.482176"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "bae0d3e4-9d98-4972-9fa7-a4a4c28cbe56",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "귀무가설 기각, 4개의 교육훈련 기법 간 차이가 있다.\n"
     ]
    }
   ],
   "source": [
    "pval = anova_lm_result['PR(>F)'][0]\n",
    "if pval < 0.05:\n",
    "    print('귀무가설 기각, 4개의 교육훈련 기법 간 차이가 있다.')\n",
    "else:\n",
    "    print('귀무가설 채택, 4개의 교육훈련 기법 간 차이가 없다.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "98fac050-efb1-4ae3-adca-66a2d3edfe17",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Multiple Comparison of Means - Tukey HSD, FWER=0.05</caption>\n",
       "<tr>\n",
       "  <th>group1</th> <th>group2</th> <th>meandiff</th>  <th>p-adj</th>   <th>lower</th>   <th>upper</th>  <th>reject</th>\n",
       "</tr>\n",
       "<tr>\n",
       "     <td>a</td>      <td>b</td>    <td>-13.625</td> <td>0.0298</td> <td>-26.1919</td> <td>-1.0581</td>  <td>True</td> \n",
       "</tr>\n",
       "<tr>\n",
       "     <td>a</td>      <td>c</td>    <td>-11.875</td> <td>0.0692</td> <td>-24.4419</td> <td>0.6919</td>   <td>False</td>\n",
       "</tr>\n",
       "<tr>\n",
       "     <td>a</td>      <td>d</td>     <td>-9.5</td>   <td>0.1896</td> <td>-22.0669</td> <td>3.0669</td>   <td>False</td>\n",
       "</tr>\n",
       "<tr>\n",
       "     <td>b</td>      <td>c</td>     <td>1.75</td>   <td>0.9809</td> <td>-10.8169</td> <td>14.3169</td>  <td>False</td>\n",
       "</tr>\n",
       "<tr>\n",
       "     <td>b</td>      <td>d</td>     <td>4.125</td>  <td>0.8068</td>  <td>-8.4419</td> <td>16.6919</td>  <td>False</td>\n",
       "</tr>\n",
       "<tr>\n",
       "     <td>c</td>      <td>d</td>     <td>2.375</td>  <td>0.9546</td> <td>-10.1919</td> <td>14.9419</td>  <td>False</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## 사후 분석 : ANOVA 검증 결과 유의미하다는 결론을 얻었을 때,\n",
    "#             구체적으로 어떤 수준(들)에서 평균 차이가 나는지를 검증\n",
    "from statsmodels.stats.multicomp import pairwise_tukeyhsd\n",
    "\n",
    "HSD = pairwise_tukeyhsd(data['점수'], data['교육기관'], alpha=0.05)\n",
    "HSD.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c84aa8cf-fec4-4198-8cae-61ffeb7473fa",
   "metadata": {},
   "source": [
    "***[결론 ]***\n",
    "- 교육기관에 따른 weight의 평균 차이는 유의미하였다(F(3, 28) = 4.221093, p < 0.05)\n",
    "- Tukey의 HSD를 이용하여 사후 분석을 실시한 결과, a조건과 b조건에서 유의미한 평균 차이가 있었다 ( p < 0.05 )\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "725ad5f3-6be9-4435-9628-e5785e69caf2",
   "metadata": {},
   "source": [
    "### 이원분산분석( Two-way ANOVA )\n",
    "- 종속변수가 1개, 독립변수가 2개이면서 독립변수 각각이 factor(요인)로서 factor 내에 lavel(수준)을 가지고 있을 때, 각 집단의 평균이 유의미한 차이가 있는지 유무를 확인\n",
    "---\n",
    "- 증명하고자 하는 가설\n",
    "  - 귀무가설 : 강도에 대한 각 판유리와 온도 간에는 상호작용 효과가 없다.\n",
    "  - 대립가설 : 강도에 대한 각 판유리와 온도 간에는 상호작용 효과가 있다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "926098b9-12e5-4032-a54d-57834fbe7a5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 종속변수 : 강도, 독립변수 : 판유리/온도\n",
    "dat = pd.DataFrame({'판유리':['유리1', '유리1', '유리1', '유리2', '유리2', '유리2', '유리3', '유리3',\n",
    "                             '유리3', '유리1', '유리1', '유리1', '유리2', '유리2', '유리2', '유리3',\n",
    "                             '유리3', '유리3', '유리1', '유리1', '유리1', '유리2', '유리2', '유리2',\n",
    "                             '유리3', '유리3', '유리3'],\n",
    "                    '온도':[100, 100, 100, 100, 100, 100, 100, 100, 100, 125, 125, 125, 125, 125,\n",
    "                            125, 125, 125, 125, 150, 150, 150, 150, 150, 150, 150, 150, 150],\n",
    "                    '강도':[580, 568, 570, 550, 530, 579, 546, 575, 599, 1090, 1087, 1085, 1070,\n",
    "                           1035, 1000, 1045, 1053, 1066, 1392, 1380, 1385, 1328, 1312, 1299, 867,\n",
    "                           904, 889]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "f2888240-f59f-4c5a-a63e-4534e76e9f49",
   "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>판유리</th>\n",
       "      <th>온도</th>\n",
       "      <th>강도</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>유리1</td>\n",
       "      <td>100</td>\n",
       "      <td>580</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>유리1</td>\n",
       "      <td>100</td>\n",
       "      <td>568</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>유리1</td>\n",
       "      <td>100</td>\n",
       "      <td>570</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>유리2</td>\n",
       "      <td>100</td>\n",
       "      <td>550</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>유리2</td>\n",
       "      <td>100</td>\n",
       "      <td>530</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   판유리   온도   강도\n",
       "0  유리1  100  580\n",
       "1  유리1  100  568\n",
       "2  유리1  100  570\n",
       "3  유리2  100  550\n",
       "4  유리2  100  530"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dat.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "a0cb5ac4-9815-445a-a976-28cd2423420c",
   "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></th>\n",
       "      <th>강도</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>판유리</th>\n",
       "      <th>온도</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">유리1</th>\n",
       "      <th>100</th>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>125</th>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">유리2</th>\n",
       "      <th>100</th>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>125</th>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">유리3</th>\n",
       "      <th>100</th>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>125</th>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         강도\n",
       "판유리 온도     \n",
       "유리1 100   3\n",
       "    125   3\n",
       "    150   3\n",
       "유리2 100   3\n",
       "    125   3\n",
       "    150   3\n",
       "유리3 100   3\n",
       "    125   3\n",
       "    150   3"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dat.groupby(['판유리', '온도']).agg(len)\n",
    "# 판유리 온도 요인으로 구분한 각 집단별 표본 수는 모두 3으로 동일"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "70f31719-72f9-46b5-9e64-0051ddf83b6d",
   "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>df</th>\n",
       "      <th>sum_sq</th>\n",
       "      <th>mean_sq</th>\n",
       "      <th>F</th>\n",
       "      <th>PR(&gt;F)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>C(판유리)</th>\n",
       "      <td>2.0</td>\n",
       "      <td>1.507143e+05</td>\n",
       "      <td>75357.148148</td>\n",
       "      <td>206.144174</td>\n",
       "      <td>3.922982e-13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(온도)</th>\n",
       "      <td>2.0</td>\n",
       "      <td>1.969825e+06</td>\n",
       "      <td>984912.259259</td>\n",
       "      <td>2694.288855</td>\n",
       "      <td>5.025167e-23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(판유리):C(온도)</th>\n",
       "      <td>4.0</td>\n",
       "      <td>2.903206e+05</td>\n",
       "      <td>72580.148148</td>\n",
       "      <td>198.547518</td>\n",
       "      <td>1.264106e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Residual</th>\n",
       "      <td>18.0</td>\n",
       "      <td>6.580000e+03</td>\n",
       "      <td>365.555556</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                df        sum_sq        mean_sq            F        PR(>F)\n",
       "C(판유리)         2.0  1.507143e+05   75357.148148   206.144174  3.922982e-13\n",
       "C(온도)          2.0  1.969825e+06  984912.259259  2694.288855  5.025167e-23\n",
       "C(판유리):C(온도)   4.0  2.903206e+05   72580.148148   198.547518  1.264106e-14\n",
       "Residual      18.0  6.580000e+03     365.555556          NaN           NaN"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 이원분산분석\n",
    "from statsmodels.formula.api import ols\n",
    "from statsmodels.stats.anova import anova_lm\n",
    "\n",
    "model = ols('강도 ~ C(판유리) * C(온도)', dat).fit()\n",
    "anova_result = pd.DataFrame(anova_lm(model))\n",
    "anova_result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55884b23-0c49-459f-a6b7-07f3eb19a5df",
   "metadata": {},
   "source": [
    "***[결과 해석 ]***\n",
    "- 판유리 : F(2, 18) = 206.370554, p < 0.05로 유의미. 판유리의 수준에 따라 평균에 차이가 난다고 볼 수 있음\n",
    "- 온도 : F(2, 18) = 2695.259499, p < 0.05로 유의미. 온도의 수준에 따라 평균에 차이가 난다고 볼 수 있음\n",
    "- 판유리:온도 : F(4, 18) = 198.725707, p < 0.05로 유의미. 상호작용 효과는 발견"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "efcdb869-c912-4939-b9e5-1d65ff8c7ef7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Multiple Comparison of Means - Tukey HSD, FWER=0.05    \n",
      "===========================================================\n",
      " group1  group2  meandiff p-adj    lower     upper   reject\n",
      "-----------------------------------------------------------\n",
      "유리1/100 유리1/125  514.6667    0.0  459.9678  569.3655   True\n",
      "유리1/100 유리1/150     813.0    0.0  758.3012  867.6988   True\n",
      "유리1/100 유리2/100  -19.6667 0.9307  -74.3655   35.0322  False\n",
      "유리1/100 유리2/125  462.3333    0.0  407.6345  517.0322   True\n",
      "유리1/100 유리2/150  740.3333    0.0  685.6345  795.0322   True\n",
      "유리1/100 유리3/100    0.6667    1.0  -54.0322   55.3655  False\n",
      "유리1/100 유리3/125     482.0    0.0  427.3012  536.6988   True\n",
      "유리1/100 유리3/150     314.0    0.0  259.3012  368.6988   True\n",
      "유리1/125 유리1/150  298.3333    0.0  243.6345  353.0322   True\n",
      "유리1/125 유리2/100 -534.3333    0.0 -589.0322 -479.6345   True\n",
      "유리1/125 유리2/125  -52.3333  0.067 -107.0322    2.3655  False\n",
      "유리1/125 유리2/150  225.6667    0.0  170.9678  280.3655   True\n",
      "유리1/125 유리3/100    -514.0    0.0 -568.6988 -459.3012   True\n",
      "유리1/125 유리3/125  -32.6667 0.5066  -87.3655   22.0322  False\n",
      "유리1/125 유리3/150 -200.6667    0.0 -255.3655 -145.9678   True\n",
      "유리1/150 유리2/100 -832.6667    0.0 -887.3655 -777.9678   True\n",
      "유리1/150 유리2/125 -350.6667    0.0 -405.3655 -295.9678   True\n",
      "유리1/150 유리2/150  -72.6667 0.0048 -127.3655  -17.9678   True\n",
      "유리1/150 유리3/100 -812.3333    0.0 -867.0322 -757.6345   True\n",
      "유리1/150 유리3/125    -331.0    0.0 -385.6988 -276.3012   True\n",
      "유리1/150 유리3/150    -499.0    0.0 -553.6988 -444.3012   True\n",
      "유리2/100 유리2/125     482.0    0.0  427.3012  536.6988   True\n",
      "유리2/100 유리2/150     760.0    0.0  705.3012  814.6988   True\n",
      "유리2/100 유리3/100   20.3333  0.918  -34.3655   75.0322  False\n",
      "유리2/100 유리3/125  501.6667    0.0  446.9678  556.3655   True\n",
      "유리2/100 유리3/150  333.6667    0.0  278.9678  388.3655   True\n",
      "유리2/125 유리2/150     278.0    0.0  223.3012  332.6988   True\n",
      "유리2/125 유리3/100 -461.6667    0.0 -516.3655 -406.9678   True\n",
      "유리2/125 유리3/125   19.6667 0.9307  -35.0322   74.3655  False\n",
      "유리2/125 유리3/150 -148.3333    0.0 -203.0322  -93.6345   True\n",
      "유리2/150 유리3/100 -739.6667    0.0 -794.3655 -684.9678   True\n",
      "유리2/150 유리3/125 -258.3333    0.0 -313.0322 -203.6345   True\n",
      "유리2/150 유리3/150 -426.3333    0.0 -481.0322 -371.6345   True\n",
      "유리3/100 유리3/125  481.3333    0.0  426.6345  536.0322   True\n",
      "유리3/100 유리3/150  313.3333    0.0  258.6345  368.0322   True\n",
      "유리3/125 유리3/150    -168.0    0.0 -222.6988 -113.3012   True\n",
      "-----------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "## 사후 분석\n",
    "dat['온도'] = dat['온도'].astype('str')\n",
    "dat['combi'] = dat.판유리 + '/' + dat.온도\n",
    "print(pairwise_tukeyhsd(dat['강도'], dat['combi']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c4b1010-0aaf-4a30-9012-8b1aafa9c5cd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}