Keras_GLMM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalized linear mixed model (GLMM) với Keras\n",
    "\n",
    "### BS. Lê Ngọc Khả Nhi\n",
    "\n",
    "# Giới thiệu\n",
    "\n",
    "Thân chào các bạn, trong bài thực hành hôm nay Nhi làm một thí nghiệm với mục tiêu dùng Neural network như công cụ suy diễn thống kê với tính năng như một mô hình Generalized linear mixed model (GLMM).\n",
    "\n",
    "## Tóm tắt về GLMM\n",
    "\n",
    "Trước hết, Nhi tóm tắt lại khái niệm về mô hình GLMM:\n",
    "\n",
    "Như ta biết, mô hình tuyến tính tổng quát (GLM) có dạng: $\\mathbf{y} = \\mathbf{X}\\beta + \\epsilon$, trong đó $\\mathbf{X}$ tương ứng với tập hợp các hiệu ứng chính của các hiệp biến/biến độc lập đối với variance của biến kết quả/phụ thuộc $\\mathbf{y}$. Những hiệu ứng chính này được khảo sát qua giá trị của các coefficients $\\beta$, còn $\\epsilon$ tương ứng với giá trị nhiễu/sai số ngẫu nhiên, với giả định là $\\epsilon\\sim\\mathcal{N}(\\mathbf{0}, \\mathbf{I})$\n",
    "\n",
    "Mô hình tuyến tính hỗn hợp (Mixed linear model hay GLMM) là một sự mở rộng của mô hình GLM, bằng cách đưa thêm vào một nhóm hiệu ứng ngẫu nhiên (random effect) $\\mathbf{Z}\\eta$, trong đó $\\mathbf{Z}$ tương ứng với tập hợp các yếu tố ngẫu nhiên (phân nhóm), đo bằng coefficients $\\eta$. Lúc này nội dung mô hình GLMM trở thành: \\begin{aligned} \\eta &\\sim \\mathcal{N}(\\mathbf{0}, \\sigma^2 \\mathbf{I}), \\\\ \\mathbf{y} &= \\mathbf{X}\\beta + \\mathbf{Z}\\eta + \\epsilon.\\end{aligned}\n",
    "\n",
    "Cách suy luận thống kê với GLMM cũng tương tự GLM, đó là dựa vào dữ liệu, ta sẽ xác định giá trị của $\\beta$, $\\eta$, và $\\sigma^2$, trong đó $\\beta$ được gọi là 'hiệu ứng chinh' (fixed effects), còn $\\eta$ được gọi là hiệu ứng ngẫu nhiên hay bất định (random/latent effect). Trung bình của biến kết quả $\\mathbf{y}$ được ước lượng chủ yếu bằng các hiệu ứng chính $\\mathbf{X}\\beta$, còn nhiệm vụ của $\\eta$ là giải thích càng nhiều càng tốt sai số ngẫu nhiên, vốn được cho là nhiễu trong $\\epsilon$ của mô hình GLM. \n",
    "\n",
    "## Dùng neural network mô phỏng GLMM\n",
    "\n",
    "Trong các bài hoc nhập môn Deep learning (với tensorflow hay Keras), chúng ta đã thấy có thể dùng Neural network như một mô hình GLM, bao gồm logistic model. Như vậy, ta hoàn toàn có thể dùng Neural network như một mô hình GLMM vì so với GLM, GLMM chỉ có thêm 1 bộ phận Z bên cạnh X. \n",
    "\n",
    "Trong bài này, ta sẽ dùng Keras (tensorflow backend) để dựng một neural network có tính năng như một mô hình GLMM\n",
    "\n",
    "Ta dùng dữ liệu minh họa từ nghiên cứu 'University Lecture/Instructor Evaluations by Students at ETH', được trích xuất từ Rpackage lme4.\n",
    "\n",
    "Nội dung dữ liệu như sau:\n",
    "\n",
    "'''\n",
    "University lecture evaluations by students at ETH Zurich, anonymized for privacy protection. This is an interesting “medium” sized example of a partially nested mixed effect model.\n",
    "\n",
    "A data frame with 73421 observations on the following 7 variables.\n",
    "\n",
    "s = a factor with levels 1:2972 denoting individual students.\n",
    "\n",
    "d = a factor with 1128 levels from 1:2160, denoting individual professors or lecturers.\n",
    "\n",
    "studage = an ordered factor with levels 2 < 4 < 6 < 8, denoting student's “age” measured in the semester number the student has been enrolled.\n",
    "\n",
    "lectage\n",
    "an ordered factor with 6 levels, 1 < 2 < ... < 6, measuring how many semesters back the lecture rated had taken place.\n",
    "\n",
    "service = a binary factor with levels 0 and 1; a lecture is a “service”, if held for a different department than the lecturer's main one.\n",
    "\n",
    "dept = a factor with 14 levels from 1:15, using a random code for the department of the lecture.\n",
    "\n",
    "y = a numeric vector of ratings of lectures by the students, using the discrete scale 1:5, with meanings of ‘poor’ to ‘very good’.\n",
    "\n",
    "Each observation is one student's rating for a specific lecture (of one lecturer, during one semester in the past).\n",
    "\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from patsy import dmatrices\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import keras\n",
    "from keras.models import Model\n",
    "from keras.layers import Input, Add, Dense\n",
    "from keras import backend as K\n",
    "from keras.callbacks import EarlyStopping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ta tải dữ liệu về và tạo dataframe df có nội dung như sau:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 73421 entries, 1 to 73421\n",
      "Data columns (total 7 columns):\n",
      " #   Column    Non-Null Count  Dtype   \n",
      "---  ------    --------------  -----   \n",
      " 0   Student   73421 non-null  category\n",
      " 1   Lecturer  73421 non-null  category\n",
      " 2   Stud_Age  73421 non-null  object  \n",
      " 3   Lect_Age  73421 non-null  object  \n",
      " 4   Service   73421 non-null  object  \n",
      " 5   Dept      73421 non-null  category\n",
      " 6   Score     73421 non-null  float64 \n",
      "dtypes: category(3), float64(1), object(3)\n",
      "memory usage: 3.3+ MB\n"
     ]
    }
   ],
   "source": [
    "df = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/lme4/InstEval.csv', index_col=0)\n",
    "\n",
    "df.columns = [\"Student\",\"Lecturer\",\"Stud_Age\",\"Lect_Age\",\"Service\",\"Dept\",\"Score\"]\n",
    "\n",
    "df[\"Student\"] = df[\"Student\"].astype('category')\n",
    "df[\"Lecturer\"] = df[\"Lecturer\"].astype('category')\n",
    "df[\"Dept\"] = df[\"Dept\"].astype('category')\n",
    "\n",
    "df[\"Stud_Age\"] = df[\"Stud_Age\"].astype('object')\n",
    "df[\"Lect_Age\"] = df[\"Lect_Age\"].astype('object')\n",
    "df[\"Service\"] = df[\"Service\"].astype('object')\n",
    "\n",
    "df[\"Score\"] = df[\"Score\"].astype('float')\n",
    "\n",
    "df.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "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>Student</th>\n",
       "      <th>Lecturer</th>\n",
       "      <th>Stud_Age</th>\n",
       "      <th>Lect_Age</th>\n",
       "      <th>Service</th>\n",
       "      <th>Dept</th>\n",
       "      <th>Score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1002</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>1050</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>6</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>1582</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>2050</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>2</td>\n",
       "      <td>115</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Student Lecturer Stud_Age Lect_Age Service Dept  Score\n",
       "1       1     1002        2        2       0    2    5.0\n",
       "2       1     1050        2        1       1    6    2.0\n",
       "3       1     1582        2        2       0    2    5.0\n",
       "4       1     2050        2        2       1    3    3.0\n",
       "5       2      115        2        1       0    5    2.0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mục tiêu giả định là dựng một mô hình GLMM với nội dung :\n",
    "    \n",
    "'Score ~ Stud_Age + Lect_Age + Service + (1|Student) + (1|Lecturer) + (1|Dept), data = df'\n",
    "\n",
    "Sau đây là kết quả của mô hình dựng bằng lme4 trong R:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### R codes\n",
    "\n",
    "library('lme4')\n",
    "\n",
    "mod = lmer(Score ~ Stud_Age + Lect_Age + Service + (1|Student) + (1|Lecturer) + (1|Dept), data = df)\n",
    "\n",
    "summary(mod)\n",
    "\n",
    "\n",
    "#### OUTPUT\n",
    "\n",
    "Linear mixed model fit by REML ['lmerMod']\n",
    "Formula: Score ~ Stud_Age + Lect_Age + Service + (1 | Student) + (1 |      Lecturer) + (1 | Dept)\n",
    "   Data: data\n",
    "\n",
    "REML criterion at convergence: 237615.8\n",
    "\n",
    "Scaled residuals: \n",
    "     Min       1Q   Median       3Q      Max \n",
    "-3.09054 -0.74589  0.04067  0.77162  3.15582 \n",
    "\n",
    "Random effects:\n",
    " Groups   Name        Variance Std.Dev.\n",
    " Student  (Intercept) 0.106724 0.32669 \n",
    " Lecturer (Intercept) 0.260789 0.51067 \n",
    " Dept     (Intercept) 0.006926 0.08323 \n",
    " Residual             1.383443 1.17620 \n",
    "Number of obs: 73421, groups:  Student, 2972; Lecturer, 1128; Dept, 14\n",
    "\n",
    "Fixed effects:\n",
    "            Estimate Std. Error t value\n",
    "(Intercept)  3.29077    0.03325  98.971\n",
    "Stud_Age4    0.05197    0.02318   2.242\n",
    "Stud_Age6    0.07216    0.02398   3.010\n",
    "Stud_Age8    0.13634    0.02639   5.166\n",
    "Lect_Age2   -0.08074    0.01538  -5.249\n",
    "Lect_Age3   -0.11019    0.01672  -6.588\n",
    "Lect_Age4   -0.18917    0.01960  -9.650\n",
    "Lect_Age5   -0.16443    0.02141  -7.681\n",
    "Lect_Age6   -0.24604    0.02049 -12.009\n",
    "Service1    -0.07278    0.01348  -5.399\n",
    "\n",
    "Correlation of Fixed Effects:\n",
    "          (Intr) Std_A4 Std_A6 Std_A8 Lct_A2 Lct_A3 Lct_A4 Lct_A5 Lct_A6\n",
    "Stud_Age4 -0.322                                                        \n",
    "Stud_Age6 -0.353  0.517                                                 \n",
    "Stud_Age8 -0.346  0.491  0.575                                          \n",
    "Lect_Age2 -0.195 -0.009 -0.023 -0.025                                   \n",
    "Lect_Age3 -0.034 -0.197 -0.236 -0.259  0.342                            \n",
    "Lect_Age4 -0.057 -0.170 -0.209 -0.234  0.473  0.403                     \n",
    "Lect_Age5  0.000 -0.136 -0.274 -0.283  0.300  0.417  0.374              \n",
    "Lect_Age6  0.004 -0.144 -0.297 -0.373  0.406  0.447  0.498  0.490       \n",
    "Service1  -0.142  0.017  0.056  0.054 -0.033 -0.044 -0.063 -0.072 -0.121"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GLMM neural network bằng Keras"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nhi viết một class Keras_GLMM với tính năng dựng một neural network bằng keras và tensorflow từ 2 công thức của fixed_effect và random_effects, sau đó fit model này trên dữ liệu đầu vào. Mô hình sau khi huấn luyện có thể xuất ra nội dung của betas trong fixed effects với cấu trúc gần giống như một mô hình GLMM mà ta thấy trong R.\n",
    "\n",
    "Cơ chế hoạt động như sau:\n",
    "\n",
    "Từ 2 formulas cho Fixed và Random effects, ta dùng hàm dmatrices của package patsy để dựng design matrices X,Y,Z và danh sách tên của fixed effects trong X.\n",
    "\n",
    "Mô hình neural network GLMM gồm 2 input layers riêng, tiếp nhận X và Z, mỗi input layer này nối trực tiếp với 1 output neuron để ước lượng fixed và random effects, 2 neuron này sau đó sẽ được cộng chung với nhau để xuất kết quả Y qua 1 output neuron sau cùng.\n",
    "\n",
    "Cấu trúc này được huấn luyện theo cơ chế 'SGD' nhằm tối thiểu giá trị hàm loss là 'mean_squared_error', \n",
    "\n",
    "Lịch sử huấn luyện có thể xem lại qua biểu đồ learning curve\n",
    "\n",
    "Cuối cùng, method Summary cho phép xuất giá trị betas và etas\n",
    "\n",
    "(Code này dựa trên một bài giảng của anh Junpeng Lao, một kỹ sư tại Google)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "run_control": {
     "marked": false
    }
   },
   "outputs": [],
   "source": [
    "class Keras_GLMM:\n",
    "    \n",
    "    def __init__(self, data, fe_form, re_form):\n",
    "        self.data = data\n",
    "        self.fe_form = fe_form\n",
    "        self.re_form = re_form\n",
    "        \n",
    "        self.Y, self.X = dmatrices(self.fe_form, \n",
    "                                   data=self.data, \n",
    "                                   return_type='matrix')\n",
    "        \n",
    "        self.Terms = self.X.design_info.column_names\n",
    "        \n",
    "        self.X = np.asarray(self.X)\n",
    "        self.Y = np.asarray(self.Y).flatten()\n",
    "        \n",
    "        _, self.Z = dmatrices(self.re_form, \n",
    "                              data=self.data, \n",
    "                              return_type='matrix')\n",
    "        \n",
    "        self.Z = np.asarray(self.Z)\n",
    "        \n",
    "        self.n_fixed = np.shape(self.X)\n",
    "        self.n_randm = np.shape(self.Z)\n",
    "        \n",
    "    def Model_build(self):\n",
    "        K.clear_session()\n",
    "\n",
    "        fixedpred = np.argmax(self.X,axis=1)\n",
    "        randmpred = np.argmax(self.Z,axis=1)\n",
    "\n",
    "        X_input = Input(batch_shape=(None, self.n_fixed[1]-1))\n",
    "        fixed_comp = Dense(1, input_dim = self.n_fixed[1]-1, name = 'Fixed_effect')(X_input)\n",
    "\n",
    "        Z_input  = Input(batch_shape=(None, self.n_randm[1]))\n",
    "        randm_comp = Dense(1, input_dim = self.n_randm[1], use_bias=None, name = 'Random_Effect')(Z_input)\n",
    "\n",
    "        mixed = keras.layers.add([fixed_comp, randm_comp])\n",
    "\n",
    "        self.model = Model([X_input, Z_input], mixed)\n",
    "\n",
    "        self.model.compile(loss='mean_squared_error', \n",
    "                           optimizer='sgd')\n",
    "        return self.model\n",
    "    \n",
    "    def Model_fit(self, n_epoch = 100, batch_size = 500, patience = 30):\n",
    "        try:\n",
    "            early_stop = EarlyStopping(monitor = 'loss', patience = patience)\n",
    "\n",
    "            self.hist = self.model.fit([self.X[:,1:], self.Z], self.Y.flatten(), \n",
    "                                  epochs = n_epoch, \n",
    "                                  batch_size = batch_size, \n",
    "                                  verbose=0,\n",
    "                                  shuffle=True,\n",
    "                                  callbacks=[early_stop])\n",
    "        except:\n",
    "            print('Mô hình chưa được compile, hãy dùng method Model_build')\n",
    "    \n",
    "    def history_plot(self):\n",
    "        try:\n",
    "            plt.figure(figsize=(10,5))\n",
    "            plt.plot(self.hist.history['loss'],'b-',label='train')\n",
    "            plt.title('Model loss')\n",
    "            plt.ylabel('Mean squared error')\n",
    "            plt.xlabel('epoch')\n",
    "            plt.legend()\n",
    "            plt.show()\n",
    "        except:\n",
    "            print('Mô hình chưa được tạo ra hoặc chưa được huấn luyện')\n",
    "            \n",
    "    def Summary(self):\n",
    "        try:\n",
    "            Ypredict = self.model.predict([self.X[:,1:], self.Z])\n",
    "            betas = np.hstack((self.model.get_weights()[1], \n",
    "                               self.model.get_weights()[0].flatten()))\n",
    "\n",
    "            etas = self.model.get_weights()[2].flatten()\n",
    "\n",
    "            fe_params = pd.DataFrame(betas, index = self.Terms, columns=['Keras_GLMM'])\n",
    "            re_params = pd.Series(etas)            \n",
    "            \n",
    "            self.model.summary()\n",
    "            return fe_params, re_params\n",
    "        \n",
    "        except:\n",
    "            print('Mô hình chưa được tạo ra hoặc chưa được huấn luyện')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ta áp dụng class Keras_GLMM như sau:\n",
    "    \n",
    "Trước hết khai báo 2 công thức fe_form, re_form và dataframe df, ta tạo ra được objet mod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = Keras_GLMM(data = df,  \n",
    "                 fe_form = 'Score ~ Stud_Age + Lect_Age + Service', \n",
    "                 re_form = 'Score ~ -1 + Student+Lecturer+Dept')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do mô hình chưa được khởi tạo, các method history_plot, Summary chưa thể hoạt động:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mô hình chưa được tạo ra hoặc chưa được huấn luyện\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 720x360 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mod.history_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mô hình chưa được tạo ra hoặc chưa được huấn luyện\n"
     ]
    }
   ],
   "source": [
    "mod.Summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ta khởi tạo mô hình bằng method Model_build(), bạn cũng có thể trực tiếp dùng Model_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<keras.engine.training.Model at 0x25b65143048>"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod.Model_build()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Khi mới vừa tạo ra, mô hình chưa được huấn luyện nên nội dung betas là ngẫu nhiên và không sử dụng được:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"model_1\"\n",
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "input_1 (InputLayer)            (None, 9)            0                                            \n",
      "__________________________________________________________________________________________________\n",
      "input_2 (InputLayer)            (None, 4112)         0                                            \n",
      "__________________________________________________________________________________________________\n",
      "Fixed_effect (Dense)            (None, 1)            10          input_1[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "Random_Effect (Dense)           (None, 1)            4112        input_2[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "add_1 (Add)                     (None, 1)            0           Fixed_effect[0][0]               \n",
      "                                                                 Random_Effect[0][0]              \n",
      "==================================================================================================\n",
      "Total params: 4,122\n",
      "Trainable params: 4,122\n",
      "Non-trainable params: 0\n",
      "__________________________________________________________________________________________________\n"
     ]
    },
    {
     "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>Keras_GLMM</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Stud_Age[T.4]</th>\n",
       "      <td>-0.462801</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Stud_Age[T.6]</th>\n",
       "      <td>-0.357625</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Stud_Age[T.8]</th>\n",
       "      <td>-0.198635</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.2]</th>\n",
       "      <td>0.301004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.3]</th>\n",
       "      <td>-0.605316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.4]</th>\n",
       "      <td>0.607983</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.5]</th>\n",
       "      <td>0.147077</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.6]</th>\n",
       "      <td>0.206578</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Service[T.1]</th>\n",
       "      <td>-0.029131</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               Keras_GLMM\n",
       "Intercept        0.000000\n",
       "Stud_Age[T.4]   -0.462801\n",
       "Stud_Age[T.6]   -0.357625\n",
       "Stud_Age[T.8]   -0.198635\n",
       "Lect_Age[T.2]    0.301004\n",
       "Lect_Age[T.3]   -0.605316\n",
       "Lect_Age[T.4]    0.607983\n",
       "Lect_Age[T.5]    0.147077\n",
       "Lect_Age[T.6]    0.206578\n",
       "Service[T.1]    -0.029131"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fe_df, re_df = mod.Summary()\n",
    "\n",
    "fe_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ta bắt đầu huấn luyện (mặc định là 100 epochs với batch_size = 500). Đây là 1 dữ liệu rất lớn với 73,421 trường hợp, nên hy vọng bạn có GPU trong máy của mình. Nếu bạn chỉ có CPU, có thể chạy notebook này trên google colab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.Model_fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mô hình sau huấn luyện được lưu trong mod.model, ta có thể dùng method summary của keras"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"model_1\"\n",
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "input_1 (InputLayer)            (None, 9)            0                                            \n",
      "__________________________________________________________________________________________________\n",
      "input_2 (InputLayer)            (None, 4112)         0                                            \n",
      "__________________________________________________________________________________________________\n",
      "Fixed_effect (Dense)            (None, 1)            10          input_1[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "Random_Effect (Dense)           (None, 1)            4112        input_2[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "add_1 (Add)                     (None, 1)            0           Fixed_effect[0][0]               \n",
      "                                                                 Random_Effect[0][0]              \n",
      "==================================================================================================\n",
      "Total params: 4,122\n",
      "Trainable params: 4,122\n",
      "Non-trainable params: 0\n",
      "__________________________________________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "mod.model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Method history_plot vẽ learning curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAFNCAYAAABFbcjcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de7hddX3n8fc3J8nJFQJJuOVCELGCCEHTiMUpoB0KqIV5qi2o6DjatH2gQmsv6FQdecbWqdVapiCi0oJyGSugjIKCDoqUckkwcgvUCMEEIgmBXEkgl+/8sdb2bE7OOdn7nL3O3ifn/Xqe9ay1f2utvb/HpfDx9/uttSIzkSRJUmcY0+4CJEmS1MNwJkmS1EEMZ5IkSR3EcCZJktRBDGeSJEkdxHAmSZLUQQxnkkatiJgXERkRYxs49r9GxJ1D/R5J2hPDmaQRISJWRMRLETGjV/vSMhjNa09lktRahjNJI8kTwNm1DxHxWmBi+8qRpNYznEkaSb4KvLfu8/uAq+oPiIh9I+KqiFgbEU9GxF9HxJhyX1dE/H1EPBsRjwNv7ePcr0TE6oh4KiL+Z0R0NVtkRBwSETdFxHMRsTwi/qBu38KIWBwRGyPimYj4XNk+ISK+FhHrImJ9RNwXEQc2+9uSRj7DmaSR5G5gn4g4sgxNvw98rdcx/xvYF3gFcCJFmHt/ue8PgLcBxwELgHf0OvdKYAfwyvKYU4APDqLOa4FVwCHlb/xNRLyl3PePwD9m5j7A4cDXy/b3lXXPAaYDfwRsHcRvSxrhDGeSRppa79l/Bh4FnqrtqAtsH8nMTZm5AvgscE55yO8Bn8/MlZn5HPC3deceCJwGXJCZWzJzDfAPwFnNFBcRc4A3AX+Vmdsycynw5boatgOvjIgZmbk5M++ua58OvDIzd2bmkszc2MxvS9o7GM4kjTRfBd4F/Fd6DWkCM4DxwJN1bU8Cs8rtQ4CVvfbVHAqMA1aXw4rrgS8CBzRZ3yHAc5m5qZ8aPgC8Cni0HLp8W93f9T3guoh4OiL+LiLGNfnbkvYChjNJI0pmPklxY8DpwA29dj9L0QN1aF3bXHp611ZTDBvW76tZCbwIzMjMaeWyT2a+pskSnwb2j4ipfdWQmT/LzLMpQt//Ar4REZMzc3tmfjIzjwJ+g2L49b1IGnUMZ5JGog8Ab87MLfWNmbmTYg7XpyJiakQcCvwZPfPSvg58KCJmR8R+wIV1564GbgU+GxH7RMSYiDg8Ik5sprDMXAncBfxtOcn/mLLeqwEi4j0RMTMzdwHry9N2RsTJEfHacmh2I0XI3NnMb0vaOxjOJI04mfnzzFzcz+4/AbYAjwN3AtcAV5T7vkQxdPhT4H5273l7L8Ww6CPA88A3gIMHUeLZwDyKXrQbgU9k5m3lvlOBhyNiM8XNAWdl5jbgoPL3NgLLgB+x+80OkkaByMx21yBJkqSSPWeSJEkdxHAmSZLUQQxnkiRJHcRwJkmS1EEMZ5IkSR1kbLsLaKUZM2bkvHnz2l2GJEnSHi1ZsuTZzJzZu32vCmfz5s1j8eL+Hn0kSZLUOSLiyb7aHdaUJEnqIIYzSZKkDmI4kyRJ6iB71ZwzSZI0Mmzfvp1Vq1axbdu2dpdSuQkTJjB79mzGjRvX0PGGM0mSNOxWrVrF1KlTmTdvHhHR7nIqk5msW7eOVatWcdhhhzV0jsOakiRp2G3bto3p06fv1cEMICKYPn16Uz2EhjNJktQWe3swq2n27zScSZKkUWf9+vVceumlTZ93+umns379+goq6mE4kyRJo05/4Wznzp0DnnfzzTczbdq0qsoCvCGgKTfcAFOmwCmntLsSSZI0FBdeeCE///nPmT9/PuPGjWPKlCkcfPDBLF26lEceeYQzzzyTlStXsm3bNs4//3wWLVoE9LyNaPPmzZx22mm86U1v4q677mLWrFl861vfYuLEiUOuzZ6zJlx0EVxySburkCRJQ/XpT3+aww8/nKVLl/KZz3yGe++9l0996lM88sgjAFxxxRUsWbKExYsXc/HFF7Nu3brdvuNnP/sZ5557Lg8//DDTpk3j+uuvb0lt9pw1obsbXnyx3VVIkrR3ueACWLq0td85fz58/vONH79w4cKXPeri4osv5sYbbwRg5cqV/OxnP2P69OkvO+ewww5j/vz5ALz+9a9nxYoVQ64bDGdNGT8eXnqp3VVIkqRWmzx58q+2f/jDH/L973+ff//3f2fSpEmcdNJJfT4Ko7u7+1fbXV1dbN26tSW1GM6a0N0NLfrPXZIklZrp4WqVqVOnsmnTpj73bdiwgf32249Jkybx6KOPcvfddw9rbZWFs4iYANwBdJe/843M/ESvY94N/FX5cTPwx5n503LfCmATsBPYkZkLqqq1Ud3dUPHds5IkaRhMnz6dE044gaOPPpqJEydy4IEH/mrfqaeeymWXXcYxxxzDr/3ar3H88ccPa21V9py9CLw5MzdHxDjgzoi4JTPr4+cTwImZ+XxEnAZcDryhbv/JmflshTU2xTlnkiTtPa655po+27u7u7nlllv63FebVzZjxgweeuihX7X/+Z//ecvqqiycZWZS9IYBjCuX7HXMXXUf7wZmV1VPKzjnTJIkVa3SR2lERFdELAXWALdl5j0DHP4BoD6mJnBrRCyJiEUD/MaiiFgcEYvXrl3bmsL7Yc+ZJEmqWqXhLDN3ZuZ8ih6xhRFxdF/HRcTJFOHsr+qaT8jM1wGnAedGxG/28xuXZ+aCzFwwc+bMFv8FL2c4kyRJVRuWh9Bm5nrgh8CpvfdFxDHAl4EzMnNd3TlPl+s1wI3AwuGodSAOa0qS1DrFDKi9X7N/Z2XhLCJmRsS0cnsi8FvAo72OmQvcAJyTmf9R1z45IqbWtoFTgIdoM3vOJElqjQkTJrBu3bq9PqBlJuvWrWPChAkNn1Pl3ZoHA1dGRBdFCPx6Zn47Iv4IIDMvAz4OTAcujQjoeWTGgcCNZdtY4JrM/G6FtTbEcCZJUmvMnj2bVatWUfV88U4wYcIEZs9u/J7HKu/WfAA4ro/2y+q2Pwh8sI9jHgeOraq2weruhh07YNcuGONbSSVJGrRx48a97HVJ6mHEaML48cXaeWeSJKkqhrMm1F6h5dCmJEmqiuGsCYYzSZJUNcNZExzWlCRJVTOcNcGeM0mSVDXDWRMMZ5IkqWqGsyYYziRJUtUMZ01wzpkkSaqa4awJ9pxJkqSqGc6aYDiTJElVM5w1wWFNSZJUNcNZE+w5kyRJVTOcNcFwJkmSqmY4a4LhTJIkVc1w1gTnnEmSpKoZzppgz5kkSaqa4awJhjNJklQ1w1kTHNaUJElVM5w1oRbO7DmTJElVMZw1IaIIaIYzSZJUFcNZk7q7DWeSJKk6hrMmjR/vnDNJklQdw1mT7DmTJElVMpw1yXAmSZKqZDhrksOakiSpSoazJtlzJkmSqlRZOIuICRFxb0T8NCIejohP9nFMRMTFEbE8Ih6IiNfV7Ts1Ih4r911YVZ3NMpxJkqQqVdlz9iLw5sw8FpgPnBoRx/c65jTgiHJZBHwBICK6gEvK/UcBZ0fEURXW2jDDmSRJqlJl4SwLm8uP48olex12BnBVeezdwLSIOBhYCCzPzMcz8yXguvLYtnPOmSRJqlKlc84ioisilgJrgNsy855eh8wCVtZ9XlW29dfedvacSZKkKlUazjJzZ2bOB2YDCyPi6F6HRF+nDdC+m4hYFBGLI2Lx2rVrh1ZwAwxnkiSpSsNyt2Zmrgd+CJzaa9cqYE7d59nA0wO09/Xdl2fmgsxcMHPmzJbV3B+HNSVJUpWqvFtzZkRMK7cnAr8FPNrrsJuA95Z3bR4PbMjM1cB9wBERcVhEjAfOKo9tO3vOJElSlcZW+N0HA1eWd16OAb6emd+OiD8CyMzLgJuB04HlwAvA+8t9OyLiPOB7QBdwRWY+XGGtDTOcSZKkKlUWzjLzAeC4Ptovq9tO4Nx+zr+ZIrx1FMOZJEmqkm8IaJJzziRJUpUMZ02y50ySJFXJcNak7m7Yvh127Wp3JZIkaW9kOGvS+PHFevv29tYhSZL2ToazJnV3F2uHNiVJUhUMZ00ynEmSpCoZzppkOJMkSVUynDWpNufMx2lIkqQqGM6aZM+ZJEmqkuGsSYYzSZJUJcNZkxzWlCRJVTKcNcmeM0mSVCXDWZMMZ5IkqUqGsyYZziRJUpUMZ01yzpkkSaqS4axJ9pxJkqQqGc6aZDiTJElVMpw1yWFNSZJUJcNZk+w5kyRJVTKcNclwJkmSqmQ4a1JtWNNwJkmSqmA4a5JzziRJUpUMZ00aMwbGjbPnTJIkVcNwNgjd3YYzSZJUDcPZIIwf77CmJEmqhuFsEOw5kyRJVRlb1RdHxBzgKuAgYBdweWb+Y69j/gJ4d10tRwIzM/O5iFgBbAJ2Ajsyc0FVtTbLcCZJkqoyYDiLiABmZ+bKQXz3DuDDmXl/REwFlkTEbZn5SO2AzPwM8Jnyt94O/GlmPlf3HSdn5rOD+O1KOawpSZKqMuCwZmYm8M3BfHFmrs7M+8vtTcAyYNYAp5wNXDuY3xpu9pxJkqSqNDLn7O6I+PWh/EhEzAOOA+7pZ/8k4FTg+rrmBG6NiCURsWgov99qhjNJklSVRuacnQz8YUQ8CWwBgqJT7ZhGfiAiplCErgsyc2M/h70d+LdeQ5onZObTEXEAcFtEPJqZd/Tx/YuARQBz585tpKQhM5xJkqSqNBLOThvsl0fEOIpgdnVm3jDAoWfRa0gzM58u12si4kZgIbBbOMvMy4HLARYsWJCDrbUZzjmTJElV2eOwZmY+CUyj6N16OzCtbBtQeTPBV4Blmfm5AY7bFzgR+FZd2+TyJgIiYjJwCvDQnn5zuNhzJkmSqrLHnrOIOB/4A6DW8/W1iLg8M//3Hk49ATgHeDAilpZtHwXmAmTmZWXbfwFuzcwtdeceCNxY5DvGAtdk5ncb+HuGheFMkiRVpZFhzQ8Ab6iFp4j4X8C/AwOGs8y8k2J+2oAy81+Af+nV9jhwbAO1tYXDmpIkqSqN3K0ZFA+CrdlJA6Frb2bPmSRJqkojPWdXAPeUk/IBzqSYSzZqGc4kSVJV9vSGgDEUzyb7EfAmih6z92fmT4ahto5lOJMkSVUZMJxl5q6I+GxmvhG4f5hq6njOOZMkSVVpZM7ZrRHxu+WjMYQ9Z5IkqTqNzDn7M2AysCMittHzhoB9Kq2sg3V3Fz1nmWBklSRJrTRgz1k55+zUzByTmeMzc5/MnDqagxkUw5oA27e3tw5JkrT3GTCcZeYu4O+HqZYRo7u7WDu0KUmSWs05Z4NgOJMkSVVpZs7ZzojYinPODGeSJKkyewxnmTl1OAoZSWpzznychiRJarU9DmtG4T0R8bHy85yIWFh9aZ3LnjNJklSVRuacXQq8EXhX+XkzcEllFY0AhjNJklSVRuacvSEzXxcRPwHIzOcjYnzFdXU0hzUlSVJVGuk52x4RXUACRMRMYFelVXU4e84kSVJVGglnFwM3AgdExKeAO4G/qbSqDmc4kyRJVWnkbs2rI2IJ8BaKx2icmZnLKq+sgxnOJElSVRqZc0ZmPgo8WnEtI4ZzziRJUlUaGdZUL/acSZKkqhjOBsFwJkmSqmI4GwSHNSVJUlX6nXMWEZsoH5/RF9+tac+ZJElqvX7DWe2dmhFxEfBL4KsUd2u+GxjV79s0nEmSpKo0Mqz525l5aWZuysyNmfkF4HerLqyTGc4kSVJVGglnOyPi3RHRFRFjIuLdwM6qC+tkzjmTJElVaSScvQv4PeCZcnknPS9BH5XGjIGxY+05kyRJrdfIGwJWAGdUX8rI0t1tOJMkSa23x56ziHhVRPwgIh4qPx8TEX/dwHlzIuL2iFgWEQ9HxPl9HHNSRGyIiKXl8vG6fadGxGMRsTwiLmz2D6va+PEOa0qSpNZrZFjzS8BHgO0AmfkAcFYD5+0APpyZRwLHA+dGxFF9HPfjzJxfLhcBREQXcAlwGnAUcHY/57aNPWeSJKkKjYSzSZl5b6+2HXs6KTNXZ+b95fYmYBkwq8G6FgLLM/PxzHwJuI4OG1o1nEmSpCo0Es6ejYjDKR9IGxHvAFY38yMRMQ84Drinj91vjIifRsQtEfGasm0WsLLumFX0E+wiYlFELI6IxWvXrm2mrCExnEmSpCrs8YYA4FzgcuDVEfEU8ATFg2gbEhFTgOuBCzJzY6/d9wOHZubmiDgd+CZwBMXDbnvr820FmXl5WR8LFizo940GreacM0mSVIUBw1k59+uPM/O3ImIyMKYcomxIRIyjCGZXZ+YNvffXh7XMvDkiLo2IGRQ9ZXPqDp0NPN3o7w4He84kSVIVBhzWzMydwOvL7S1NBrMAvgIsy8zP9XPMQeVxRMTCsp51wH3AERFxWESMp7gB4aZGf3s4GM4kSVIVGhnW/ElE3AT8K7Cl1thXT1gvJwDnAA9GxNKy7aPA3PL8y4B3AH8cETuArcBZmZnAjog4D/ge0AVckZkPN/5nVc9hTUmSVIVGwtn+FL1Zb65rS2DAcJaZd9L33LH6Y/4J+Kd+9t0M3NxAfW3R3Q3PP9/uKiRJ0t6mkTcEvH84ChlpHNaUJElV2GM4i4gJwAeA1wATau2Z+d8qrKvjGc4kSVIVGnnO2VeBg4DfBn5EcedkwzcG7K2ccyZJkqrQSDh7ZWZ+DNiSmVcCbwVeW21Znc+eM0mSVIVGwtn2cr0+Io4G9gXmVVbRCGE4kyRJVWjkbs3LI2I/4GMUzxqbAny80qpGAIc1JUlSFRq5W/PL5eaPgFdUW87IYc+ZJEmqQiN3a/bZS5aZF7W+nJGjFs4yIQZ8mpskSVLjGplztqVu2QmchnPO6O4u1tu3D3ycJElSMxoZ1vxs/eeI+Hs67D2X7TB+fLF+6aWebUmSpKFqpOest0k49+xXPWfOO5MkSa3UyJyzBynepQnFS8hnAqN6vhkYziRJUjUaeZTG2+q2dwDPZOaOiuoZMeqHNSVJklqlkXDW+1VN+0Td7YmZ+VxLKxoh7DmTJElVaCSc3Q/MAZ4HApgG/KLcl4zS+WeGM0mSVIVGbgj4LvD2zJyRmdMphjlvyMzDMnNUBjMwnEmSpGo0Es5+PTNvrn3IzFuAE6sraWRwzpkkSapCI8Oaz0bEXwNfoxjGfA+wrtKqRgB7ziRJUhUa6Tk7m+LxGTcC3yy3z66yqJHAcCZJkqrQyBsCngPOB4iILmByZm6surBO57CmJEmqwh57ziLimojYJyImAw8Dj0XEX1RfWmez50ySJFWhkWHNo8qesjOBm4G5wDmVVjUCGM4kSVIVGgln4yJiHEU4+1ZmbqfndU6jluFMkiRVoZFw9kVgBTAZuCMiDgWcc+acM0mSVIE9hrPMvDgzZ2Xm6ZmZFG8HOLn60jqbPWeSJKkKjTzn7GXKgDbqX3xuOJMkSVVoZFhzUCJiTkTcHhHLIuLhiDi/j2PeHREPlMtdEXFs3b4VEfFgRCyNiMVV1TlYDmtKkqQqNN1z1oQdwIcz8/6ImAosiYjbMvORumOeAE7MzOcj4jTgcuANdftPzsxnK6xx0Lq6isWeM0mS1EoNhbOI+A1gXv3xmXnVQOdk5mpgdbm9KSKWAbOAR+qOuavulLuB2Y0W3gm6uw1nkiSptfYYziLiq8DhwFJgZ9mcwIDhrNd3zAOOA+4Z4LAPALfUfU7g1ohI4IuZeXmjvzdcDGeSJKnVGuk5W0DxINpBPdssIqYA1wMX9Pfap4g4mSKcvamu+YTMfDoiDgBui4hHM/OOPs5dBCwCmDt37mBKHLTx451zJkmSWquRGwIeAg4azJeXD6+9Hrg6M2/o55hjgC8DZ2Tmulp7Zj5drtdQvHR9YV/nZ+blmbkgMxfMnDlzMGUOmj1nkiSp1RrpOZsBPBIR9wK/iiKZ+TsDnRQRAXwFWJaZn+vnmLnADcA5mfkfde2TgTHlXLXJwCnARQ3UOqwMZ5IkqdUaCWf/Y5DffQLFOzgfjIilZdtHKd7NSWZeBnwcmA5cWmQ5dmTmAuBA4MaybSxwTWZ+d5B1VMZhTUmS1Gp7DGeZ+aPBfHFm3gnEHo75IPDBPtofB47d/YzOYs+ZJElqtT3OOYuI4yPivojYHBEvRcTOiBj179YEw5kkSWq9Rm4I+CfgbOBnwESKnq5/qrKokcJwJkmSWq2h1zdl5nKgKzN3ZuY/AydVWtUI4ZwzSZLUao3cEPBCRIwHlkbE31E89X9ytWWNDPacSZKkVmuk5+yc8rjzgC3AHOB3qyxqpDCcSZKkVmvkbs0nI2IicHBmfnIYahoxHNaUJEmt1sjdmm+neK/md8vP8yPipqoLGwnsOZMkSa3WyLDm/6B4ddJ6gMxcCsyrrqSRw3AmSZJarZFwtiMzN1ReyQhkOJMkSa3WyN2aD0XEu4CuiDgC+BBwV7VljQzOOZMkSa3WSM/ZnwCvoXjp+bXARuCCKosaKWo9Z5ntrkSSJO0tGrlb8wXgv5eL6nR3F8Fsxw4YN67d1UiSpL1Bv+FsT3dkZubvtL6ckWX8+GL90kuGM0mS1BoD9Zy9EVhJMZR5DxDDUtEI0t1drF98ESb7zgRJktQCA4Wzg4D/TPHS83cB3wGuzcyHh6OwkaA+nEmSJLVCvzcElC85/25mvg84HlgO/DAi/mTYqutwhjNJktRqA94QEBHdwFspes/mARcDN1Rf1shQP+dMkiSpFQa6IeBK4GjgFuCTmfnQsFU1QthzJkmSWm2gnrNzgC3Aq4APRfzqfoAAMjP3qbi2jmc4kyRJrdZvOMvMRh5QO6o5rClJklrNADYE9pxJkqRWM5wNgeFMkiS1muFsCGrDmoYzSZLUKoazIaj1nDnnTJIktYrhbAgc1pQkSa1mOBsCw5kkSWq1ysJZRMyJiNsjYllEPBwR5/dxTETExRGxPCIeiIjX1e07NSIeK/ddWFWdQ+GjNCRJUqtV2XO2A/hwZh5J8W7OcyPiqF7HnAYcUS6LgC8AREQXcEm5/yjg7D7ObTt7ziRJUqtVFs4yc3Vm3l9ubwKWAbN6HXYGcFUW7gamRcTBwEJgeWY+npkvAdeVx3YUw5kkSWq1YZlzFhHzgOOAe3rtmgWsrPu8qmzrr72j+CgNSZLUapWHs4iYAlwPXJCZG3vv7uOUHKC9r+9fFBGLI2Lx2rVrh1Zsk8aOhTFjnHMmSZJap9JwFhHjKILZ1Zl5Qx+HrALm1H2eDTw9QPtuMvPyzFyQmQtmzpzZmsKb0N1tz5kkSWqdKu/WDOArwLLM/Fw/h90EvLe8a/N4YENmrgbuA46IiMMiYjxwVnlsxzGcSZKkVhpb4XefAJwDPBgRS8u2jwJzATLzMuBm4HRgOfAC8P5y346IOA/4HtAFXJGZD1dY66CNH++wpiRJap3Kwllm3knfc8fqj0ng3H723UwR3jqaPWeSJKmVfEPAEBnOJElSKxnOhmj8eMOZJElqHcPZEHV3O+dMkiS1juFsiBzWlCRJrWQ4GyLDmSRJaiXD2RD5KA1JktRKhrMhmjAB1q9vdxWSJGlvYTgbore8BR59FO66q92VSJKkvYHhbIg++EGYPh3+9m/bXYkkSdobGM6GaPJk+NCH4NvfhgcfbHc1kiRppDOctcB558GUKfDpT7e7EkmSNNIZzlpg//3hD/8QrrsOHn+83dVIkqSRzHDWIn/2ZzB2LHzmM+2uRJIkjWSGsxY55BB43/vgn/8ZfvnLdlcjSZJGKsNZC/3lX8L27fAP/9DuSiRJ0khlOGuhV74S3vlO+MIXfDCtJEkaHMNZi114IWzaBJde2u5KJEnSSGQ4a7H58+G00+Dv/g4WL253NZIkaaQxnFXg0kthv/2KVzv927+1uxpJkjSSGM4qMG8e/PjHcNBBcMop8IMftLsiSZI0UhjOKjJ7NtxxB7ziFfDWt8J3vtPuiiRJ0khgOKvQgQfCD38Ir30tnHkm/Ou/trsiSZLU6QxnFZs+Hb7/fXjDG+D3fg/OPhtWrGh3VZIkqVMZzobBvvvCrbfCxz4G3/oWvPrV8JGPwMaN7a5MkiR1GsPZMJk0CS66CB57rOhB+/Sni4fWXnopbN3a7uokSVKnMJwNszlz4Kqr4L774Mgj4dxz4eCD4bzz4Kc/bXd1kiSp3QxnbbJgQXGzwO23F3dzfvnLxQNsFy6EL30J1q1rd4WSJKkdKgtnEXFFRKyJiIf62f8XEbG0XB6KiJ0RsX+5b0VEPFju22ufsx8BJ50EV18NTz8Nn/88vPACLFoEBxwAJ54In/scLF/e7kolSdJwicys5osjfhPYDFyVmUfv4di3A3+amW8uP68AFmTms8385oIFC3LxCH9nUiYsWVLcOHDTTfDAA0X7kUcWD7Q9+eQitE2b1t46JUnS0ETEksxcsFt7VeGs/NF5wLcbCGfXALdn5pfKzysYpeGstyeegP/7f4vlzjth2zYYMwaOO64IaiecUAyFHnJIuyuVJEnN6NhwFhGTgFXAKzPzubLtCeB5IIEvZublA5y/CFgEMHfu3Nc/+eSTLau/07z4Itx9dzFP7fbbi+2XXir2zZpVhLRf//ViPtv8+TBzZnvrlSRJ/evkcPb7wHsy8+11bYdk5tMRcQBwG/AnmXnHnn5vb+w5G8jWrbB0Kdx7b89SPz/tkEOKkHbccXDMMcXQ6KteBd3d7atZkiQV+gtnY9tRTC9nAdfWN2Tm0+V6TUTcCCwE9hjORpuJE+GNbyyWmueeKwJb/fK978HOncX+ri44/PAiqNXCWm2ZMaO4SUGSJLVPW8NZROwLnAi8p65tMjAmMzeV26cAF7WpxBFn//3hzW8ulppt24qH3y5bBo880rN85zuwY0fPcdOmwRFHFOGtfnnFK4peuDE+eEWSpMpVFs4i4lrgJGBGRKwCPgGMA8jMy8rD/gtwa2ZuqTv1QODGKLpwxgLXZOZ3q6pzNJgwAY49tljq7dhRvOfzP/7j5cs99xQvaa/1tgGMGwdz58K8eXDoocV67tziobpz58Ls2cXvSJKkoal0ztlwG21zzqq0fTv84hfw858Xy5NPFsuKFcV69erdzznggCKk9V5mzSp63mbNgqlTh/1PkSSpI3XynDN1oHHjeoY1+7JtGzz1VBHgfvELWLmyWDtMEdYAAAzFSURBVK9aVTz+48c/huef3/28KVOKoHbIIcVrq2pL7fNBBxXrffd1/pskaXQynGlQJkwYOLxB8baDVauKtx889VSxrm2vXl0Mn65e3feL3ydMKILaQQfBgQcWvXIHHvjy7dp6v/0McpKkvYfhTJWZNKnnTtD+ZMLGjUVIW70afvnL3ddPPFE8023tWti1a/fvGDu2CGp9LTNnFnehzpzZs22vnCSpkxnO1FYRRVjad1949asHPnbnzuJRIc88Uyxr1vRsP/NMEd7WrCluanjmmb575KAYsq0Ftvp1X8v06cV64sTW/+2SJPXFcKYRo6urpwfs6AFfCFbYsqUIbGvXwrPP9r99//3Fuq85cjUTJ/aEtVpgq23XL/vv37OeNs3Hj0iSmmc4015r8uRimTevseN37Ch65p59tljWrXv5ura9bl0R6NatKwJdfzc8RxTz4eoDW22p/7zffj3H1bbHjWvZfwySpBHGcCaV6ueuNWrnTli/vghqzz23+7p+e82a4mHA69bBhg0Df++UKT1hrb91bZk2rWe9774GO0ka6Qxn0hB0dfUMaTZjx46eUPf88z3Lc8+9fLv2+dFHe9pefHHg7548ue/wVtuuD3O9lylTHIqVpHYznEltMHZsz00Hzdq6tSeorV/fs9RCXf32888Xd7vef3/RW7dp08DfHQH77NPTC1e/rl/q9/Xe7u4e3H8mkqSC4UwaYSZOLJZDDmn+3B07ipBWC3AbNvR87qttw4bijRAPPFB83rix/zl2Nd3dPUFt332LsFf/uXd7/f599imWKVN83Imk0ctwJo0iY8cObhi2ZteuovetPsT1DnO9l40bi0eb1LY3btzz79R68AZa6sPcPvsUrwbrfczUqcXQsySNJIYzSQ0bM6anl2uwdu2CzZtfHt7qt+uX2lDsxo1FIHzyyZ62zZsb+73Jk3cPb7XtqVN33x5o8WYLScPBcCZpWI0Z0xOS5swZ/Pfs3FkEtN4hrq9gV79v06ZiHl79vpdeauw3u7v7Dm1Tpgy8Xb+u3+7udvhW0u4MZ5JGpK6unl68oYQ8KO6A3bRpcMu6dUWPXq03b9OmIjg2+jcMFOJ6L33tmzz55Z8nTDDwSSOd4UzSqNfdXSyDuXu2t0zYtq0nqNUvW7b0hLja/vp1bXvdup7PmzcX5zVqzJieBzDXh7f+1r1DXu3cSZN2/x7n70nDw3AmSS0U0XNH7cyZrfnOXbvghRdeHthqYa8W3nqHud7bGzbAU0+9vG3btubq6O5+eYDra6kPdb0DX/3niROLY2vrCRN8xp5UYziTpA43ZkxPD1cr7dy5e8DbsqX/pfdxL7xQrJ95Zvdj9/Sw5L5MnDhw4OtrXduuhbz6wNd72wCokcJwJkmjVFdXz80ZrbZjR09462vZurXYX1v3d+wLL+we/rZuLZbBqPVq1gJb/VIf+OqD3UDr3m217bH+21VD4H99JEktN3ZsdcEPiqHeWrCrD3v1ga++rfe++kC4dWuxXru2JxDWjh1MDyAUf/9A4a1+e8KEl6/7C4G9j50woWe/j3nZuxjOJEkjTv2ND62a29eXnTuLuXn1ga13+Ktv6+uY3tsbNsDq1T1t27YNrTcQil7QgYJdLczVL70DYX/b9efX7+vudpi4KoYzSZL60dXVEwKrllk8c6+vnr9agNu27eVhrvfSVyDcvLm4A7h2bu/vGorx43cPbgNtD3bp7t59e29+TqDhTJKkDhDR81iXadOG5zczi6HbWqDrHd7q23uve4e92vfUH7NuXdHe17GNPvx5IPWhbaBA11ewG+i4iRPhrW8den2DZTiTJGmUiugJJPvtN7y/vWvX7oGtrxC3bVvPek9tvc/dsKG4oaQWHHufm9l3bZMmNfd8wVYznEmSpGE3ZkzPjRHtkFncVdxXaNu+vT011RjOJEnSqBNR3OXaiXe6ep+FJElSB6ksnEXEFRGxJiIe6mf/SRGxISKWlsvH6/adGhGPRcTyiLiwqholSZI6TZU9Z/8CnLqHY36cmfPL5SKAiOgCLgFOA44Czo6IoyqsU5IkqWNUFs4y8w7guUGcuhBYnpmPZ+ZLwHXAGS0tTpIkqUO1e87ZGyPipxFxS0S8pmybBaysO2ZV2SZJkrTXa+fdmvcDh2bm5og4HfgmcATQ1/N++3kSCUTEImARwNy5c6uoU5Ikadi0recsMzdm5uZy+2ZgXETMoOgpm1N36Gzg6QG+5/LMXJCZC2ZW+YI1SZKkYdC2cBYRB0UUb8WKiIVlLeuA+4AjIuKwiBgPnAXc1K46JUmShlNlw5oRcS1wEjAjIlYBnwDGAWTmZcA7gD+OiB3AVuCszExgR0ScB3wP6AKuyMyHq6pTkiSpk0T292KpEWjBggW5ePHidpchSZK0RxGxJDMX7Na+N4WziFgLPFnxz8wAnq34N9Q8r0vn8tp0Jq9L5/LadKYqrsuhmbnbhPm9KpwNh4hY3FfKVXt5XTqX16YzeV06l9emMw3ndWn3c84kSZJUx3AmSZLUQQxnzbu83QWoT16XzuW16Uxel87ltelMw3ZdnHMmSZLUQew5kyRJ6iCGswZFxKkR8VhELI+IC9tdz2gWEXMi4vaIWBYRD0fE+WX7/hFxW0T8rFzv1+5aR6OI6IqIn0TEt8vPXpcOEBHTIuIbEfFo+b+dN3pt2i8i/rT859hDEXFtREzwurRHRFwREWsi4qG6tn6vRUR8pMwEj0XEb7eyFsNZAyKiC7gEOA04Cjg7Io5qb1Wj2g7gw5l5JHA8cG55PS4EfpCZRwA/KD9r+J0PLKv77HXpDP8IfDczXw0cS3GNvDZtFBGzgA8BCzLzaIq34pyF16Vd/gU4tVdbn9ei/HfOWcBrynMuLbNCSxjOGrMQWJ6Zj2fmS8B1wBltrmnUyszVmXl/ub2J4l8ysyiuyZXlYVcCZ7anwtErImYDbwW+XNfsdWmziNgH+E3gKwCZ+VJmrsdr0wnGAhMjYiwwCXgar0tbZOYdwHO9mvu7FmcA12Xmi5n5BLCcIiu0hOGsMbOAlXWfV5VtarOImAccB9wDHJiZq6EIcMAB7ats1Po88JfArro2r0v7vQJYC/xzOeT85YiYjNemrTLzKeDvgV8Aq4ENmXkrXpdO0t+1qDQXGM4aE320eZtrm0XEFOB64ILM3Njueka7iHgbsCYzl7S7Fu1mLPA64AuZeRywBYfK2q6cv3QGcBhwCDA5It7T3qrUoEpzgeGsMauAOXWfZ1N0PatNImIcRTC7OjNvKJufiYiDy/0HA2vaVd8odQLwOxGxgmLo/80R8TW8Lp1gFbAqM+8pP3+DIqx5bdrrt4AnMnNtZm4HbgB+A69LJ+nvWlSaCwxnjbkPOCIiDouI8RSTAG9qc02jVkQExdyZZZn5ubpdNwHvK7ffB3xruGsbzTLzI5k5OzPnUfxv5P9l5nvwurRdZv4SWBkRv1Y2vQV4BK9Nu/0COD4iJpX/XHsLxRxar0vn6O9a3AScFRHdEXEYcARwb6t+1IfQNigiTqeYT9MFXJGZn2pzSaNWRLwJ+DHwID1zmz5KMe/s68Bcin/ovTMze0/u1DCIiJOAP8/Mt0XEdLwubRcR8ylu1BgPPA68n+L/oHtt2igiPgn8PsVd6D8BPghMwesy7CLiWuAkYAbwDPAJ4Jv0cy0i4r8D/43i2l2Qmbe0rBbDmSRJUudwWFOSJKmDGM4kSZI6iOFMkiSpgxjOJEmSOojhTJIkqYMYziRpiCLipIj4drvrkLR3MJxJkiR1EMOZpFEjIt4TEfdGxNKI+GJEdEXE5oj4bETcHxE/iIiZ5bHzI+LuiHggIm4s34NIRLwyIr4fET8tzzm8/PopEfGNiHg0Iq4un/guSU0znEkaFSLiSIonsZ+QmfOBncC7gcnA/Zn5OuBHFE8FB7gK+KvMPIbibRS19quBSzLzWIr3IK4u248DLgCOAl5B8a5RSWra2HYXIEnD5C3A64H7yk6tiRQvMd4F/J/ymK8BN0TEvsC0zPxR2X4l8K8RMRWYlZk3AmTmNoDy++7NzFXl56XAPODO6v8sSXsbw5mk0SKAKzPzIy9rjPhYr+MGeqfdQEOVL9Zt78R/vkoaJIc1JY0WPwDeEREHAETE/hFxKMU/B99RHvMu4M7M3AA8HxH/qWw/B/hRZm4EVkXEmeV3dEfEpGH9KyTt9fx/dpJGhcx8JCL+Grg1IsYA24FzgS3AayJiCbCBYl4awPuAy8rw9Tjw/rL9HOCLEXFR+R3vHMY/Q9IoEJkD9eBL0t4tIjZn5pR21yFJNQ5rSpIkdRB7ziRJkjqIPWeSJEkdxHAmSZLUQQxnkiRJHcRwJkmS1EEMZ5IkSR3EcCZJktRB/j9ZYQnIuCxaxwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mod.history_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bây giờ ta có thể dùng method Summary để lấy nội dung betas và etas:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"model_1\"\n",
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "input_1 (InputLayer)            (None, 9)            0                                            \n",
      "__________________________________________________________________________________________________\n",
      "input_2 (InputLayer)            (None, 4112)         0                                            \n",
      "__________________________________________________________________________________________________\n",
      "Fixed_effect (Dense)            (None, 1)            10          input_1[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "Random_Effect (Dense)           (None, 1)            4112        input_2[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "add_1 (Add)                     (None, 1)            0           Fixed_effect[0][0]               \n",
      "                                                                 Random_Effect[0][0]              \n",
      "==================================================================================================\n",
      "Total params: 4,122\n",
      "Trainable params: 4,122\n",
      "Non-trainable params: 0\n",
      "__________________________________________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "fe_df, re_df = mod.Summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nội dung betas như sau:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "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>Keras_GLMM</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td>3.186542</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Stud_Age[T.4]</th>\n",
       "      <td>0.053660</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Stud_Age[T.6]</th>\n",
       "      <td>0.062380</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Stud_Age[T.8]</th>\n",
       "      <td>0.141339</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.2]</th>\n",
       "      <td>-0.067621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.3]</th>\n",
       "      <td>-0.107238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.4]</th>\n",
       "      <td>-0.167474</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.5]</th>\n",
       "      <td>-0.154622</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.6]</th>\n",
       "      <td>-0.230838</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Service[T.1]</th>\n",
       "      <td>-0.094988</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               Keras_GLMM\n",
       "Intercept        3.186542\n",
       "Stud_Age[T.4]    0.053660\n",
       "Stud_Age[T.6]    0.062380\n",
       "Stud_Age[T.8]    0.141339\n",
       "Lect_Age[T.2]   -0.067621\n",
       "Lect_Age[T.3]   -0.107238\n",
       "Lect_Age[T.4]   -0.167474\n",
       "Lect_Age[T.5]   -0.154622\n",
       "Lect_Age[T.6]   -0.230838\n",
       "Service[T.1]    -0.094988"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fe_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ta thấy rằng kết quả này rất gần với fixed effect của R với package lme4\n",
    "\n",
    "Còn đây là 4112 giá trị của random effects (tương ứng với mỗi hàng trong matrix Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0       0.023151\n",
       "1      -0.015272\n",
       "2       0.017752\n",
       "3       0.011739\n",
       "4      -0.008241\n",
       "          ...   \n",
       "4107   -0.113134\n",
       "4108    0.034460\n",
       "4109    0.157137\n",
       "4110    0.017055\n",
       "4111    0.145813\n",
       "Length: 4112, dtype: float32"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "re_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trong python, ta có thể dựng GLMM bằng package statsmodel, tuy nhiên có giới hạn là chỉ chấp nhận 1 random effect duy nhất, ta không thể code các loại cross, nested random effects cho nhiều biến"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "       <td>Model:</td>       <td>MixedLM</td> <td>Dependent Variable:</td>     <td>Score</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <td>No. Observations:</td>  <td>73421</td>        <td>Method:</td>           <td>REML</td>    \n",
       "</tr>\n",
       "<tr>\n",
       "     <td>No. Groups:</td>     <td>2972</td>         <td>Scale:</td>           <td>1.6654</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <td>Min. group size:</td>     <td>1</td>      <td>Log-Likelihood:</td>   <td>-124203.2237</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <td>Max. group size:</td>    <td>92</td>        <td>Converged:</td>           <td>Yes</td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <td>Mean group size:</td>   <td>24.7</td>            <td></td>                 <td></td>      \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "        <td></td>         <th>Coef.</th> <th>Std.Err.</th>    <th>z</th>    <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>      <td>3.283</td>   <td>0.016</td>  <td>198.950</td> <td>0.000</td>  <td>3.250</td>  <td>3.315</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Stud_Age[T.4]</th>  <td>0.074</td>   <td>0.022</td>   <td>3.288</td>  <td>0.001</td>  <td>0.030</td>  <td>0.118</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Stud_Age[T.6]</th>  <td>0.093</td>   <td>0.022</td>   <td>4.179</td>  <td>0.000</td>  <td>0.050</td>  <td>0.137</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Stud_Age[T.8]</th>  <td>0.175</td>   <td>0.024</td>   <td>7.340</td>  <td>0.000</td>  <td>0.128</td>  <td>0.222</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Lect_Age[T.2]</th> <td>-0.087</td>   <td>0.014</td>  <td>-6.214</td>  <td>0.000</td> <td>-0.114</td> <td>-0.059</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Lect_Age[T.3]</th> <td>-0.150</td>   <td>0.017</td>  <td>-9.009</td>  <td>0.000</td> <td>-0.182</td> <td>-0.117</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Lect_Age[T.4]</th> <td>-0.198</td>   <td>0.018</td>  <td>-11.065</td> <td>0.000</td> <td>-0.233</td> <td>-0.163</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Lect_Age[T.5]</th> <td>-0.188</td>   <td>0.021</td>  <td>-9.159</td>  <td>0.000</td> <td>-0.229</td> <td>-0.148</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Lect_Age[T.6]</th> <td>-0.288</td>   <td>0.018</td>  <td>-15.992</td> <td>0.000</td> <td>-0.323</td> <td>-0.253</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Service[T.1]</th>  <td>-0.092</td>   <td>0.010</td>  <td>-9.017</td>  <td>0.000</td> <td>-0.112</td> <td>-0.072</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Group Var</th>      <td>0.101</td>   <td>0.004</td>     <td></td>       <td></td>       <td></td>       <td></td>   \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary2.Summary'>\n",
       "\"\"\"\n",
       "          Mixed Linear Model Regression Results\n",
       "==========================================================\n",
       "Model:            MixedLM Dependent Variable: Score       \n",
       "No. Observations: 73421   Method:             REML        \n",
       "No. Groups:       2972    Scale:              1.6654      \n",
       "Min. group size:  1       Log-Likelihood:     -124203.2237\n",
       "Max. group size:  92      Converged:          Yes         \n",
       "Mean group size:  24.7                                    \n",
       "----------------------------------------------------------\n",
       "               Coef.  Std.Err.    z    P>|z| [0.025 0.975]\n",
       "----------------------------------------------------------\n",
       "Intercept       3.283    0.016 198.950 0.000  3.250  3.315\n",
       "Stud_Age[T.4]   0.074    0.022   3.288 0.001  0.030  0.118\n",
       "Stud_Age[T.6]   0.093    0.022   4.179 0.000  0.050  0.137\n",
       "Stud_Age[T.8]   0.175    0.024   7.340 0.000  0.128  0.222\n",
       "Lect_Age[T.2]  -0.087    0.014  -6.214 0.000 -0.114 -0.059\n",
       "Lect_Age[T.3]  -0.150    0.017  -9.009 0.000 -0.182 -0.117\n",
       "Lect_Age[T.4]  -0.198    0.018 -11.065 0.000 -0.233 -0.163\n",
       "Lect_Age[T.5]  -0.188    0.021  -9.159 0.000 -0.229 -0.148\n",
       "Lect_Age[T.6]  -0.288    0.018 -15.992 0.000 -0.323 -0.253\n",
       "Service[T.1]   -0.092    0.010  -9.017 0.000 -0.112 -0.072\n",
       "Group Var       0.101    0.004                            \n",
       "==========================================================\n",
       "\n",
       "\"\"\""
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = \"Score ~ Stud_Age + Lect_Age + Service\"\n",
    "\n",
    "md = smf.mixedlm(formula, df, groups=df[\"Student\"])\n",
    "\n",
    "mdf = md.fit()\n",
    "\n",
    "mdf.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ta có thể so sánh nội dung fixed effect betas của 2 mô hình Neural network Keras và GLMM của statsmodel, chúng tương tự nhau:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "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>Keras_GLMM</th>\n",
       "      <th>GLMM</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td>3.186542</td>\n",
       "      <td>3.282506</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Stud_Age[T.4]</th>\n",
       "      <td>0.053660</td>\n",
       "      <td>0.073679</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Stud_Age[T.6]</th>\n",
       "      <td>0.062380</td>\n",
       "      <td>0.093308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Stud_Age[T.8]</th>\n",
       "      <td>0.141339</td>\n",
       "      <td>0.175263</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.2]</th>\n",
       "      <td>-0.067621</td>\n",
       "      <td>-0.086601</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.3]</th>\n",
       "      <td>-0.107238</td>\n",
       "      <td>-0.149557</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.4]</th>\n",
       "      <td>-0.167474</td>\n",
       "      <td>-0.198324</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.5]</th>\n",
       "      <td>-0.154622</td>\n",
       "      <td>-0.188271</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lect_Age[T.6]</th>\n",
       "      <td>-0.230838</td>\n",
       "      <td>-0.287835</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Service[T.1]</th>\n",
       "      <td>-0.094988</td>\n",
       "      <td>-0.091856</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               Keras_GLMM      GLMM\n",
       "Intercept        3.186542  3.282506\n",
       "Stud_Age[T.4]    0.053660  0.073679\n",
       "Stud_Age[T.6]    0.062380  0.093308\n",
       "Stud_Age[T.8]    0.141339  0.175263\n",
       "Lect_Age[T.2]   -0.067621 -0.086601\n",
       "Lect_Age[T.3]   -0.107238 -0.149557\n",
       "Lect_Age[T.4]   -0.167474 -0.198324\n",
       "Lect_Age[T.5]   -0.154622 -0.188271\n",
       "Lect_Age[T.6]   -0.230838 -0.287835\n",
       "Service[T.1]    -0.094988 -0.091856"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fe_df['GLMM'] = mdf.fe_params\n",
    "\n",
    "fe_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bài thực hành đến đây là hết, xin tạm biệt và hẹn gặp lại :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}