6.1-MultinomialLogitAndProbitModels/02-logit-estimators-results.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to use different estimators and how to interpret the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas  as pd\n",
    "import numpy as np\n",
    "import biogeme.database as db\n",
    "import biogeme.biogeme as bio\n",
    "import biogeme.models as models\n",
    "import biogeme.optimization as opt\n",
    "from biogeme.expressions import Beta, DefineVariable\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Import Swissmetro data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "pandas = pd.read_csv(\"../Data/swissmetro.dat\",sep='\\t')\n",
    "database = db.Database(\"data/swissmetro\", pandas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Use collumn names as variables**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "globals().update(database.variables)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exclude some unwanted entries**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "exclude = (( PURPOSE != 1 ) * ( PURPOSE != 3 ) + ( CHOICE == 0 )) > 0\n",
    "\n",
    "database.remove(exclude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Define some dummy variables**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "SM_COST = SM_CO * ( GA == 0 )\n",
    "TRAIN_COST = TRAIN_CO * ( GA == 0 )\n",
    "\n",
    "CAR_AV_SP = DefineVariable ('CAR_AV_SP', CAR_AV * ( SP !=0 ), database)\n",
    "TRAIN_AV_SP = DefineVariable ('TRAIN_AV_SP', TRAIN_AV * ( SP != 0 ), database)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Rescale some data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "TRAIN_TT_SCALED   = DefineVariable('TRAIN_TT_SCALED',   TRAIN_TT / 100.0, database)\n",
    "TRAIN_COST_SCALED = DefineVariable('TRAIN_COST_SCALED', TRAIN_COST / 100, database)\n",
    "SM_TT_SCALED      = DefineVariable('SM_TT_SCALED',      SM_TT / 100.0   , database)\n",
    "SM_COST_SCALED    = DefineVariable('SM_COST_SCALED',    SM_COST / 100   , database)\n",
    "CAR_TT_SCALED     = DefineVariable('CAR_TT_SCALED',     CAR_TT / 100    , database)\n",
    "CAR_CO_SCALED     = DefineVariable('CAR_CO_SCALED',     CAR_CO / 100    , database)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create parameters to be estimated**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ASC_CAR = Beta('ASC_CAR',0,None ,None ,0)\n",
    "ASC_TRAIN = Beta('ASC_TRAIN',0,None ,None ,0)\n",
    "ASC_SM = Beta('ASC_SM',0,None ,None ,1)\n",
    "B_TIME = Beta('B_TIME',0,None ,None ,0)\n",
    "B_COST = Beta('B_COST',0,None ,None ,0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Define the utility functions**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "V1 = ASC_TRAIN + \\\n",
    "     B_TIME * TRAIN_TT_SCALED + \\\n",
    "     B_COST * TRAIN_COST_SCALED\n",
    "V2 = ASC_SM + \\\n",
    "     B_TIME * SM_TT_SCALED + \\\n",
    "     B_COST * SM_COST_SCALED\n",
    "V3 = ASC_CAR + \\\n",
    "     B_TIME * CAR_TT_SCALED + \\\n",
    "     B_COST * CAR_CO_SCALED"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Associate utility functions with alternatives and associate availability of alternatives**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = {1: V1,\n",
    "     2: V2,\n",
    "     3: V3}\n",
    "\n",
    "av = {1: TRAIN_AV_SP,\n",
    "      2: SM_AV,\n",
    "      3: CAR_AV_SP}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Define the model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "logprob = models.loglogit(V, av, CHOICE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Define the Biogeme object**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "biogeme  = bio.BIOGEME(database, logprob)\n",
    "\n",
    "biogeme.modelName = \"swissmetro_logit_estimators\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Define the algorithms to estimat the maximum likelihood**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "algos = {'CFSQP                ': None,\n",
    "         'scipy                ': opt.scipy,\n",
    "         'Trust region (dogleg)': opt.newtonTrustRegionForBiogeme,\n",
    "         'Trust region (cg)    ': opt.newtonTrustRegionForBiogeme,\n",
    "         'TR-BFGS              ': opt.bfgsTrustRegionForBiogeme,\n",
    "         'Simple bounds Newton ': opt.simpleBoundsNewtonAlgorithmForBiogeme,\n",
    "         'Simple bounds BFGS   ': opt.simpleBoundsNewtonAlgorithmForBiogeme,\n",
    "         'Simple bounds hybrid ': opt.simpleBoundsNewtonAlgorithmForBiogeme}\n",
    "\n",
    "algoParameters = {'Trust region (dogleg)': {'dogleg':True},\n",
    "                  'Trust region (cg)': {'dogleg':False},\n",
    "                  'Simple bounds Newton ': {'proportionAnalyticalHessian': 1.0},\n",
    "                  'Simple bounds BFGS   ': {'proportionAnalyticalHessian': 0.0},\n",
    "                  'Simple bounds hybrid ': {'proportionAnalyticalHessian': 0.5}}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Estimate the model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "biogeme.generateHtml = False\n",
    "biogeme.generatePickle = False\n",
    "\n",
    "results = {}\n",
    "for name,algo in algos.items():\n",
    "    p = algoParameters.get(name)\n",
    "    results[name] = biogeme.estimate(algorithm=algo,algoParameters=p)\n",
    "    g = results[name].data.g\n",
    "    biogeme.createLogFile()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Print results for the different algorithms**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm\t\tloglike\t\tnormg\ttime\n",
      "+++++++++\t\t+++++++\t\t+++++\t++++\n",
      "CFSQP                \t-5331.25\t0.0004\t0:00:00.279782\n",
      "scipy                \t-5331.25\t6.5e-05\t0:00:00.085049\n",
      "Trust region (dogleg)\t-5331.25\t0.03\t0:00:00.070050\n",
      "Trust region (cg)    \t-5331.25\t0.03\t0:00:00.068124\n",
      "TR-BFGS              \t-5331.25\t0.013\t0:00:00.251144\n",
      "Simple bounds Newton \t-5331.25\t0.00063\t0:00:00.069949\n",
      "Simple bounds BFGS   \t-5331.25\t0.015\t0:00:00.269823\n",
      "Simple bounds hybrid \t-5331.25\t0.00038\t0:00:00.081564\n"
     ]
    }
   ],
   "source": [
    "print(\"Algorithm\\t\\tloglike\\t\\tnormg\\ttime\")\n",
    "print(\"+++++++++\\t\\t+++++++\\t\\t+++++\\t++++\")\n",
    "for name,algo in algos.items():\n",
    "    print(f'{name}\\t{results[name].data.logLike:.2f}\\t'\n",
    "          f'{results[name].data.gradientNorm:.2g}\\t'\n",
    "          f'{results[name].data.optimizationMessages[\"Optimization time\"]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Print results**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ASC_CAR   =\t-0.15463245\t-0.15463283\n",
      "ASC_TRAIN =\t-0.70118664\t-0.70118664\n",
      "B_COST    =\t-1.0837907\t-1.0837909\n",
      "B_TIME    =\t-1.2778603\t-1.2778601\n"
     ]
    }
   ],
   "source": [
    "results_1 = results['scipy                ']\n",
    "results_2 = results['CFSQP                ']\n",
    "\n",
    "betas_1 = results_1.getBetaValues()\n",
    "betas_2 = results_2.getBetaValues()\n",
    "\n",
    "for k in betas_1.keys():\n",
    "    b_1 = betas_1[k]\n",
    "    b_2 = betas_2[k]\n",
    "    print(f\"{k:10}=\\t{b_1:.8g}\\t{b_2:.8g}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Explore the results**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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>Value</th>\n",
       "      <th>Std err</th>\n",
       "      <th>t-test</th>\n",
       "      <th>p-value</th>\n",
       "      <th>Rob. Std err</th>\n",
       "      <th>Rob. t-test</th>\n",
       "      <th>Rob. p-value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ASC_CAR</th>\n",
       "      <td>-0.154632</td>\n",
       "      <td>0.043235</td>\n",
       "      <td>-3.576518</td>\n",
       "      <td>0.000348</td>\n",
       "      <td>0.058163</td>\n",
       "      <td>-2.658586</td>\n",
       "      <td>0.007847</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ASC_TRAIN</th>\n",
       "      <td>-0.701187</td>\n",
       "      <td>0.054874</td>\n",
       "      <td>-12.778137</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.082562</td>\n",
       "      <td>-8.492846</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_COST</th>\n",
       "      <td>-1.083791</td>\n",
       "      <td>0.051830</td>\n",
       "      <td>-20.910412</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.068225</td>\n",
       "      <td>-15.885522</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B_TIME</th>\n",
       "      <td>-1.277860</td>\n",
       "      <td>0.056883</td>\n",
       "      <td>-22.464577</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.104254</td>\n",
       "      <td>-12.257125</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              Value   Std err     t-test   p-value  Rob. Std err  Rob. t-test  \\\n",
       "ASC_CAR   -0.154632  0.043235  -3.576518  0.000348      0.058163    -2.658586   \n",
       "ASC_TRAIN -0.701187  0.054874 -12.778137  0.000000      0.082562    -8.492846   \n",
       "B_COST    -1.083791  0.051830 -20.910412  0.000000      0.068225   -15.885522   \n",
       "B_TIME    -1.277860  0.056883 -22.464577  0.000000      0.104254   -12.257125   \n",
       "\n",
       "           Rob. p-value  \n",
       "ASC_CAR        0.007847  \n",
       "ASC_TRAIN      0.000000  \n",
       "B_COST         0.000000  \n",
       "B_TIME         0.000000  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "results = results['scipy                ']\n",
    "\n",
    "Results = results.getEstimatedParameters()\n",
    "display(Results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Get the general statistics**\n",
    "* Number of estimated parameters($K$)\n",
    "* Sample size ($N$)\n",
    "* Number of excluded observations\n",
    "* Log likelihood of the sample with the default values for the parameters ($\\mathcal{L}^i)$)\n",
    "* Log likelihood for the final estimated model ($\\mathcal{L}^*)$)\n",
    "* Likelihood ratio:\n",
    "\\begin{align}\n",
    "-2 (\\mathcal{L}^i-\\mathcal{L}^*)\n",
    "\\end{align}\n",
    "* Rho-square for the init model\n",
    "\\begin{align}\n",
    "\\rho^2 = 1- \\frac{\\mathcal{L}^*}{\\mathcal{L}^i}\n",
    "\\end{align}\n",
    "* Rho-square adjusted for the init model\n",
    "\\begin{align}\n",
    "\\rho^2 = 1- \\frac{\\mathcal{L}^* - K}{\\mathcal{L}^i}\n",
    "\\end{align}\n",
    "* Akaike Information Criterion\n",
    "\\begin{align}\n",
    "2 K - 2 \\mathcal{L}^*\n",
    "\\end{align}\n",
    "* Bayesian Information Criterion\n",
    "\\begin{align}\n",
    "K N- 2 \\mathcal{L}^*\n",
    "\\end{align}\n",
    "* Final gradient norm\n",
    "* Iterations\n",
    "* Optimization Time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of estimated parameters               = 4\n",
      "Sample size                                  = 6768\n",
      "Excluded observations                        = 3960\n",
      "Init log likelihood                          = -6964.662979192295\n",
      "Final log likelihood                         = -5331.252006915832\n",
      "Likelihood ratio test for the init. model    = 3266.8219445529267\n",
      "Rho-square for the init. model               = 0.23452835796311466\n",
      "Rho-square-bar for the init. model           = 0.23395403010088356\n",
      "Akaike Information Criterion                 = 10670.504013831664\n",
      "Bayesian Information Criterion               = 10697.783857436807\n",
      "Final gradient norm                          = 6.477959424450974e-05\n",
      "Nbr of threads                               = 8\n"
     ]
    }
   ],
   "source": [
    "gs = results.getGeneralStatistics()\n",
    "\n",
    "for k,v in gs.items():\n",
    "    print(\"{}= {}\".format(k.ljust(45),v[0]))"
   ]
  }
 ],
 "metadata": {
  "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}