Box-Tidwell-Test-in-R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic Regression Assumptions\n",
    "### Box-Tidwell Test\n",
    "- Using R to perform Box-Tidwell Test (or Assumption 2: Linearity of independent variables and log odds)\n",
    "- Details on how to setup R kernel in Jupyter notebook: https://dzone.com/articles/using-r-on-jupyternbspnotebook\n",
    "- We are using R here because Python does not have any native packages that executes Box Tidwell transformation directly\n",
    "___\n",
    "\n",
    "- One of the important assumptions of logistic regression is the linearity of the logit over the covariates. When this assumption is not met, a Box-Tidwell transformation might be helpful. \n",
    "- The Box-Tidwell Test is used to check this assumption by testing whether the logit transform is a linear function of the predictor, effectively by adding the non-linear transform of the original predictor as an interaction term to test if this addition made no better prediction.\n",
    "- A significant p-value in the Box-Tidwell transformation means that the linearity assumption is **violated**\n",
    "- If one variable is indeed found to be non-linear, then we may need to incorporate higher order polynomial terms for that variable in the regression analysis to capture the non-linearity\n",
    "\n",
    "#### References\n",
    "- https://stats.stackexchange.com/questions/503906/solution-in-case-of-violation-of-the-linearity-assumption-in-the-logistic-regres\n",
    "- https://www.youtube.com/watch?v=sciPFNcYqi8&ab_channel=MikeCrowson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Installing package into 'C:/Users/kenne/Documents/R/win-library/4.1'\n",
      "(as 'lib' is unspecified)\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "package 'car' successfully unpacked and MD5 sums checked\n",
      "\n",
      "The downloaded binary packages are in\n",
      "\tC:\\Users\\kenne\\AppData\\Local\\Temp\\RtmpsdMeYP\\downloaded_packages\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Installing package into 'C:/Users/kenne/Documents/R/win-library/4.1'\n",
      "(as 'lib' is unspecified)\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "package 'tidyverse' successfully unpacked and MD5 sums checked\n",
      "\n",
      "The downloaded binary packages are in\n",
      "\tC:\\Users\\kenne\\AppData\\Local\\Temp\\RtmpsdMeYP\\downloaded_packages\n"
     ]
    }
   ],
   "source": [
    "# Installing car package to make use of the box tidwell function within\n",
    "install.packages(\"car\")\n",
    "install.packages(\"tidyverse\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: carData\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(car)\n",
    "library(ggplot2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 6 × 8</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>Survived</th><th scope=col>Age</th><th scope=col>TravelAlone</th><th scope=col>Pclass_2</th><th scope=col>Pclass_3</th><th scope=col>Embarked_Q</th><th scope=col>Embarked_S</th><th scope=col>Sex_male</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>0</td><td>22</td><td>0</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>1</td><td>38</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>1</td><td>26</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>1</td><td>35</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>0</td><td>35</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>0</td><td>28</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>1</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 6 × 8\n",
       "\\begin{tabular}{r|llllllll}\n",
       "  & Survived & Age & TravelAlone & Pclass\\_2 & Pclass\\_3 & Embarked\\_Q & Embarked\\_S & Sex\\_male\\\\\n",
       "  & <int> & <dbl> & <int> & <int> & <int> & <int> & <int> & <int>\\\\\n",
       "\\hline\n",
       "\t1 & 0 & 22 & 0 & 0 & 1 & 0 & 1 & 1\\\\\n",
       "\t2 & 1 & 38 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n",
       "\t3 & 1 & 26 & 1 & 0 & 1 & 0 & 1 & 0\\\\\n",
       "\t4 & 1 & 35 & 0 & 0 & 0 & 0 & 1 & 0\\\\\n",
       "\t5 & 0 & 35 & 1 & 0 & 1 & 0 & 1 & 1\\\\\n",
       "\t6 & 0 & 28 & 1 & 0 & 1 & 1 & 0 & 1\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 6 × 8\n",
       "\n",
       "| <!--/--> | Survived &lt;int&gt; | Age &lt;dbl&gt; | TravelAlone &lt;int&gt; | Pclass_2 &lt;int&gt; | Pclass_3 &lt;int&gt; | Embarked_Q &lt;int&gt; | Embarked_S &lt;int&gt; | Sex_male &lt;int&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|\n",
       "| 1 | 0 | 22 | 0 | 0 | 1 | 0 | 1 | 1 |\n",
       "| 2 | 1 | 38 | 0 | 0 | 0 | 0 | 0 | 0 |\n",
       "| 3 | 1 | 26 | 1 | 0 | 1 | 0 | 1 | 0 |\n",
       "| 4 | 1 | 35 | 0 | 0 | 0 | 0 | 1 | 0 |\n",
       "| 5 | 0 | 35 | 1 | 0 | 1 | 0 | 1 | 1 |\n",
       "| 6 | 0 | 28 | 1 | 0 | 1 | 1 | 0 | 1 |\n",
       "\n"
      ],
      "text/plain": [
       "  Survived Age TravelAlone Pclass_2 Pclass_3 Embarked_Q Embarked_S Sex_male\n",
       "1 0        22  0           0        1        0          1          1       \n",
       "2 1        38  0           0        0        0          0          0       \n",
       "3 1        26  1           0        1        0          1          0       \n",
       "4 1        35  0           0        0        0          1          0       \n",
       "5 0        35  1           0        1        0          1          1       \n",
       "6 0        28  1           0        1        1          0          1       "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Import processed dataset (generated from Logistic_Regression_Assumptions.ipynb)\n",
    "df_titanic <- read.csv(file = 'data/train_processed.csv')\n",
    "df_titanic <- df_titanic[ , -which(names(df_titanic) %in% c(\"Fare\"))]\n",
    "head(df_titanic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Box-Tidwell Test with `car` library's `boxTidwell`\n",
    "- We can run Box-Tidwell test directly using the boxTidwell function within the `car` package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " MLE of lambda Score Statistic (z) Pr(>|z|)  \n",
       "      -0.30553              2.4529  0.01417 *\n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "iterations =  7 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Running Box-Tidwell test with `boxTidwell`\n",
    "boxTidwell(formula = Survived ~ Age, \n",
    "           other.x = ~ TravelAlone + Pclass_2 + Pclass_3 + Embarked_Q + Embarked_S + Sex_male, \n",
    "           data = df_titanic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- p-value = 0.01417, which is <0.05 (meaning that it is statistically significant)\n",
    "- This means that there is non-linearity in the Age feature, and the assumption has been violated\n",
    "- We can resolve this by including a polynomial term (e.g. Age^2) to account for the non-linearity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### 2. Box-Tidwell Test (add interaction terms directly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "glm(formula = Survived ~ Age + (Age:log(Age)) + TravelAlone + \n",
       "    Pclass_2 + Pclass_3 + Embarked_Q + Embarked_S + Sex_male, \n",
       "    family = binomial(link = \"logit\"), data = df_titanic)\n",
       "\n",
       "Deviance Residuals: \n",
       "    Min       1Q   Median       3Q      Max  \n",
       "-2.7718  -0.6597  -0.3871   0.6388   2.4610  \n",
       "\n",
       "Coefficients:\n",
       "             Estimate Std. Error z value Pr(>|z|)    \n",
       "(Intercept)   4.67512    0.54968   8.505  < 2e-16 ***\n",
       "Age          -0.18594    0.06988  -2.661 0.007791 ** \n",
       "TravelAlone   0.21111    0.19919   1.060 0.289224    \n",
       "Pclass_2     -0.96542    0.27101  -3.562 0.000368 ***\n",
       "Pclass_3     -2.34673    0.25961  -9.039  < 2e-16 ***\n",
       "Embarked_Q    0.01839    0.37310   0.049 0.960678    \n",
       "Embarked_S   -0.53215    0.23650  -2.250 0.024444 *  \n",
       "Sex_male     -2.62753    0.20036 -13.114  < 2e-16 ***\n",
       "Age:log(Age)  0.03508    0.01599   2.194 0.028232 *  \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "(Dispersion parameter for binomial family taken to be 1)\n",
       "\n",
       "    Null deviance: 1171.07  on 875  degrees of freedom\n",
       "Residual deviance:  779.05  on 867  degrees of freedom\n",
       "AIC: 797.05\n",
       "\n",
       "Number of Fisher Scoring iterations: 5\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lreg <- glm(Survived ~ Age + (Age:log(Age)) + TravelAlone + Pclass_2 + Pclass_3 + Embarked_Q + Embarked_S + Sex_male, \n",
    "            data=df_titanic, \n",
    "            family=binomial(link=\"logit\"))\n",
    "logodds <- lreg$linear.predictors\n",
    "summary(lreg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAANlBMVEUAAAAzMzNNTU1oaGh8\nfHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enr6+vw8PD////agy6EAAAACXBIWXMA\nABJ0AAASdAHeZh94AAAgAElEQVR4nO2dgXbcKBJFNU42mdmZjcf//7Nrx+6WugWoqnhAFbx3\nzmayiXRFI65EIbmzvTEMU51tdAMYZoZQJIYBhCIxDCAUiWEAoUgMAwhFYhhAKBLDAEKRGAYQ\nk0iv5Vz9vTkEE+wLTJEIJhgApkgEEwwAUySCCQaAKRLBBAPAFIlgggFgikQwwQAwRSKYYACY\nIhFMMABMkQgmGACmSAQTDABTJIIJBoApEsEEA8AWkV4oEsEEP22iF+mFIhFM8PMmapFeeEci\nmODTJlqRXji1I5jg8yZmkf74iGwfhlknMpFe3nhHIpjg8yY6kV7uv1Akggk+bKIU6TMUiWCC\nHzfRibTfligSwQQfNqFIBBNcD6ZIBBMMAFtEOqS+AbYQTLAvMEUimGAAmCIRTDAATJEIJhgA\npkgEEwwAUySCCQaAKRLBBAPAFIlgggHgziJt22b8LBdgVAgm2ATuK9K2gUya8lwQHBjcVaRt\nQ5k05bkgODCYIhFMMABMkQgmGADuKhJrJIJnBfcViat2BE8K7iwSKgQT7AtMkQgmGACmSAQT\nDABTJIIJBoApEsEEA8AUiWCCAWCKRDDBADBFIphgAJgiEUwwAEyRCCYYAKZIBBMMAFMkggkG\ngCkSwQQDwBSJYIIBYIpEMMEAMEUimGAAmCIRTDAATJEIJhgApkgEEwwAUySCCQaAKRLBBAPA\nFIlgggFgikQwwQAwRSKYYACYIhFMMABMkQgmGACmSAQTDABTJIIJBoApEsEEA8AUiWCCAWCK\nRDDBADBFIphgAJgiEUwwAEyRCCYYAKZIBBMMAFMkggkGgCkSwQQDwBSJYIIBYIpEMMEAMEUi\nmGAAmCIRTDAATJEIJhgApkgEEwwAUySCCQaAKRLBBAPAFIlgggFgikQwwQAwRSKYYAB4lEjb\ntuk+ixRcGYIJNoEHibRtdSZNeS4IDgweI9K2VZpk6jLJEac8yQS3By8kkuiQU55kgtuD1xFJ\ndswpTzLB7cFjRBpRI1EkghuCB4k0YNWOIhHcEDxKpMrU1UgFn1y1mOA44IVEej14lDXJV4sJ\nDgNeSaSvFOd4LltMsH8wRcKBiyF4bjBFwoGLIXhu8IIiFVcdfLaYYPfgFUUqrTo4bTHB3sFL\nivSV1BzPd4sJdgumSBSJYACYIlEkggHglUXaa6RdJwHY9naT864guBK8tEivu0c3Oa7Bxvdt\nvXcFwXXgtUX6zHGKdwm2/gRIjK4g2AqmSBSJYACYIlEkggFgivTKGongejBF+ghX7QiuBFMk\nggkGgCkSwQQDwBQpDa78Tok8GB6CXYApUhJc+y1HWTA+BLsAU6QUuPp793LgBiHYBZgipcAU\niWAluFKkSfMl0tdvR7eGCRXekQ7g43vhmDtT2K4gmFO7CvDRI71JiV3idgXBFKkabBMptU/4\nriC4vMlCIkmMyIsk9ylpn6+uIBgNXkgk0b3lGWwplijSguB1RJLN0k5gQ7FEkRYEUyQZWFUs\nsUZaD0yRZGDdqgNX7ZYDryNS5b8hW/lACf4S7D2u+nhd8EIiVf6r5rUetTLJVx8vC15JpHHg\nBu/u3ROsK2YFdxYJNZqCnQuKND24r0irvrlGkaYHdxUJN54En8x2oEbngjXS7OBpRTIeqdW5\nODdm0VnurOBZRbIeqtu5WHWWOyu4q0gdR493kbrOcgluD+4rUr/5DEUiuCu4s0ioaGokwXjd\nN6FIBJvA04pU+heXE5vKv/vbGNZIc4PnFekrgkv/cZN+Leaq3VRgijRKJIKnAlOkFiIlDqgB\nq25WDvqY4AVEGlEjpQ6JbbENrAvBGvD8Iu3X9/zwxK7aJW+C0HuoDawMwRrwAiLdIhqevkSS\n+OSqj9cFryOS7ELvSqRu6hNcDaZIaXDV4jSoRuqnPsHVYIqUBFc+LsWs2lGkQOB1RJK9NPS2\nbzr8n3VJtKLSUFUI1oAXEkn00pAnkc4NrZwz6kKwBrySSF8pWuJKpOcbUOUqhjIEa8AUKQkW\nTQNVQaxiUCS/YIqUBkumgaogVjEokl/wcJFs46qqy0q3m7fEpsAfG6rksUZyCx4tknFc1XVZ\n4XbjWySu2rkFDxbJOrAgXZY6+ACRKvmQR8glMDxTgilSUaQONVLtETCPkPNgfKYEU6SySM1X\n7arveadbHbrF8EwJHizS2O9xTIy97ucCLhL6HorPlODRIo39ZuHdo+5ffnJoAlQkdFXXIFOC\nh4tkCxR8HHv9WwyukSjSGDBFGi0SeNUOty4Y4+R5AfsRSXW2ZxIJDIatC8bvip5gNyLpzja2\ny0bWSHgwal1wgq7oCPYikvJsg7ts4KpdezBF6gGmSNODKVIPMEUCgCUNH9hi1kgdwF5Eivyt\niKKmj2wxV+3ag92IFPd7emU3U08tJhgP9iOSKp7AFIlgigQA9xQpcZxWXQF6Yy8RTycPBl5S\npMIY8V0jpY7UqI9RLxol4nRU1IEditR8Daw0Ri4fb2aAl9sgTvLh3tf6yRfulb1zKFIXkZpf\n34tj5PqFG2OwIjV/F4MiKcHuRGpfcWhEwo0nqEjJWxM0FEkJnlekLCSuSImflWg23Fkj6cDT\nilSgKGok3KUfumqXujWh08wjitRFJFCNVBxhhQNkayTbgN33wZ7k860JnynHezOwP5Ewa2DW\nEZZbtbPxDvuAT3JVs2SZcrw3AzsUSZLurzybeMedVn3ck2geReoskmIGltwZ+cqzqVjCVVjF\nY7QBh3uETJEyUawJZHY3HDQLthRLu0jubxyNwMkbuesWW8F+RSpOpgZ0maUqcbImAL6mKCIT\nCdU1FCkVbyJ9ReeEizUB8CxXE5FIsL6hSKn4F8lULMHTfd1FFUGNhOscipRMZY1ki+L6rrwz\nDauRhookWLWjSI1FQv+0gyjyisM0x2sQ5yJdgylSa5F8g5Xn38M9FA02hjXSKJEup9m4zCiS\n+y8/4apdH5EkK6ioKMC6C+nluro1fAHBBTiASKlLv4tzoTLg+kmvMeDFtb0xLvo4DJgiPaY7\nuLrUxj7uOTRmnj7uAaZIj1lcpENrBq4zRgQHEMlrjYQBexUJtpR2zjwn77hJAJF8rtqhwIkR\niym+NHnQJ/ejtzMsrjUDhxApHFi1BpbyCLIcqMlhrSH3mu0Uj3uagSlSA3BqyEV5QLVXSE+N\nqJ6D3uP75BnBPURqMNl2fS6SQy6KSA+tOP0BZLLn+uRZwR1EalG2Dugy+YeYQqTT58VN9iiS\nSSTclOAJ3CSI56aVIg2pkQRg2D2KIi0skupT1NVII1btRGDUPYoiUSQhJbHlhKPHuiAxYVf0\nq5HAJjkXSQGuzkCwaLKX2quqXflMLlILk1zXSDpwbUaCBZO91D4zvjJBkaTgyk9wA8ef5ZYi\n9GjcFyo1A1OkvuAJZrnFSCZ2FMko0vSjRw7GD6NoFQdFqhBp8vmMAgwfRvEqjngtdiQSPFHB\naJEiXt+j3UNlYIrUF9zihhRMpDnBFKkz+Dzsa0SgSF7AFGk0uM6ESSuOeGCKNBhce0+Zs+KI\nB6ZIg8F892gOMEUaDKZIc4DVIr28ZyWRTIM8sRNf4psbrBXp5f7LGiKZhnlqJ81LfKpDxu/j\nKcAUqQg2TbySOylarDtk+D6eA6wVabeJIuVSKdJhd8mhw/fxHGC7SH98RL5PzHwN6g47pXav\nwjCDIj9jKy02dKiRknvnvu30nPh9PAWYIl2Am6/aJffOfdtpIpiuqGuxLlOCDSIdPFpAJFBU\nPu4VUi+R9gPtBwzXx8FEOnrkUKTKV2Yatdj2sEi0F6LFu7KHQ0453puB1SI9eORPJNHYK2zS\npsXW1xd6rdoly7Ipx3szsFakl5eHVxvqG2BL1Qs3pU18iVQkfv2OIrkAq+9Ij6lvgC01IhW3\nCSISegaWWt+Ycrw3A1MkGbgyLTyCjvfDWgNrJAt4MpEkI3aESOAfG2p44+CqnQ08m0iSEdu/\nRkKDOQNzB55OJEm6r9rBwZyBeQMvKdIEYM7AnIEpEsEEA8AUieCm4MQ02nmLbWCK9JE5J0qV\nK4Xgl/jA4FQokjpOS3dXqxi1z66wL/GBwclQJHV8Lia7WlevfpuCImnAFAkn0pAnvSqRVGJR\nJA14uEi2i6ZzkTxU2IkxrLtFsUbSgEeLZJx+uKyR7kPXx+hJeyTvba7aacCDRbJO5H2u2h09\ncjCfSUzs+ou0CpgiIcH3iZ0LkZ6TnuwBwMpMCaZIeLBXkTKTPQBYlynBg0XyUSMlwLY1kH1n\nBzVSIqqiyUWLw4BHi+Rh1S4BNvp9S5AKmyLhwMNFsqUx2DrjvAQ3CEVyAaZIKfAiIslqJHBH\nOO2KWjBFSoFXEUmyagfvCaddUQn2I5LqdDmvkfJgfPpdU1BdErUrypt4EUk3dH2v2pXA8HQT\naT9DlZ0TtSvKmzgRacxTd6er1J7A9xOTMqoG3CAUadx7YB6fm/oCP+iTfQPKAMaHIo0RKXXM\nKU9yHXivkChSYRMnIg17xZ8iScGnWxMKDAtF+kj/VTuKpAMf1hrk5yqx6QRdkdiki0hO18AS\nC7tTnmQ0WOfRGgs6XURy+1TmOP0Xg22fZcrRc5nkNNB1i63gHiJVT6pzYFCO7bsGGz/LlKPn\nMkqRKkcJRVJnnEjWD+N6vDcD60SqHSYUSR2KFASsqZFCLwf2EMlvjfQVTY1EkXRRrNpRpMsG\nOF21u0ezascaqRWYItU0wJaRYK7atQKzRqpogC0ETwnmqp29AbYQTLAvMEUimGAAmCIRTDAA\nTJGQ4MIc32mLCQaBKRIQXFp18tliglHgziKhHij16zLtm85u/n0kgruC+4oEe8WhW5dpWryL\npHicXx2CXYC7ioR76a5Xl6lafN9Y84JZdQh2AaZIRbCuxUePlvghHIL3TShSCaxs8X1iR5FW\nA3cVafIa6WEf/gz7WuC+Is29anfcx/Iz7LYQ7ALcWSRUnIOPHol/ht0Ygl2AKZIabFkPx/9E\n1j1T9nE8MEXSgi3r4Q1+RvieGfs4IJgiKcH6J0vanZSZsI8jgv2IpBpoQURKFUvoTDks44FH\niXQaVpKR5mEx2VT3UKTpwX1EutZGMtR8LCan1rblO+Ez5bCMB+4ikkAbgUheFpNN07XzlvEe\nqRFc2qSDSBJtAol0bo7EifgveRBc3MSJSIKB5VYkkRPpl/gALjnoCoJHiZQae9eDykeNdM+D\nR5eNz4lUbZKHriB42FcWm0bQvpNiOVAXw4J9pUi1Jk05LOOBu4jU8SuLa4em4VzIiqVkjUSR\npgH3EQmeHLh6bFpafPRI/uUnZ5FsDZ9yWMYDUyQZuBjJHC8FPnnEb+cPCx4uEvYyPEYkycGT\n4JNHhlfFpxyW8cCjRUJfhgfUSI/HNn+L0MPeis8w5bCMBx4skvUO4mLVLnXsrAnyrlD2yZTD\nMh54OpEqUwc+eKT/8pMHDe2vTMBCsAZMkfDgpAmWl/hEfeO6K9YBDxZpyqWqpAm9XpnAhWAN\neLRIHR6e9F8DS5nQ65UJXAjWgIeLZIvh+g4Hl495P7T9y08oUiDw9CKNXAOzrmg/7P/8J88b\nOehjgimSHSw+uL4Rx/2TvGMc9DHBnkRqVMoMFQk0xzvQzp9lymEZD+xGpFaljHIMtzkXtXO8\nE+aQKYdlPLAXkdrdOI4eDXsqUzvHO1CuqyZQphzvzcDzi+RjMVnViBJFUDWBMuV4bwZeQCTV\nERqfi0qR0qsPjUyacrw3A3sRqf3jHhcioe8gFMkL2I1I7V9AGFkj7eBzG9DreCixphzvzcB+\nRFLFBB63alcCV8/1WlVNvk6ed/BKIt1TGGndW+y3anJ68pyCHYrU/MZRGmkjRbKN/ycwRRoD\n9idS81KmONQGimQUgCK5ALsTqf3iWvEeMK5GshrwDGaNNAS8tEiJIw04F7vUEJFKM8TmK6Pr\ngucVKQ8p3QPGnQuYSBeHaAHWZUqwO5FQNVIJU7gHDDwXiRZj112Uqk453puB/YmEGT2SQeNM\npPMHB6+7UKSGYIciSSIVqTxsEosOnk4yulykSA3B04t0ZdJ9YyHYGgcisUZqCK4UqT7vp7YV\n9yvSTZs0oyb4dnn8lNNltu+12w4RbSkF2+Ly2TTBMPBgkZTTdjn4AFctOvg6ydhVO/yxbfHV\nxyDwzCLpTNKADYkGNk4VJInWFRRJuOTgddVuJNh6YiQJ1hUy8GCRGtZIxwNwzVcbiqQEjxYJ\n87MDFweYec230WinSErwcJFs0YBV48FFixVpNtxZI+nAfUSqPCfn3Zt3GXwUNZ+BoVus4k19\nsfIjUuXVLbF76y7DX49bizS0xZNPn72IVDnfTu3euMsaVAiNRWpwa6pa0Ck1giJRpFowPA83\nJOytqUakYiMo0ooi9b++67I383RrqkyFSOVGUCSjSJFrpAHXdxv4fGvaY/sIFTUSRWojUtxV\nuxHXdyM432LjR6h4xECRGokETy9wIJFuSXtkmKfWtJg1kmuR9pNDkfLgxMTuOO2zgzUpHYgi\nDRbpMA5A4Os5Y5gaKR/rOsSU470ZOJBIx3GAAUtWMUAejTzJpXWIKrAxU4JXFik1sKY8ybaV\nkym7ohmYIi0g0i0FjxJ/4aHFccCBRILXSMuJlJ+nphRz0eIw4EgiwVftJDUSKq7ByUmf6xa7\nA4cS6R7UCgAvw59JiYTq43Ncd4UVHFIk2Jr0OVOe5MskRGIf68ARRcI9JT1nypN8nbRH7GM5\nmCI9ZsqTLEhiYsc+1oDnFck2DKY8yYZQJCU4okii+btxHEx5ki0R9R8vVvsmEUUSnEDrFXXK\nk2yKzCNerG6bhBTp+YdwzqFI7cHs4+MmoUXKn0ie5PZg9vFxk8gilc4kpx3NwRTpuMmsIrEQ\nbg/mxeqwybQi1YAbZEowL1b7JpFFcvuFqIlWTTl6CN43CS0S/s3KNi/coMCpEOwCHFwkl+Dk\njNN1iwmuBlMkPJgiLQimSHgwRVoQnBfpr5e3t3+2lz8nFKlQWVEkgk3grEh/bdvbr5f38VA0\nqb4BtlSBS2t9FIlgEzgr0rftn/f//fW/7WU2kYpPn05gy7ogV+3WA2dFer8h/b19+/3flUVK\nbGr78YI32Z76uOzj9cBZkV62Xz+2/31USSuLlNjW/l4M/vHxF7hNCNaAsyL9+X7SXz5uSD9n\nE+kwoC9nYGeR7G9qNnih6RPcKARrwFmR3n5uL3+/35iKHsUU6fXg0cWaAEUiWAbOiyRKfQNs\nwb3Jc/FNq+mZHUUi+HkTilQGp9YaWCMRfNokKdJ2zNoiJXczHIyrdnODVxbJ+N3fdpGahGAX\n4LRIH/nP919vb7++/6fkUWyRbHck+9SuTQh2Ac6K9J/t388/LppU3wBbxolkX2xoFIJdgLMi\nfU3p/p14akeRCIaBsyJ93z6ndhPfkUw1EkUiOLlJTqTfb36/5+XXvCKZ/n0k1kgEpzbJifT2\n789v2/btz39LHolEmu3hieDzJDaZcvQQvG+SFUkUQQMWfJyf+si+W0xwLbi5SG1ejLkdGE22\n8SRvEU05egjeN0mKBHwguzUx6evAaLKN97wXRVoQ3E2kFt/jqCRfbmpr6WkvirQgOC3SR1Bv\nNrgR6XrbapE+/3P6g73FDUKwC3BWJNybDU5EEmxcK9LXf09/sLe4QQh2Ac6KBHyzoYFJhhpJ\n0op9C017H7S5m/R0yClHD8H7JjmRgG82tBNJM95FrTh4pDPp8QinP5h09BC8b5ITCfhmQ0OR\nNBHdbmpaTJEWBmdFwr3Z4OZfX5HcbqrUP5VGrJGWAedFEkXYALRHVV1WtOQk0n1LyWc4lUb7\nTlOOHoL3TbqIBE9jkRL3FMUtKnkAl11BMAycF+lzavcTMbXDp7VIpypHM9mTiYS6S7vs4/XA\nWZEm/jGK0qrDW2JTrUiil1ZhdaPPPl4OnBXpx235+8d0Ilm+IDJZ91weIP8nxuWMRJz28Wrg\nrEi3B7GL/6h5VbGUwJSObYvvPl4GTJFkXxBpm+Plj0SRJgNnRZp5avcVkUjP24JEYo00GTgr\nUpvFBl9LVYmyRyOS/MOk7PPVFQRXgrMiNVn+9nYZPpc9WXBdsQT74OdMOSzjgfMiiaJrgNPC\n4NisPLiuWGrm0ZzDMh6YIinfLbUWS85nuQRXgrMibfe8/MiXSboGUKQnjq6dMjAsBGvAApHe\nkzVJ2QCno2dvlqB5phopyjWFYCs4K9Lft+Xv/779zC+BaxvgdD5z8Ej82kLdHYkiTQbOivTt\n9p0N30oPZesbYEsbsGp0UySCj5vkRDq+2UCRyhvf9ynsHGSWS7ARnBXp/p0N39/++bgrLSaS\nqVgqyRFllkuwDZwV6fBmw7b9tYZIyoWE52KpeEOL1hUE68BZkd7+/fP2ZsP2M+fRbCKZFhIS\nIiX21bRYdbMK18dzgvMiiWJuQOXEpnGXVYqU2lnRYsmx9y2mHJbxwINEqi21PYl0mhAm95a3\nWFKpHQ4w5bCMB86L1PI7G6oXf1t3ma55TxNCmEjZVhyPMOWwjAfOitT0Oxvci2SaeqJFyvcS\nRXIHzorU9Af7/ItkCbpGokiBwFmR2v6oufMayRhMjSRZO2SN5A08SCTnq3bWFO4kFnDhcrP/\nhdOuWA2cFWmB72xoBYaJJLrcuO6KdcBZkSb+gsjm4MoaSReCXYCzIs37lcUdwIk7ifMWE1wJ\nzouUyct71hbJVN5V1oSFTNnH8cBakV7uv6wqkmnBsXaVskhuwp3z5LUDp0XajqFIx5gegVU/\nNyuTG3Bfpzx5DcEUSQk+OCEfwqadlGB4Jjx5DcH2qd0fHxHtM1X2q8v5KoPdSQeWbo47NpMO\nayQZ+HYLUN0MTDtJuY2+73XGk9cOTJE+onpPYLt8gQe1k4xrmWBKEuXk+QBTpFfjm2vKuuft\neSdQFDSK1BBMkczvUh89uh6eb087wWJSHwtWZkowRbL/UIJqunYDw9fYLOqjwbpMCdaKNOOb\nDVaRErvn46IrVBa7aHEYsFqkx9Q3wBYHNdLj3pfFUpCuINgIpkgfUa3apfY2fEEkKiawfHmk\nQaYEUyQAWFIsuRrvmuURfKYEUyQceBcpMUw9jfdAVV0YMEVKgy2La/fxmRqmnoovioQHU6Qk\n2Pa45+gR4EfNH+5wyOKLIuHBFCkFVj66vOcw7JEiwYsv1khwMEVKgWEiYZYD46xiiDIlmCKl\nwEaRTjOxA8XXcqAkBGvAFCkJtnt0XLU7ClDVYnCNJAvBGjBFSoOrVu0SIll4CXAicft4KjBF\nwoEfFwden8yqMIki+QdTJCD4tMgmWi0QcimSazBFQoJzq3a1czzWSO7BFAkPLolkvDMVdvLV\nFZJP56vFIDBFagA+L3+D5nipuOoK0adz1WIUmCK1AGdndJOLlPh4iU/rqcUwMEVqCD6Oq7pX\nJgoRtNh2RIhIqY8b4+QpwRSpIfgs0pAvPzEeEyFS8sIR4+QpwRSpJfgwivJPeivNumyx9S6I\nqJEoEkWCZB9EWZFq71GuRHr+eBSJIrUBp4rxOpN8iZQ8eAtwKhRJnajgTDXeVKSONVL64G3A\nidS8KCkClzahSD3BQ0Tqt2rnAzxiQYci9QWnH7QkyiYtuEGigsc8Yggh0kTTbIk2unGgafES\nhlKkXKZa+Lle/lYOBEWLC9wRpUwjMEXKJNUxUU9yIqV1PMl4kLe4MML2v9g3CNvHrJHSCSRS\n4QTmwKXlB9GIQIiUOiR86esertpZRKr/UHFEKg18gUiHIax4VRwqklJiW3yevEpwa5EQ5yNK\njVQc+NcinW4GaJHyJ+MsUsJvWFyevFpwY5FkY+EqZ4LLc2ESKXX/aSVS3olTI1remlyevFpw\nCJESYDgRALaJVPgp2jMugce+gFAwChaXJ68WTJGA4NJnlb8SV5hVpQ7Q6DXbxI0SFZ8nrxLc\nWKRWNavTc1H4rPJX4vIiJf+mVVeo5peqOD15deDWIjVaRR3QZaYPohrv+TnevsH5b6Z8KhMP\n3FykNnHw3FQS23jPzqoksz9cRj6ViQemSDKwaYZz3MmwuHY+5lmxRlPnj2BXMeDgRCiSOnOL\nlD/ms2Itb02411JagFOhSOosIVLhdtPj1oTo41S3tZ+MjgBTJCG4Y430sH/5b1K3JlQaieR+\nMmoDUyQp2HT6950ale7nWxMsbURq0NB7KJI6YcHwUZRdmKh2tk2NpBRJ9RkokjpRwe1uHEmP\nqo7UZtVO1wO6z0CR1IkKbjgDqxuxBTA4ao/kn4EiqRMV3LGU8SpSabZWeTGgSOqEBeMrbY1I\nziuO2oKKIqkTFwxfsTq1+HaE2qrJw7M65y0+bkKRJgPvYy/aRCngPfSwCUWaC1ywJX3Fl4Jx\nCVfVORQJNbEJP96bgVUiFQeugxoJBa6NO5FgpXb48d4MXLqsJz3KvvIa5e0RCbgy3kTCLf6G\nH+/twKUuzi6KpXaK3xU9wRRpOrD8Oc39fCRPzARd0RFMkVqAfUyUTtmeT8DDDYki1YC7irRK\njQSbKEl6Sw7etpRJh7+xgpWZEtxXpLGrdthhWTwOZliKrjsWkbbDH2WPNOV4bwbuLBIqDobl\nxYEAIslmwjUi7b93OhmNA15HJPSwVB6prsWFhhtqpBusfABPJ88/mCJVgzOHqgfvLS61XANO\n3ZC4/A0BU6RqcPpYCPDDML/9rhb8WCF9cmFVnShTgtcRqV+NhAPv867DgEeAdz6XvzHghUTq\ntmoHB9/GOX68P5uaAku6TRWXfVwLXkmkuODmM7DCrW7/K5RRPvu4EkyR7GDVwKprcfMZWLb4\nkq14qOLi5KHBFMkM1g0sT8uBsjyJJFuskYMbhCKpAwLXDEvlwHK1HCiKQCSbWN5HhQm8tEhV\nE6UhIg0AF1b2jOVT2K4obrKwSCkVKNIJfJDl+seZwq6M1oIpklUkt99w02xxLTGxe571ifrE\n96gwgimSWSSn85l+i2tnkWR3ad+jwgheWaR9zO0n33mLr8E9F9dO5ZPo4Kg75jkUSR3sqt3h\n7Htv8UfTIY4AABQlSURBVCW4p0in8kly8EQ91WwyikoEkWydCO2y4+mnSDlwMdkFidSWz/VU\nv8loQ/BokYydSJGy4FKPtq/qJBM7QWFlu1ktLJL1+kmR8uCiR8PXGUUiGW9WrfpY0gaK9Ppw\nuiYQKZtDZ4uGhr1Nl614LYpUvFkV0rjF5WNTpFe/q3bY8Z4oTzBgXW6Hzt92fIkka8RgkVzU\nSA3BdYaCx7tyfHZ7ZULypFcFxiaGSB5W7dqB6+aM8PF+mlU9/FUFWBfF5VV3nV1aJFv6XS1r\ncjwFHkR6eG52rkoqwKpoLq+q87BwjWRNr1ee6+JOpAfy1ZHCnjx4JGOBIqXAyrl5LgeMiQau\nkR7JpwPlHuXAE3VUlDehSAkwSKTqh/eiS6EemzpOomhqZVLUUVHehCIlwCiR8kUJLrgX4U9/\n4LnF3sAUKQkGj6LaOV4psBfhn/9/s8le2FFR3IQiJcHY8b4PS/x1vtuaL6zpcUdFaROK1AP8\n4BHWpF5rvrh7VLSTR5EcgdsVS73WfHH3qHAnjyIhwZXjP7mKgXCqV1ecRbJeFWYaFfsmFEkG\nrr2TpFYxIHenbl0BW9mbaFQcNqFIInD1nCyxioGZ5/XrCsnKngmMyvwiOV3zvWdvX3ORkMwc\nGJPlX91XgStFkuV3h7fDV0fQvq9B//Vb4GFdd8xVYre+UdrdkbwvVR3bd10jAS/DsWokgoub\nUCThS9qVC9hJcKRVu65gW8dQJHUGiHTeWPWZphzvrcABi68eIrl/L+bQPoVIuk8143hvBTZe\neeFLWvd4EWmCVbuHbS1v+0w43puBbSLhL9f3uBEJHgfzd4rUDGwSqUEBcU9wkQrd4mD0KIsl\nBy2OAzbfkChSIqV+aX6SBadEtyA+43hvB7ZVSBQplWLHtD7JonOimuNNOd5dgVkjpTNSJNXV\njSI5AS+wamdJRJF8V3UENwR3Fklz0RhYI+nm26JiacrRQ/C+SVeRdNPYgdd3U0OHVnUELyQS\nbmHFw6rdKfvHS+w95egheN8kgkjnfVyei/vHS31Kly0mGAaOIFJiJ5/n4uhRjBYTjAJ3Fcm2\n1J8alwO6TNLy+8RueZFKveWzxZXgviJVlh5nMPzRQfYT1a6HTzl6sin2lssW14I7i2RJQSTU\n2sUJLGpENqn18ClHTy7l3vLY4mpwAJHyNRJuFfAJnGmD5snS805Tjp5cKJJLkbKrdmNEUvmk\n30kbj8OSIvkUKQfuKJLpy0/29uHnoPe4HJaskUKJ1PVn2BPTtcs8eDTuFf8B4NKn9dniSvAo\nkSqHVf9Vu1tqi6X+LSa4B3iQSLUX6HFdVlks9byHEtwRPEak6qnOwC6zFEst53hTDst4YIqk\nBluKpYZzvCmHZTwwRbKClZ/heZ0Rdmdy0BUEs0ayg00iNZjjOegKgsOv2sGjAOtMeFpnxM3x\nPHQFwcGfI40Fqwx4AuPmeC66gmCKNAgMm+PF74opwBRpFBg1x5ugK2YAU6TB4Oo53jxdERo8\nXKQpL8OJD3W5gGm9NTnvilXAo0Wa8jKc+lCKt2Fty4HwEKwBDxbJWmq7PhfJDyXvCtsDKnwI\n1oApEh5sE8k6x3PdFeuAKRIenDTB8hKfqG9cd8U64MEiTV4jHT6d4ZUJ2VXGd1csAx4t0tSr\ndkcTDK9MUKRA4OEi2RIDbBUpsXs+MbpiejBFagiuFIk1UiAwRWoJttVID/tfbhOkK2YHU6Sm\n4POqna0mFIDhIVgDnkUk0Phs3mLjKuVzam91khCsAU8iEmh8Nm+x9bnZU2qLL1EI1oDnEAk0\nPs9gXE4i1bSXIrkDU6QiGJdnkaoaTJHcgSlSEYzLU41U2WLWSN7Ac4gUpkYyfiveOYVVuyjr\nLnOBJxEp3Ohpdw+Nc02ZCjyLSOHArcZ7nFnuXGCHIi3yOP/8MW3jnyK5APsTadkXzIwCUCQX\nYHciyQbChOfCagBrJBdgiuQFjBIp3LrLHGCK5AUMEwkVgjVgdyKxRnr8IwA4faw2YEmmBPsT\nad2TnPLo+rJSdbEqwafs42ZghyJJAgYX3hNAxT7LvTCpZvpchHvqCv9givTq9s217SF48MW7\n6J66wj/Ys0i9rpbHkerpJG8ikyAipQ7gqSv8gx2LhBo9lzWXV5FkJlXUSI9TPP13wxozJdiv\nSKjRc12uuxVJNLmrWdB5uCE9HeDy+mOOrz4GgacX6aLC2LfRgpUxg5uIdIAfjnE67HEbXPz1\nMQBMkb42UoOVcTHLlR8gVT6hjPLZx5VgvyKBRo9IpOPWYrIuztddEhO75/JJ04/FUKS+IoFG\nj9KjViYFGz1nkXRXpFKCdYUM7FkkEFjpUf13/CTjoisUOa1D4BYkonUFRdKmwTTmnmBdsUty\nMuphE8BPUOFCkdTpNZ/BJVhXHMDPRt1j7aW4XVHahCIdUrj61qZfVzRbXMsuSNSCUaFI6rRa\nxdgrpLAi7S2v/ASXLdb10r6l01FRB15RpNLpf7vexJZeXYGr865bbFwQ9TkqKsELilS8kD4X\nBrB0F6l6CVLQYqVHudewmk1GUaFIqYhEOm5dcawCGJVrkTwsrpVEgt3+KZIwmGm2SqR4J/m0\nYOJiTaAgkqh9og9AkWRBTbMFNdLDpsEe5x/WGqqWIHudPEn7ZJ+AIolSnGZrSbm/mkCke/YK\nyYFI+emEoH2yj4CahZ9DkdSZSaRbPNRIJbDwhiT6+cy6luVCkdSJXyMlYvsIjh4hi2d/jUya\nS6QxjyIezk3FiXK5gBkGLL0hUSRRstNsXIrgmjPlaVgGBIsqJIqkzhBw1amaqys8glkj6TNm\nfYYi+QZz1U6bQdeeZ5FUjQg3LAlWgSOKNGw2fPJI0YgpRw/B+yYU6Qkt/WtlK6YcPQTvm1Ck\nZ3KbVkw5egjeNwkoUrMaSeUGRfL/pLcjOKRI+wnECmUVSbLPjKPHeEGbsSuCinQDg29NRpFE\nO004eqxT7Am7IrhI8GJJhXvw6HK3CUcPRTpuQpGOUdH2iR1FgoKtCSXSy3vmFcn8QuWKIhln\n1o2eXLzGEunl/osDkZx82c+qNZLNiVZLrq8UyZAbGH5OTC0+tSLRrClFMqTBLOKeSCLtNrkQ\nySc4NVR8t7gfmCI9i/THR+T7LJSvsTK6GT4zf+fIP5ubxQZlej033e4Z/qOIDsGskaKL1G1N\nYHsMDpxMOPDqq3b3de+DR5FE6rdKvaVM8tQVBOPB+jvS0SOKVDgSRVoJrBbpwSOKVDoURVoI\nrBXp5eXh1Yb6Btjiukb6LAMeRZqzMCB430R7R3pMfQNsqVu1K4xqYIuPIu0Sw43y1cfLglcS\n6ZbSrQnb4geP5D9yoYrPPl4OvKBIxWIJ3OJ9Yvd5SFmhporLPl4PTJFwYMEhG1RNLvt4PfC8\nImUH6gCRWs7xphyW8cDTilQYqKXKv/Xj/LNRtZlyWMYDzypScaA+j2oN2Job+LlqgoHhIVgD\nXlKkwja93oalSJOBKVK9SCIlnsGwJ0tTDst44FlFsv7TVP3ePSrML1WZcljGA08rkuRCn1iK\nrhOpcMwcuHqON+WwjAeeVyRJjpW/GSxb0qZIc4PXFukzx7FcUyM9GPW8kVAkvVMx+nh6MEWq\nFyn9ItDTNlnws0dak2L08fRgioR7b+fGSc7W8i1+9kgpU4w+nh5MkV4fFx2qTcqUPaIWbwaT\ngvTx7OCVRMoPz8e5mV2m51rpEK1I4ib46uNlwQuJdD08LTeEPObxD2UtpkhRweuIJBiflhuC\n9EjCFlOkoGCKlNimVqXU7ooWs0YKCKZIqY0AJj3/kabFqmN76uOFweuIJHtN9A2w5JADtwnB\nLsALiSR6TfTt9WCSsXFZcJMQ7AK8kkhfKVry9rwJyieXXUEwDEyRkuCjR/xqBYIFm1CkJPjo\nkesfZN2+3qXAk12ePLfgBUXSfPmJe5HarIz8js+T5xW8okiKLz/BFUttugKzWp+O05PnFLyk\nSF9JDcBnMKxYokhzgynSxXc2oIolijQ3mCLJvvzEqUiskbyAVxZJ8wWR1cUSV+3mBi8tkuYr\ni2uLJe9dQXAdeG2RNODHYonfUULw4yYUSQe2FiUTdgXBx00o0jHXfmwPkZPDdQXBKjBFOkYi\nh9EkN10hbrSbFocAU6RDZHLYTPLSFfJGe2lxDDBFOkTsRlyRFK120uIgYIp0iHyUPYok2cdJ\nV1CkRmCKdIzqHqN7suSkKyhSIzBFeohqIe7xydLFxl66gjVSGzBFqgbHEomrdm3AFKkaLCuW\nPLWYYDyYIl2ABRdwUbE0QVcQXNqEIhXBopJCUiy5qOpUiX/yeoIpUhGseli0b5zYxcM6oy7h\nT15XMEUqgm0ipfYZ/eRLn/AnryuYIhXBuoF69Oh5P4o0N5gilcG6cXqf2G0nmSjS3GCKdAE2\nDNMtZRJrpLnBFKkBOGUSV+3mBlMkNVj6ZKmPSAT7AFMkLVg6mXoUiTeOycEUSQmWl/dPHr3/\npolNE/bxIaouo0jqRBDpa+vDTk0WBibs4z26LqNI6rgQSeVTM5cm7ON7lB321moKTZEagJ+m\na6IkF8RBmbGPb1GL1OhpAEVqAX6crsmwDU2aso+/ohVp4PSZIlnB+mLpUSTUCXfQFe3AyhqJ\nImnjAKx04u31cF961Q6RC3CTuACrOokiqeMBrCuW3h53wp1yD13hBswaSRkXYFWxdAOfd6o8\n9S66wg24iUcUqQtYJ9Jpp9qLqKeuWBhMkarBJpFSczybUJ66YmEwRaoHa2qkh73uO9fcmlx1\nxbpgigQAS1ftMvt+CmRdfvDVFcuCKVILcEKILLh2jue8K1YBU6QG4NStJQ+unOP57oplwBQJ\nD05O0uQv8SnneK67Yh0wRcKDjSIZHzG57op1wBQJD06aYHn3yLYciArBGjBFagA+emT4FiHd\nHM93VywDpkgtwOdJWvt3j+AhWAOmSA3BlW8tJETSrKvXhmANmCI1BNe+tZDyqNOXihOsBFOk\nlmDjivZx/xPOsBxoC8EaMEVqClZVOxJazXKgMgRrwBSpBxgukm05UBeCNWCK1AUM8Qj4sxeS\nROtjimRIOPB5vJsMOE0VQYamEq6PKZI+8cFVBtxFanlrit/HPcEUaQy4smo6z/Hwt6bwfdwV\nTJHGgGtvJc9zvAa3pvB93BVMkcaAUbeShrem8H3cFUyRBoETC3CmtLs1xe/jnmCKNAp8WoCr\n+xYh/K1pgj7uCKZIg8HVA/8Ght+a5unjHmCKNBpcO8e7njNWGgrPlGCKNBx8nuPVgVG3ppn6\nuD2YInkBw0R6AL7ab00T9nFDMEVyA04Mc8m4l9/qlKrO2MftwBTJDzjlEeLLTwpVWIk/ZR83\nA1Mkv2BZcSM39CxS0VRPXeEfXCkS0zBfA//2OxAweYA3CJ/5CO9I3sCy4kYDPk/sCssQnrrC\nP5giOQYXips68JGfM9WygiiLqz5GgSmSZ3C2uKkF7/zCMgR/YlABpkgRwKX7RS34zj8vQ/An\nBuVgihQCnPAI+uUn6RsSf2JQDqZIEcH4+0XSVMNTXEli9LESTJEigp/vFw3KGdFChy0x+lgJ\npkgRwc/3iwajHfzdlsfE6GMlmCKFBD/dL5qY9PZ4KEmbVGB8KJI6y4Mf7xcPJqGEurVY7pHw\nwFH6WAWmSLHBZ5FgtyZdixX3xGh9TJFWAH+tNTRYZQOLtP9luD6mSKuAn29IDkU6/K1yzigP\nRVKH4BO4sMpmG7HKFgs8+vp77SqGOBRJHYKz4KRHPb785Gpi9yRSg4VGiqQOwXlw0iPkl0FY\nQpEoUnSwC5HONRJFokixwD5EOq/aDamRmpWLFGl+cKcaSQ0esGrXriso0gLgLqt2EcANb84U\nieB1wBSJYIIBYIpEMMEIMGskgglGgLlqRzDBnsEUiWCCAWCKRDDBADBFIphgAJgiEUwwAEyR\nCCYYAKZIBBMMAFMkggkGgCkSwQQDwBSJYIIBYIpEMMEAMEUimGAAmCIRTDAATJEIJhgApkgE\nEwwAUySCCQaAKRLBBAPAFIlgggFgikQwwQAwRSKYYACYIhFMMABMkQgmGACmSAQTDABTJIIJ\nBoApEsEEA8AUiWCCAWCKRDDBADBFIphgAJgiEUwwAEyRCCYYAKZIBBMMAFeKdJE/GjDbhi1u\nn8lbTJE+wha3z+QtpkgfYYvbZ/IWU6SPsMXtM3mLW4jEMMuFIjEMIBSJYQChSAwDCEViGEAo\nEsMAghfp5T1waLPcWhuq1cFaHK+PDS2Gi/Ry/yVCbq0N1erfJzdOi+P1saXFFCnYSf5oJkVq\nG4pkS6iTfGttnBa/HP8bp8UUSR2K1DQvby8LTEYpUrhhGa7FX82N0+LbIgNFUiXUsLw3NEyL\nWSPVNSJMXh5/cZ6XzwRqMUWqa0SUvOy/hmk170htQ5H0eTn8J0yrKVLbuBAp0PPrt9tM6S1a\nq3//EqbFfLOBYRhZKBLDAEKRGAYQisQwgFAkhgGEIjEMIBSJYQChSAwDCEViGEAoEsMAQpEY\nBhCKxDCAUKRA+ec/2/by8+N3v75v3/7ePk7evz+27ce/gxvGUKRAeTfnI+8m/fvy+dv3P/z9\nu2+jm8ZQpDj5tv337e1/H/r8uX1/+/f75+/evfq5/TW6bcuHIkXKr7///K3Pt+3X+//5/N3H\nn2//Gd2y5UORAuX7bUK3ferz+3f3SR4zMjwDcfJj+/bX378oksvwDMTJ5ypdamrHDA/PQ5xs\n2z9fSww/t++/J3ofv/v59vbfj//LDA1FipOf92ncvvz99bv/jW7b8qFIgfJj277/83uC9/FA\n9r+fv/v9h6NbxlCkuNmCfLnVEqFIEfNRLb3P9H6MbgdzD0WKmK9q6dfodjD3UKSQ+evbtv2g\nR45CkRgGEIrEMIBQJIYBhCIxDCAUiWEAoUgMAwhFYhhAKBLDAPJ/3sgAIHYZzv8AAAAASUVO\nRK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the log odds against the x (independent) predictor variable\n",
    "plot.dat <- data.frame(logodds = logodds, age = df_titanic$Age)\n",
    "ggplot(plot.dat, aes(x=age, y=logodds)) + geom_point()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Reference: https://stackoverflow.com/questions/27278745/can-the-boxtidwell-function-handle-binary-outcome-variables\n",
    "# x<-df_titanic$Age\n",
    "# y<-df_titanic$Survived\n",
    "\n",
    "# #FIT IS DONE USING THE glm FUNCTION\n",
    "# m1res <- glm(y ~ x,family=binomial(link = \"logit\"))\n",
    "# coeff1<- coefficients(summary(m1res))\n",
    "# lnx<-x*log(x)\n",
    "# m2res <- glm(y ~ x+lnx ,family=binomial(link = \"logit\"))\n",
    "# coeff2<- coefficients(summary(m2res))\n",
    "# alpha0<-1.0\n",
    "# pvalue<-coeff2[3,4]\n",
    "\n",
    "# pvalue\n",
    "# beta1<-coeff1[2,1]\n",
    "# beta2<-coeff2[3,1]\n",
    "# iter<-0\n",
    "# err<-1\n",
    "\n",
    "# while (pvalue<0.1) {\n",
    "#     alpha <-(beta2/beta1)+alpha0 \n",
    "#     err<-abs(alpha-alpha0)\n",
    "#     alpha0<-alpha\n",
    "#     mx<-x^alpha\n",
    "#     m1res <- glm(y ~ mx,family=binomial(link = \"logit\"))\n",
    "#     coeff1<- coefficients(summary(m1res))\n",
    "#     mlnx<-mx*log(x)\n",
    "#     m2res <- glm(y ~ mx+mlnx ,family=binomial(link = \"logit\"))\n",
    "#     coeff2<- coefficients(summary(m2res))\n",
    "#     pvalue<-coeff2[3,4]\n",
    "#     beta1<-coeff1[2,1]\n",
    "#     beta2<-coeff2[3,1]\n",
    "#     iter<- iter+1\n",
    "#     }\n",
    "\n",
    "# # PRINT THE POWER TO CONSOLE\n",
    "# alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 102 × 6</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>education</th><th scope=col>income</th><th scope=col>women</th><th scope=col>prestige</th><th scope=col>census</th><th scope=col>type</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;fct&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>gov.administrators</th><td>13.11</td><td>12351</td><td>11.16</td><td>68.8</td><td>1113</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>general.managers</th><td>12.26</td><td>25879</td><td> 4.02</td><td>69.1</td><td>1130</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>accountants</th><td>12.77</td><td> 9271</td><td>15.70</td><td>63.4</td><td>1171</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>purchasing.officers</th><td>11.42</td><td> 8865</td><td> 9.11</td><td>56.8</td><td>1175</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>chemists</th><td>14.62</td><td> 8403</td><td>11.68</td><td>73.5</td><td>2111</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>physicists</th><td>15.64</td><td>11030</td><td> 5.13</td><td>77.6</td><td>2113</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>biologists</th><td>15.09</td><td> 8258</td><td>25.65</td><td>72.6</td><td>2133</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>architects</th><td>15.44</td><td>14163</td><td> 2.69</td><td>78.1</td><td>2141</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>civil.engineers</th><td>14.52</td><td>11377</td><td> 1.03</td><td>73.1</td><td>2143</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>mining.engineers</th><td>14.64</td><td>11023</td><td> 0.94</td><td>68.8</td><td>2153</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>surveyors</th><td>12.39</td><td> 5902</td><td> 1.91</td><td>62.0</td><td>2161</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>draughtsmen</th><td>12.30</td><td> 7059</td><td> 7.83</td><td>60.0</td><td>2163</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>computer.programers</th><td>13.83</td><td> 8425</td><td>15.33</td><td>53.8</td><td>2183</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>economists</th><td>14.44</td><td> 8049</td><td>57.31</td><td>62.2</td><td>2311</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>psychologists</th><td>14.36</td><td> 7405</td><td>48.28</td><td>74.9</td><td>2315</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>social.workers</th><td>14.21</td><td> 6336</td><td>54.77</td><td>55.1</td><td>2331</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>lawyers</th><td>15.77</td><td>19263</td><td> 5.13</td><td>82.3</td><td>2343</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>librarians</th><td>14.15</td><td> 6112</td><td>77.10</td><td>58.1</td><td>2351</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>vocational.counsellors</th><td>15.22</td><td> 9593</td><td>34.89</td><td>58.3</td><td>2391</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>ministers</th><td>14.50</td><td> 4686</td><td> 4.14</td><td>72.8</td><td>2511</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>university.teachers</th><td>15.97</td><td>12480</td><td>19.59</td><td>84.6</td><td>2711</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>primary.school.teachers</th><td>13.62</td><td> 5648</td><td>83.78</td><td>59.6</td><td>2731</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>secondary.school.teachers</th><td>15.08</td><td> 8034</td><td>46.80</td><td>66.1</td><td>2733</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>physicians</th><td>15.96</td><td>25308</td><td>10.56</td><td>87.2</td><td>3111</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>veterinarians</th><td>15.94</td><td>14558</td><td> 4.32</td><td>66.7</td><td>3115</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>osteopaths.chiropractors</th><td>14.71</td><td>17498</td><td> 6.91</td><td>68.4</td><td>3117</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>nurses</th><td>12.46</td><td> 4614</td><td>96.12</td><td>64.7</td><td>3131</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>nursing.aides</th><td> 9.45</td><td> 3485</td><td>76.14</td><td>34.9</td><td>3135</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>physio.therapsts</th><td>13.62</td><td> 5092</td><td>82.66</td><td>72.1</td><td>3137</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>pharmacists</th><td>15.21</td><td>10432</td><td>24.71</td><td>69.3</td><td>3151</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>...</th><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr>\n",
       "\t<tr><th scope=row>canners</th><td> 7.42</td><td> 1890</td><td>72.24</td><td>23.2</td><td>8221</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>textile.weavers</th><td> 6.69</td><td> 4443</td><td>31.36</td><td>33.3</td><td>8267</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>textile.labourers</th><td> 6.74</td><td> 3485</td><td>39.48</td><td>28.8</td><td>8278</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>tool.die.makers</th><td>10.09</td><td> 8043</td><td> 1.50</td><td>42.5</td><td>8311</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>machinists</th><td> 8.81</td><td> 6686</td><td> 4.28</td><td>44.2</td><td>8313</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>sheet.metal.workers</th><td> 8.40</td><td> 6565</td><td> 2.30</td><td>35.9</td><td>8333</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>welders</th><td> 7.92</td><td> 6477</td><td> 5.17</td><td>41.8</td><td>8335</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>auto.workers</th><td> 8.43</td><td> 5811</td><td>13.62</td><td>35.9</td><td>8513</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>aircraft.workers</th><td> 8.78</td><td> 6573</td><td> 5.78</td><td>43.7</td><td>8515</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>electronic.workers</th><td> 8.76</td><td> 3942</td><td>74.54</td><td>50.8</td><td>8534</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>radio.tv.repairmen</th><td>10.29</td><td> 5449</td><td> 2.92</td><td>37.2</td><td>8537</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>sewing.mach.operators</th><td> 6.38</td><td> 2847</td><td>90.67</td><td>28.2</td><td>8563</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>auto.repairmen</th><td> 8.10</td><td> 5795</td><td> 0.81</td><td>38.1</td><td>8581</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>aircraft.repairmen</th><td>10.10</td><td> 7716</td><td> 0.78</td><td>50.3</td><td>8582</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>railway.sectionmen</th><td> 6.67</td><td> 4696</td><td> 0.00</td><td>27.3</td><td>8715</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>electrical.linemen</th><td> 9.05</td><td> 8316</td><td> 1.34</td><td>40.9</td><td>8731</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>electricians</th><td> 9.93</td><td> 7147</td><td> 0.99</td><td>50.2</td><td>8733</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>construction.foremen</th><td> 8.24</td><td> 8880</td><td> 0.65</td><td>51.1</td><td>8780</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>carpenters</th><td> 6.92</td><td> 5299</td><td> 0.56</td><td>38.9</td><td>8781</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>masons</th><td> 6.60</td><td> 5959</td><td> 0.52</td><td>36.2</td><td>8782</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>house.painters</th><td> 7.81</td><td> 4549</td><td> 2.46</td><td>29.9</td><td>8785</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>plumbers</th><td> 8.33</td><td> 6928</td><td> 0.61</td><td>42.9</td><td>8791</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>construction.labourers</th><td> 7.52</td><td> 3910</td><td> 1.09</td><td>26.5</td><td>8798</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>pilots</th><td>12.27</td><td>14032</td><td> 0.58</td><td>66.1</td><td>9111</td><td>prof</td></tr>\n",
       "\t<tr><th scope=row>train.engineers</th><td> 8.49</td><td> 8845</td><td> 0.00</td><td>48.9</td><td>9131</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>bus.drivers</th><td> 7.58</td><td> 5562</td><td> 9.47</td><td>35.9</td><td>9171</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>taxi.drivers</th><td> 7.93</td><td> 4224</td><td> 3.59</td><td>25.1</td><td>9173</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>longshoremen</th><td> 8.37</td><td> 4753</td><td> 0.00</td><td>26.1</td><td>9313</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>typesetters</th><td>10.00</td><td> 6462</td><td>13.58</td><td>42.2</td><td>9511</td><td>bc  </td></tr>\n",
       "\t<tr><th scope=row>bookbinders</th><td> 8.55</td><td> 3617</td><td>70.87</td><td>35.2</td><td>9517</td><td>bc  </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 102 × 6\n",
       "\\begin{tabular}{r|llllll}\n",
       "  & education & income & women & prestige & census & type\\\\\n",
       "  & <dbl> & <int> & <dbl> & <dbl> & <int> & <fct>\\\\\n",
       "\\hline\n",
       "\tgov.administrators & 13.11 & 12351 & 11.16 & 68.8 & 1113 & prof\\\\\n",
       "\tgeneral.managers & 12.26 & 25879 &  4.02 & 69.1 & 1130 & prof\\\\\n",
       "\taccountants & 12.77 &  9271 & 15.70 & 63.4 & 1171 & prof\\\\\n",
       "\tpurchasing.officers & 11.42 &  8865 &  9.11 & 56.8 & 1175 & prof\\\\\n",
       "\tchemists & 14.62 &  8403 & 11.68 & 73.5 & 2111 & prof\\\\\n",
       "\tphysicists & 15.64 & 11030 &  5.13 & 77.6 & 2113 & prof\\\\\n",
       "\tbiologists & 15.09 &  8258 & 25.65 & 72.6 & 2133 & prof\\\\\n",
       "\tarchitects & 15.44 & 14163 &  2.69 & 78.1 & 2141 & prof\\\\\n",
       "\tcivil.engineers & 14.52 & 11377 &  1.03 & 73.1 & 2143 & prof\\\\\n",
       "\tmining.engineers & 14.64 & 11023 &  0.94 & 68.8 & 2153 & prof\\\\\n",
       "\tsurveyors & 12.39 &  5902 &  1.91 & 62.0 & 2161 & prof\\\\\n",
       "\tdraughtsmen & 12.30 &  7059 &  7.83 & 60.0 & 2163 & prof\\\\\n",
       "\tcomputer.programers & 13.83 &  8425 & 15.33 & 53.8 & 2183 & prof\\\\\n",
       "\teconomists & 14.44 &  8049 & 57.31 & 62.2 & 2311 & prof\\\\\n",
       "\tpsychologists & 14.36 &  7405 & 48.28 & 74.9 & 2315 & prof\\\\\n",
       "\tsocial.workers & 14.21 &  6336 & 54.77 & 55.1 & 2331 & prof\\\\\n",
       "\tlawyers & 15.77 & 19263 &  5.13 & 82.3 & 2343 & prof\\\\\n",
       "\tlibrarians & 14.15 &  6112 & 77.10 & 58.1 & 2351 & prof\\\\\n",
       "\tvocational.counsellors & 15.22 &  9593 & 34.89 & 58.3 & 2391 & prof\\\\\n",
       "\tministers & 14.50 &  4686 &  4.14 & 72.8 & 2511 & prof\\\\\n",
       "\tuniversity.teachers & 15.97 & 12480 & 19.59 & 84.6 & 2711 & prof\\\\\n",
       "\tprimary.school.teachers & 13.62 &  5648 & 83.78 & 59.6 & 2731 & prof\\\\\n",
       "\tsecondary.school.teachers & 15.08 &  8034 & 46.80 & 66.1 & 2733 & prof\\\\\n",
       "\tphysicians & 15.96 & 25308 & 10.56 & 87.2 & 3111 & prof\\\\\n",
       "\tveterinarians & 15.94 & 14558 &  4.32 & 66.7 & 3115 & prof\\\\\n",
       "\tosteopaths.chiropractors & 14.71 & 17498 &  6.91 & 68.4 & 3117 & prof\\\\\n",
       "\tnurses & 12.46 &  4614 & 96.12 & 64.7 & 3131 & prof\\\\\n",
       "\tnursing.aides &  9.45 &  3485 & 76.14 & 34.9 & 3135 & bc  \\\\\n",
       "\tphysio.therapsts & 13.62 &  5092 & 82.66 & 72.1 & 3137 & prof\\\\\n",
       "\tpharmacists & 15.21 & 10432 & 24.71 & 69.3 & 3151 & prof\\\\\n",
       "\t... & ... & ... & ... & ... & ... & ...\\\\\n",
       "\tcanners &  7.42 &  1890 & 72.24 & 23.2 & 8221 & bc  \\\\\n",
       "\ttextile.weavers &  6.69 &  4443 & 31.36 & 33.3 & 8267 & bc  \\\\\n",
       "\ttextile.labourers &  6.74 &  3485 & 39.48 & 28.8 & 8278 & bc  \\\\\n",
       "\ttool.die.makers & 10.09 &  8043 &  1.50 & 42.5 & 8311 & bc  \\\\\n",
       "\tmachinists &  8.81 &  6686 &  4.28 & 44.2 & 8313 & bc  \\\\\n",
       "\tsheet.metal.workers &  8.40 &  6565 &  2.30 & 35.9 & 8333 & bc  \\\\\n",
       "\twelders &  7.92 &  6477 &  5.17 & 41.8 & 8335 & bc  \\\\\n",
       "\tauto.workers &  8.43 &  5811 & 13.62 & 35.9 & 8513 & bc  \\\\\n",
       "\taircraft.workers &  8.78 &  6573 &  5.78 & 43.7 & 8515 & bc  \\\\\n",
       "\telectronic.workers &  8.76 &  3942 & 74.54 & 50.8 & 8534 & bc  \\\\\n",
       "\tradio.tv.repairmen & 10.29 &  5449 &  2.92 & 37.2 & 8537 & bc  \\\\\n",
       "\tsewing.mach.operators &  6.38 &  2847 & 90.67 & 28.2 & 8563 & bc  \\\\\n",
       "\tauto.repairmen &  8.10 &  5795 &  0.81 & 38.1 & 8581 & bc  \\\\\n",
       "\taircraft.repairmen & 10.10 &  7716 &  0.78 & 50.3 & 8582 & bc  \\\\\n",
       "\trailway.sectionmen &  6.67 &  4696 &  0.00 & 27.3 & 8715 & bc  \\\\\n",
       "\telectrical.linemen &  9.05 &  8316 &  1.34 & 40.9 & 8731 & bc  \\\\\n",
       "\telectricians &  9.93 &  7147 &  0.99 & 50.2 & 8733 & bc  \\\\\n",
       "\tconstruction.foremen &  8.24 &  8880 &  0.65 & 51.1 & 8780 & bc  \\\\\n",
       "\tcarpenters &  6.92 &  5299 &  0.56 & 38.9 & 8781 & bc  \\\\\n",
       "\tmasons &  6.60 &  5959 &  0.52 & 36.2 & 8782 & bc  \\\\\n",
       "\thouse.painters &  7.81 &  4549 &  2.46 & 29.9 & 8785 & bc  \\\\\n",
       "\tplumbers &  8.33 &  6928 &  0.61 & 42.9 & 8791 & bc  \\\\\n",
       "\tconstruction.labourers &  7.52 &  3910 &  1.09 & 26.5 & 8798 & bc  \\\\\n",
       "\tpilots & 12.27 & 14032 &  0.58 & 66.1 & 9111 & prof\\\\\n",
       "\ttrain.engineers &  8.49 &  8845 &  0.00 & 48.9 & 9131 & bc  \\\\\n",
       "\tbus.drivers &  7.58 &  5562 &  9.47 & 35.9 & 9171 & bc  \\\\\n",
       "\ttaxi.drivers &  7.93 &  4224 &  3.59 & 25.1 & 9173 & bc  \\\\\n",
       "\tlongshoremen &  8.37 &  4753 &  0.00 & 26.1 & 9313 & bc  \\\\\n",
       "\ttypesetters & 10.00 &  6462 & 13.58 & 42.2 & 9511 & bc  \\\\\n",
       "\tbookbinders &  8.55 &  3617 & 70.87 & 35.2 & 9517 & bc  \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 102 × 6\n",
       "\n",
       "| <!--/--> | education &lt;dbl&gt; | income &lt;int&gt; | women &lt;dbl&gt; | prestige &lt;dbl&gt; | census &lt;int&gt; | type &lt;fct&gt; |\n",
       "|---|---|---|---|---|---|---|\n",
       "| gov.administrators | 13.11 | 12351 | 11.16 | 68.8 | 1113 | prof |\n",
       "| general.managers | 12.26 | 25879 |  4.02 | 69.1 | 1130 | prof |\n",
       "| accountants | 12.77 |  9271 | 15.70 | 63.4 | 1171 | prof |\n",
       "| purchasing.officers | 11.42 |  8865 |  9.11 | 56.8 | 1175 | prof |\n",
       "| chemists | 14.62 |  8403 | 11.68 | 73.5 | 2111 | prof |\n",
       "| physicists | 15.64 | 11030 |  5.13 | 77.6 | 2113 | prof |\n",
       "| biologists | 15.09 |  8258 | 25.65 | 72.6 | 2133 | prof |\n",
       "| architects | 15.44 | 14163 |  2.69 | 78.1 | 2141 | prof |\n",
       "| civil.engineers | 14.52 | 11377 |  1.03 | 73.1 | 2143 | prof |\n",
       "| mining.engineers | 14.64 | 11023 |  0.94 | 68.8 | 2153 | prof |\n",
       "| surveyors | 12.39 |  5902 |  1.91 | 62.0 | 2161 | prof |\n",
       "| draughtsmen | 12.30 |  7059 |  7.83 | 60.0 | 2163 | prof |\n",
       "| computer.programers | 13.83 |  8425 | 15.33 | 53.8 | 2183 | prof |\n",
       "| economists | 14.44 |  8049 | 57.31 | 62.2 | 2311 | prof |\n",
       "| psychologists | 14.36 |  7405 | 48.28 | 74.9 | 2315 | prof |\n",
       "| social.workers | 14.21 |  6336 | 54.77 | 55.1 | 2331 | prof |\n",
       "| lawyers | 15.77 | 19263 |  5.13 | 82.3 | 2343 | prof |\n",
       "| librarians | 14.15 |  6112 | 77.10 | 58.1 | 2351 | prof |\n",
       "| vocational.counsellors | 15.22 |  9593 | 34.89 | 58.3 | 2391 | prof |\n",
       "| ministers | 14.50 |  4686 |  4.14 | 72.8 | 2511 | prof |\n",
       "| university.teachers | 15.97 | 12480 | 19.59 | 84.6 | 2711 | prof |\n",
       "| primary.school.teachers | 13.62 |  5648 | 83.78 | 59.6 | 2731 | prof |\n",
       "| secondary.school.teachers | 15.08 |  8034 | 46.80 | 66.1 | 2733 | prof |\n",
       "| physicians | 15.96 | 25308 | 10.56 | 87.2 | 3111 | prof |\n",
       "| veterinarians | 15.94 | 14558 |  4.32 | 66.7 | 3115 | prof |\n",
       "| osteopaths.chiropractors | 14.71 | 17498 |  6.91 | 68.4 | 3117 | prof |\n",
       "| nurses | 12.46 |  4614 | 96.12 | 64.7 | 3131 | prof |\n",
       "| nursing.aides |  9.45 |  3485 | 76.14 | 34.9 | 3135 | bc   |\n",
       "| physio.therapsts | 13.62 |  5092 | 82.66 | 72.1 | 3137 | prof |\n",
       "| pharmacists | 15.21 | 10432 | 24.71 | 69.3 | 3151 | prof |\n",
       "| ... | ... | ... | ... | ... | ... | ... |\n",
       "| canners |  7.42 |  1890 | 72.24 | 23.2 | 8221 | bc   |\n",
       "| textile.weavers |  6.69 |  4443 | 31.36 | 33.3 | 8267 | bc   |\n",
       "| textile.labourers |  6.74 |  3485 | 39.48 | 28.8 | 8278 | bc   |\n",
       "| tool.die.makers | 10.09 |  8043 |  1.50 | 42.5 | 8311 | bc   |\n",
       "| machinists |  8.81 |  6686 |  4.28 | 44.2 | 8313 | bc   |\n",
       "| sheet.metal.workers |  8.40 |  6565 |  2.30 | 35.9 | 8333 | bc   |\n",
       "| welders |  7.92 |  6477 |  5.17 | 41.8 | 8335 | bc   |\n",
       "| auto.workers |  8.43 |  5811 | 13.62 | 35.9 | 8513 | bc   |\n",
       "| aircraft.workers |  8.78 |  6573 |  5.78 | 43.7 | 8515 | bc   |\n",
       "| electronic.workers |  8.76 |  3942 | 74.54 | 50.8 | 8534 | bc   |\n",
       "| radio.tv.repairmen | 10.29 |  5449 |  2.92 | 37.2 | 8537 | bc   |\n",
       "| sewing.mach.operators |  6.38 |  2847 | 90.67 | 28.2 | 8563 | bc   |\n",
       "| auto.repairmen |  8.10 |  5795 |  0.81 | 38.1 | 8581 | bc   |\n",
       "| aircraft.repairmen | 10.10 |  7716 |  0.78 | 50.3 | 8582 | bc   |\n",
       "| railway.sectionmen |  6.67 |  4696 |  0.00 | 27.3 | 8715 | bc   |\n",
       "| electrical.linemen |  9.05 |  8316 |  1.34 | 40.9 | 8731 | bc   |\n",
       "| electricians |  9.93 |  7147 |  0.99 | 50.2 | 8733 | bc   |\n",
       "| construction.foremen |  8.24 |  8880 |  0.65 | 51.1 | 8780 | bc   |\n",
       "| carpenters |  6.92 |  5299 |  0.56 | 38.9 | 8781 | bc   |\n",
       "| masons |  6.60 |  5959 |  0.52 | 36.2 | 8782 | bc   |\n",
       "| house.painters |  7.81 |  4549 |  2.46 | 29.9 | 8785 | bc   |\n",
       "| plumbers |  8.33 |  6928 |  0.61 | 42.9 | 8791 | bc   |\n",
       "| construction.labourers |  7.52 |  3910 |  1.09 | 26.5 | 8798 | bc   |\n",
       "| pilots | 12.27 | 14032 |  0.58 | 66.1 | 9111 | prof |\n",
       "| train.engineers |  8.49 |  8845 |  0.00 | 48.9 | 9131 | bc   |\n",
       "| bus.drivers |  7.58 |  5562 |  9.47 | 35.9 | 9171 | bc   |\n",
       "| taxi.drivers |  7.93 |  4224 |  3.59 | 25.1 | 9173 | bc   |\n",
       "| longshoremen |  8.37 |  4753 |  0.00 | 26.1 | 9313 | bc   |\n",
       "| typesetters | 10.00 |  6462 | 13.58 | 42.2 | 9511 | bc   |\n",
       "| bookbinders |  8.55 |  3617 | 70.87 | 35.2 | 9517 | bc   |\n",
       "\n"
      ],
      "text/plain": [
       "                          education income women prestige census type\n",
       "gov.administrators        13.11     12351  11.16 68.8     1113   prof\n",
       "general.managers          12.26     25879   4.02 69.1     1130   prof\n",
       "accountants               12.77      9271  15.70 63.4     1171   prof\n",
       "purchasing.officers       11.42      8865   9.11 56.8     1175   prof\n",
       "chemists                  14.62      8403  11.68 73.5     2111   prof\n",
       "physicists                15.64     11030   5.13 77.6     2113   prof\n",
       "biologists                15.09      8258  25.65 72.6     2133   prof\n",
       "architects                15.44     14163   2.69 78.1     2141   prof\n",
       "civil.engineers           14.52     11377   1.03 73.1     2143   prof\n",
       "mining.engineers          14.64     11023   0.94 68.8     2153   prof\n",
       "surveyors                 12.39      5902   1.91 62.0     2161   prof\n",
       "draughtsmen               12.30      7059   7.83 60.0     2163   prof\n",
       "computer.programers       13.83      8425  15.33 53.8     2183   prof\n",
       "economists                14.44      8049  57.31 62.2     2311   prof\n",
       "psychologists             14.36      7405  48.28 74.9     2315   prof\n",
       "social.workers            14.21      6336  54.77 55.1     2331   prof\n",
       "lawyers                   15.77     19263   5.13 82.3     2343   prof\n",
       "librarians                14.15      6112  77.10 58.1     2351   prof\n",
       "vocational.counsellors    15.22      9593  34.89 58.3     2391   prof\n",
       "ministers                 14.50      4686   4.14 72.8     2511   prof\n",
       "university.teachers       15.97     12480  19.59 84.6     2711   prof\n",
       "primary.school.teachers   13.62      5648  83.78 59.6     2731   prof\n",
       "secondary.school.teachers 15.08      8034  46.80 66.1     2733   prof\n",
       "physicians                15.96     25308  10.56 87.2     3111   prof\n",
       "veterinarians             15.94     14558   4.32 66.7     3115   prof\n",
       "osteopaths.chiropractors  14.71     17498   6.91 68.4     3117   prof\n",
       "nurses                    12.46      4614  96.12 64.7     3131   prof\n",
       "nursing.aides              9.45      3485  76.14 34.9     3135   bc  \n",
       "physio.therapsts          13.62      5092  82.66 72.1     3137   prof\n",
       "pharmacists               15.21     10432  24.71 69.3     3151   prof\n",
       "...                       ...       ...    ...   ...      ...    ... \n",
       "canners                    7.42      1890  72.24 23.2     8221   bc  \n",
       "textile.weavers            6.69      4443  31.36 33.3     8267   bc  \n",
       "textile.labourers          6.74      3485  39.48 28.8     8278   bc  \n",
       "tool.die.makers           10.09      8043   1.50 42.5     8311   bc  \n",
       "machinists                 8.81      6686   4.28 44.2     8313   bc  \n",
       "sheet.metal.workers        8.40      6565   2.30 35.9     8333   bc  \n",
       "welders                    7.92      6477   5.17 41.8     8335   bc  \n",
       "auto.workers               8.43      5811  13.62 35.9     8513   bc  \n",
       "aircraft.workers           8.78      6573   5.78 43.7     8515   bc  \n",
       "electronic.workers         8.76      3942  74.54 50.8     8534   bc  \n",
       "radio.tv.repairmen        10.29      5449   2.92 37.2     8537   bc  \n",
       "sewing.mach.operators      6.38      2847  90.67 28.2     8563   bc  \n",
       "auto.repairmen             8.10      5795   0.81 38.1     8581   bc  \n",
       "aircraft.repairmen        10.10      7716   0.78 50.3     8582   bc  \n",
       "railway.sectionmen         6.67      4696   0.00 27.3     8715   bc  \n",
       "electrical.linemen         9.05      8316   1.34 40.9     8731   bc  \n",
       "electricians               9.93      7147   0.99 50.2     8733   bc  \n",
       "construction.foremen       8.24      8880   0.65 51.1     8780   bc  \n",
       "carpenters                 6.92      5299   0.56 38.9     8781   bc  \n",
       "masons                     6.60      5959   0.52 36.2     8782   bc  \n",
       "house.painters             7.81      4549   2.46 29.9     8785   bc  \n",
       "plumbers                   8.33      6928   0.61 42.9     8791   bc  \n",
       "construction.labourers     7.52      3910   1.09 26.5     8798   bc  \n",
       "pilots                    12.27     14032   0.58 66.1     9111   prof\n",
       "train.engineers            8.49      8845   0.00 48.9     9131   bc  \n",
       "bus.drivers                7.58      5562   9.47 35.9     9171   bc  \n",
       "taxi.drivers               7.93      4224   3.59 25.1     9173   bc  \n",
       "longshoremen               8.37      4753   0.00 26.1     9313   bc  \n",
       "typesetters               10.00      6462  13.58 42.2     9511   bc  \n",
       "bookbinders                8.55      3617  70.87 35.2     9517   bc  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Standard example for Box-Tidwell (on Prestige dataset)\n",
    "# Prestige"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "          MLE of lambda Score Statistic (z)  Pr(>|z|)    \n",
       "income         -0.34763             -4.4824 7.381e-06 ***\n",
       "education       1.25383              0.2170    0.8282    \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "iterations =  8 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# boxTidwell(prestige ~ income + education, ~ type + poly(women, 2), data=Prestige)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}