Logistic Regression using R - How to handle class bias?

Doing this blog as a revision into some primitive topics in Machine Learning.

Since, almost all my current work is done using R, I will implement logistic regression using this language.

We will be using adult data set. After all is set, now we dive straight into the code. Further explanations shall not be provided as most of the code is self-explanatory in nature.

> idata <- read.csv("adult_csv.csv")

> summary(idata)

      age                   workclass         fnlwgt               education     education.num  
 Min.   :0.000   Private         :33906   Min.   :  12285   HS-grad     :15784   Min.   : 1.00  
 1st Qu.:1.000   Self-emp-not-inc: 3862   1st Qu.: 117551   Some-college:10878   1st Qu.: 9.00  
 Median :2.000   Local-gov       : 3136   Median : 178145   Bachelors   : 8025   Median :10.00  
 Mean   :1.771                   : 2799   Mean   : 189664   Masters     : 2657   Mean   :10.08  
 3rd Qu.:3.000   State-gov       : 1981   3rd Qu.: 237642   Assoc-voc   : 2061   3rd Qu.:12.00  
 Max.   :4.000   Self-emp-inc    : 1695   Max.   :1490400   11th        : 1812   Max.   :16.00  
                 (Other)         : 1463                     (Other)     : 7625                  
               marital.status            occupation            relationship                   race      
 Divorced             : 6633   Prof-specialty : 6172   Husband       :19716   Amer-Indian-Eskimo:  470  
 Married-AF-spouse    :   37   Craft-repair   : 6112   Not-in-family :12583   Asian-Pac-Islander: 1519  
 Married-civ-spouse   :22379   Exec-managerial: 6086   Other-relative: 1506   Black             : 4685  
 Married-spouse-absent:  628   Adm-clerical   : 5611   Own-child     : 7581   Other             :  406  
 Never-married        :16117   Sales          : 5504   Unmarried     : 5125   White             :41762  
 Separated            : 1530   Other-service  : 4923   Wife          : 2331                             
 Widowed              : 1518   (Other)        :14434                                                    
     sex         capitalgain      capitalloss      hoursperweek         native.country    class      
 Female:16192   Min.   :0.0000   Min.   :0.0000   Min.   :0.000   United-States:43832   <=50K:37155  
 Male  :32650   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:2.000   Mexico       :  951   >50K :11687  
                Median :0.0000   Median :0.0000   Median :2.000                :  857                
                Mean   :0.2003   Mean   :0.1149   Mean   :1.951   Philippines  :  295                
                3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:2.000   Germany      :  206                
                Max.   :4.0000   Max.   :4.0000   Max.   :4.000   Puerto-Rico  :  184                
                                                                  (Other)      : 2517  

> head(idata)

  age        workclass fnlwgt education education.num     marital.status        occupation  relationship  race
1   2        State-gov  77516 Bachelors            13      Never-married      Adm-clerical Not-in-family White
2   3 Self-emp-not-inc  83311 Bachelors            13 Married-civ-spouse   Exec-managerial       Husband White
3   2          Private 215646   HS-grad             9           Divorced Handlers-cleaners Not-in-family White
4   3          Private 234721      11th             7 Married-civ-spouse Handlers-cleaners       Husband Black
5   1          Private 338409 Bachelors            13 Married-civ-spouse    Prof-specialty          Wife Black
6   2          Private 284582   Masters            14 Married-civ-spouse   Exec-managerial          Wife White
     sex capitalgain capitalloss hoursperweek native.country class
1   Male           1           0            2  United-States <=50K
2   Male           0           0            0  United-States <=50K
3   Male           0           0            2  United-States <=50K
4   Male           0           0            2  United-States <=50K
5 Female           0           0            2           Cuba <=50K
6 Female           0           0            2  United-States <=50K

Checking class bias now, on column class,

> table(idata$class)

<=50K  >50K 
37155 11687

There is a class bias, a condition observed when the proportion of one event is much smaller than proportion of other event.
Thus, we must sample the observations in approximately equal proportions to get better models.

CREATING TRAINING AND TESTING SAMPLES

#training data
input_ones <- idata[which(idata$class == ">50K"), ]  # all <=50k values
input_zeros <- idata[which(idata$class == "<=50K"), ]  # all <=50k values
set.seed(100)  # for repeatability of samples
input_ones_training_rows <- sample(1:nrow(input_ones), 0.7*nrow(input_ones))  # >50k values for training
input_zeros_training_rows <- sample(1:nrow(input_zeros), 0.7*nrow(input_ones))  # <=50k values training. 
training_ones <- input_ones[input_ones_training_rows, ]  
training_zeros <- input_zeros[input_zeros_training_rows, ]
trainingData <- rbind(training_ones, training_zeros)  # row bind  

# Create Test Data
test_ones <- input_ones[-input_ones_training_rows, ]
test_zeros <- input_zeros[-input_zeros_training_rows, ]
testData <- rbind(test_ones, test_zeros)  # row bind 

One simple optional step is to create Weight of Evidence (WOE) for categorical variables.
This is a very simple, yet useful method - Read more here.

idata_woe <- idata
str(idata_woe)

table(idata_woe$class)

library(car)
library(dplyr)

idata_woe$class <- dplyr::recode(idata_woe$class,"<=50K" = 0,">50K" = 1)
#idata_woe$class <- car::recode(idata_woe$class,"<=50K" = 0,">50K" = 1) - this is not working
idata_woe

i <- 10
q <- quantile(idata_woe$class,
              probs = c(1:(i-1)/i),
              na.rm = TRUE,
              type=3)
cuts <- unique(q)


WOE_Class <- table(findInterval(idata_woe$class,vec = cuts,rightmost.closed = FALSE),idata_woe$class)
WOE_Class <- as.data.frame.matrix(WOE_Class)
WOE_Class$`0`<-rowSums(WOE_Class)
WOE_Class$WOE <- log((WOE_Class$`1`*sum(WOE_Class$`0`))/(WOE_Class$`0`*sum(WOE_Class$`1`)))

> WOE_Class
      0     1      WOE
1 37155     0     -Inf
2 11687 11687 1.430113

Now that we have already calculated WOE, it is only logical that we calculate the Information Value (IV) also,

library(scorecard)

iv(idata_woe, y="class")

> iv(idata_woe, y="class")
          variable info_value
 1:   relationship 1.52771659
 2: marital.status 1.34867700
 3:            age 0.89506214
 4:    capitalgain 0.89314987
 5:     occupation 0.76787908
 6:      education 0.73690518
 7:  education.num 0.73690518
 8:   hoursperweek 0.40824961
 9:         fnlwgt 0.33755308
10:            sex 0.30052143
11:      workclass 0.17056215
12:    capitalloss 0.11592250
13: native.country 0.07509682
14:           race 0.06784214

Now, we have the information value (meaning the strength of relation between “class” and other parameters.
We now start building the logistic regression model,

logitMod <- glm(class ~ relationship + age + capitalgain + occupation + education.num, data=trainingData, family=binomial(link="logit"))
predicted <- plogis(predict(logitMod, testData))  # predicted scores
# or
predicted <- predict(logitMod, testData, type="response")  # predicted scores
summary(logitMod)

Result is as follows,

> summary(logitMod)

Call:
glm(formula = class ~ relationship + age + capitalgain + occupation + 
    education.num, family = binomial(link = "logit"), data = trainingData)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.3583  -0.5408   0.0015   0.6281   3.3139  

Coefficients:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)                 -4.43416    0.18263 -24.280  < 2e-16 ***
relationshipNot-in-family   -2.26671    0.05846 -38.776  < 2e-16 ***
relationshipOther-relative  -2.52271    0.19614 -12.862  < 2e-16 ***
relationshipOwn-child       -3.43974    0.14067 -24.452  < 2e-16 ***
relationshipUnmarried       -2.66316    0.09625 -27.670  < 2e-16 ***
relationshipWife             0.31766    0.09130   3.479 0.000503 ***
age                          0.34161    0.02041  16.734  < 2e-16 ***
capitalgain                  0.95230    0.03972  23.974  < 2e-16 ***
occupationAdm-clerical       1.08011    0.13682   7.895 2.91e-15 ***
occupationArmed-Forces       2.89421    1.34843   2.146 0.031844 *  
occupationCraft-repair       1.29746    0.13042   9.949  < 2e-16 ***
occupationExec-managerial    2.08246    0.13025  15.988  < 2e-16 ***
occupationFarming-fishing    0.23595    0.17661   1.336 0.181539    
occupationHandlers-cleaners  0.58579    0.18067   3.242 0.001185 ** 
occupationMachine-op-inspct  0.83323    0.14885   5.598 2.17e-08 ***
occupationOther-service      0.16679    0.16294   1.024 0.306012    
occupationPriv-house-serv   -0.02904    0.74544  -0.039 0.968925    
occupationProf-specialty     1.70451    0.13280  12.836  < 2e-16 ***
occupationProtective-serv    1.72939    0.17889   9.667  < 2e-16 ***
occupationSales              1.56500    0.13284  11.781  < 2e-16 ***
occupationTech-support       1.48003    0.16098   9.194  < 2e-16 ***
occupationTransport-moving   1.21363    0.14650   8.284  < 2e-16 ***
education.num                0.29476    0.01128  26.127  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 22680  on 16359  degrees of freedom
Residual deviance: 13115  on 16337  degrees of freedom
AIC: 13161

Number of Fisher Scoring iterations: 6

Cheers!

Markov Chains - A Mathematical Treatise

Markov Chains

We have a set of stated, $\mathcal{S}$ = {${s_1, s_2, ..., s_m}$}. A Markov process starts in one of these states and moves successively from one state to another. Each move is called a step . If the chain is currently in state $s_i$, then it moves to state $s_j$ at the next step with a probability denoted by $p_{ij}$, and this probability does not depend upon which states the chain was in before the current state.These probabilities $p_{ij}$ are called transition probabilities. The process can remain in the state it is in, and this occurs with probability $p_{ii}$ . Initially a probability is defined by specifying a particular state as the starting state.

Say, we stay in Patna (Capital of the State of Bihar, India), and at my home we usually eat either Roti ($\text{Ro}$), Rice ($\text{Ri}$), Parantha ($\text{Pa}$)) or Bread ($\text{Br}$) for dinner. There are also some days when we eat one particular dinner on two or three or more consecutive days.

Now, we define transition probabilities of using a Transition Matrix which is given as,
$$\text{T}=\begin{array}{c c} & \begin{array}{c c c} \text{Ri} & \text{Ro} & \text{Pa} & \text{Br} \\ \end{array} \\ \begin{array}{c c c} \text{Ri} \\ \text{Ro} \\ \text{Pa} \\ \text{Br}\\ \end{array} & \left[ \begin{array}{c c c} \frac{1}{2} & \frac{1}{4} & \frac{1}{4} & 0 \\ 0 & \frac{1}{4} & \frac{1}{4} & \frac{1}{2} \\ 0 & \frac{1}{2} & \frac{1}{4} & \frac{1}{4} \\ \frac{1}{2} & 0 & \frac{1}{4} & \frac{1}{4} \\ \end{array} \right] \end{array}$$
From the above matrix, it can be seen that $P_{(Rice \implies Rice)}=\frac{1}{2}$ and similarly, $P_{(Roti \implies Parantha)}=\frac{1}{4}$.

In the similar lines of the above concept, a Transition Diagram can also be drawn.

Let us dig further into the concept now. Imagine today we eat Rice at home and two days from today we eat Bread. This is written as $P^{2}_{(Rice \implies Bread)}$. We see that if we ate Rice today then the event of us eating Bread two days from today is a disjoint union of the following,

1. We ate Rice tomorrow and Bread day after.
2. We ate something else tomorrow and Bread day after.
3. We ate Bread tomorrow and Bread day after.

Therefore, we can write $P^{2}_{(Rice \implies Bread)}$ as ,
$$P^{2}_{(Rice \implies Bread)} = [P_{(Rice \implies Rice)}]\cdot[P_{(Rice \implies Bread)}]+[P_{(Rice \implies Roti)}]\cdot[P_{(Roti \implies Bread)}]+ \\ [P_{(Rice \implies Parantha)}]\cdot[P_{(Parantha \implies Bread)}]+[P_{(Rice \implies Bread)}]\cdot[P_{(Bread \implies Bread)}]$$
The has pointed towards a more generalized concept of Dot product of two vectors.

Considering there are $r$ states in the Markov chain,
$$P^{2}_{ij}=\sum_{k=1}^{r}P_{ik}P_{kj}$$
This study was for $n=2$ , meaning we are yet to define this in a more generalized manner. Obviously, looking at the present scenario of high-end computing software packages available, this attempt seems to be a waste of efforts, but then knowing the math is always good fun!

Before doing that, we will again go into the basic definition as defined above and try to make an iterative approach into the problem.

We have assumed that our meal starts at state Rice, so, for us,
$$\text{p}_0=\begin{bmatrix} 1 \\ 0 \\ 0 \\ 0\\ \end{bmatrix}$$
$\text{p}_0$ is the state of our system at $0^{th}$ state or the beginning of the system.

Similarly, if we want to find $\text{p}_1$ we will do it as,
$$\text{p}_1 = \text{T} \cdot \text{p}_0$$
In the similar lines,
$$\boxed{\text{p}_5 = \text{T} \cdot \text{p}_4 = \text{T}^{2} \cdot \text{p}_3 = \text{T}^3 \cdot \text{p}_2 = \text{T}^4 \cdot \text{p}_1 = \text{T}^{5} \cdot \text{p}_0}$$
This means, for generalization’s sake,
$$\text{p}_n = \text{T}^{n}\cdot \text{p}_0$$
Using some very basic Linear Algebra methods to compute $\text{T}^n$,

Recall : $\text{T}=\mathcal{u}\cdot D \cdot\mathcal{u}^{-1}$

Here, $D$ is a Diagonal matrix and $\mathcal{u}$ is a matrix whose columns correspond to the Eigen-vectors of $\text{T}$.

I leave the computation to your own practice.

There for,
$$\text{T}^n = \mathcal{u} \cdot D^n . \mathcal{u}^{-1}$$
Since, $D$ is a diagonal matrix, which is obviously of form,
$$D = \begin{bmatrix} a & 0 & 0 & 0 \\ 0 & b & 0 & 0 \\ 0 & 0 & c & 0 \\ 0 & 0 & 0 & d \end{bmatrix}$$
Therefore, $D^n$ will be,
$$D = \begin{bmatrix} a^n & 0 & 0 & 0 \\ 0 & b^n & 0 & 0 \\ 0 & 0 & c^n & 0 \\ 0 & 0 & 0 & d^n \end{bmatrix}$$
Hence, we can easily compute $T^n$.

The Python implementation of this concept is relatively simple once we understand the math behind it.

All the best with that.

Cheers!

Resources:

1. Matrix Diagonalization

2. Eigen Values and Eigen Vectors

3. Markov Chains

4. Markov Chains

Study of Brachistochrone Problem.

This blog will deal with one of the most elegant topics in mathematics. The famous Brachistochrone curve. Vsuace did a great job at explaining it in this video. This looks good to me, but it doesn’t feel good; because of course the video does not contain tons of fancy mathematical vocabulary and millions of lines of $\LaTeX$.

To get started, we will just put forward a small fact that, this topic is from a larger field of mathematics called the Calculus of Variations.

Before we get into the Analysis of Brachistochrone Curve (which happens to be the heart of this post), let us brush-up on some concepts.

---

Partial Differentiation (advanced preliminary):

Let us look at how we define the derivative of a given function, say we have $g(x)$ and we have,
$$\boxed{\frac{\mathrm{d}}{\mathrm{d}x}g(x)=g'(x)=\lim_{c \to 0}\frac{f(x+ c)-f(x)}{c}}$$
Provided this limit exists.

Now, for a function of many variables it is not easy to compute the total derivative (the usual derivative). Therefore, we calculate the Partial Derivative. You might have already guessed it by now, we define partial derivative as, say we have a function $h(x,y,z)$, then,
$$\boxed{\frac {\partial}{\partial x}h(x,y,z) = \lim_{c \to 0} \frac{h(x+c, y,z)-f(x,y,z)}{c}}$$
Let us look into a small example now, say we have,
$$f(x,y,z) = x^2y^3+x+y+z+\sin(x) \cos(y)$$
Now,
$$\frac{\partial}{\partial x}f(x,y,z) = 2xy^3+1+0+0+\cos(x) \cos(y)$$
similarly,
$$\frac{\partial}{\partial y}f(x,y,z) = 3x^2y^2+0+1+0-\sin(x) \sin(y)$$
and,
$$\frac{\partial}{\partial z}f(x,y,z) = 0 +0+0+1+0$$

---

Idea of Speed and Time (basic preliminary):

Let us consider, the following scenario, (Made using Geogebra)

We have, in the figure we have denoted the path ACB as $\pi_1$ and the path ADB as $\pi_2$.
Consider that, The time taken for a body to move from A to B using path $\pi_1$ is $t_{AB_{\pi_1}}$ similarly, the time taken for a body to move from A to B using path $\pi_2$ is $t_{AB_{\pi_2}}$.
Also, consider that the distance for $\pi_1$ be $l_1$ and the distance for $\pi_2$ be $l_2$.

Then,
Speed ($\pi_1$) =$v_1$= $\frac{l_1}{t_{AB_{\pi_1}}}$ and Speed ($\pi_2$) =$v_2$ =$\frac{l_2}{t_{AB_{\pi_2}}}$.

3. Idea of Distance(basic preliminary):

Distance between two points $(a,b)$ and $(c,d)$ is given by,
$$l = \sqrt{(a-c)^2+(b-d)^2}$$
In case, the distance is measured from origin ($O(0,0)$) to some point $A(x,y)$, we get,
$$l = \sqrt{(x-0)^2+(y-0)^2} = \sqrt{x^2 + y^2}$$

---

Now that the notion is clear, we shall proceed into the analysis.
[NOTE : Some topics which are new to me (or you) will be explained as and when it is required].

---

Time needed for a body to travel from $A$ to $B$ is given by,
$t_{AB}=\frac{l}{v}$ for a linear path.

For a curve or a non-linear path, we consider piece-wise linear distance ($\mathrm{d}l$) and speed ($v$). We can define, $t_{AB}$ as,
$$\boxed{t_{AB} = \int_{A}^{B}\frac{\mathrm{d}l}{v}} \tag*{(1)}$$
Now, we must understand a fact that, all distance is composed of $x$ and $y$ coordinates. We consider the point $A$ as the origin ($O(0,0)$) and $B$ as $(x,y)$. Hence,
$$l = \sqrt{x^2 + y^2}\\ \implies \mathrm{d}l = \sqrt{\mathrm{d}x^2 + \mathrm{d}y^2}\\ \implies \mathrm{d}l = \sqrt{(1 + \frac{\mathrm{d}y^2}{\mathrm{d}x^2})\mathrm{d}x^2}$$
$$\boxed{\implies \mathrm{d}l = \sqrt{(1 + \frac{\mathrm{d}y^2}{\mathrm{d}x^2})}\mathrm{d}x} \tag*{(2)}$$

---

Also, since the translation is both in the $x$ and $y$ axes, we can say that the speed is gained by $KE$ and $PE$ equality,
$$KE = PE\\ \implies \frac{1}{2}mv^2 = mgy \\ $$
$$\implies \boxed{v = \sqrt{2gy}} \tag*{(3)}$$

---

Use $(2)$ and $(3)$ in $(1)$, we have,
$$\boxed{t_{AB} = \int_{A}^{B}\frac{\sqrt{(1 + \frac{\mathrm{d}y^2}{\mathrm{d}x^2})}\mathrm{d}x}{\sqrt{2gy}}} \tag*{(4)}$$
We write
$$\frac{\mathrm{d}y}{\mathrm{d}x}=y'$$
Now, $(4)$ becomes,
$$t_{AB} = \int_{A}^{B}\frac{\sqrt{(1 +y'^2 )}}{\sqrt{2gy}}\mathrm{d}x \tag*{(4.1)}$$
Therefore, the function is,
$$\theta = \frac{\sqrt{(1 +y'^2 )}}{\sqrt{2gy}}\\ \boxed{\theta = (1 +y'^2 )^{1/2}(2gy)^{-1/2}} \tag*{(5)}$$

---

Please notice that, the integral $t_{AB}$ can be written as,
$$t_{AB} = \int_{A}^{B} f(y, y', x) \mathrm{d}x \tag*{(6)}$$
In our case, $f(y, y', x)$ is $\theta$.

Equation$(6)$ closely follows the The Euler-Lagrange equation.

It is pretty simple to define,

---

The Euler-Lagrange equation

For any $I$ such that,
$$I = \int_{a}^{b}f(y,y',x)\mathrm{d}x$$
Where, $y' = \frac{\mathrm{d}y}{\mathrm{d}x}$, then $I$ has an stationary point, if the Euler-Lagrange Equation, given by,
$$\frac{\partial f}{\partial y}-\frac{\mathrm{d}}{\mathrm{d}x}\frac{\partial f}{\partial y'}=0$$
is satisfied.

Therefore, for our analysis,
$$\boxed{\frac{\partial \theta}{\partial y} - \frac{\mathrm{d}}{\mathrm{d}x}\frac{\partial \theta}{\partial y'}=0} \tag*{(7)}$$

---

Stationary value

This is a value at the stationary point.
A stationary point $a$ is the point at which the first derivative of a function becomes zero,
$$f'(x)|_{x=a} = 0 ; a = \text{Stationary point}$$
Side note :
Find the stationary points of $f(x) = x^3-6x^2+9x$.

---

Beltrami Identity

Our $\theta$ (from $(5)$) is such a cool expression, because $x$ does not appear in that; therefore, of course $\boxed{\frac{\partial \theta}{\partial x} = 0}$; this again leads us to another beautiful form, Beltrami Identity, which is given by,
$$\boxed{\theta - y' \frac{\partial \theta}{\partial y'}=C ; C \text{ is a constant}}.\tag*{(8)}$$

---

We now find, $\frac{\partial \theta}{\partial y'}$,
$$\frac{\partial \theta}{\partial y'}=\frac{\partial}{\partial y'}[(1 +y'^2 )^{1/2}(2gy)^{-1/2}] \\ \implies \frac{\partial \theta}{\partial y'} =(2gy)^{-1/2} \frac{1}{2}(1 +y'^2 )^{-1/2}\cdot 2y'$$
$$\therefore \boxed{\frac{\partial \theta}{\partial y'} = (2gy)^{-1/2} (1 +y'^2 )^{-1/2}\cdot y'} \tag*{(9)}$$

---

Use $\theta$ and $\frac{\partial \theta}{\partial y'}$ in equation $(8)$, we get,
$$(1 +y'^2 )^{1/2}(2gy)^{-1/2}-y' \cdot (2gy)^{-1/2} (1 +y'^2 )^{-1/2}\cdot y' = C$$
$$\implies \frac{\sqrt{(1 +y'^2 )}}{\sqrt{2gy}}-\frac{y'^2}{\sqrt{2gy}\sqrt{1+y'^2}} = C$$
$$\implies \frac{1+y'^2-y'^2}{\sqrt{2gy}\sqrt{1+y'^2}}=C$$
$$\implies \boxed{\frac{1}{\sqrt{2gy}\sqrt{1+y'^2}}=C} \tag*{(10)}$$
Rearranging a little gives,
$$\sqrt{1+y'^2}=\frac{1}{\sqrt{2gy}C}$$
$$\implies \boxed{(1+\frac{\mathrm{d}y}{\mathrm{d}x})y = \frac{1}{2gC^2} = M^2} \tag*{(11)}$$
This is the equation of the cycloid as per [4], I am yet to figure out how this happened.
The solution of $(11)$ can be found using a parametric equation,
$$\boxed{x(t) = \frac{1}{2}M^2(t-\sin t)\\ y(t)=\frac{1}{2}M^2(1-\cos t)} \tag*{(12)}$$
To derive the above equation, please refer Math-stackexchange|Derive the parametric equations of a cycloid.

I have even plotted the solutions at $(12)$, we can see that it is a cycloid,

If you wish to play with the visualisation please go to Cycloid | Parametric @ Desmos by Pragyaditya

Cheers!

---

References:
1. Stationary Points
2. Euler Lagrange Equation
3. Derivation of Beltrami Identity
4. Introduction to calculus of variations
5. Brachistochrone @ Wolfram