Logistic Regression in R

Use case machine learning with R

Note:

-An R Notebook is an R Markdown document with chunks that can be executed independently and interactively, with output visible immediately beneath the input.

-Notebook output are available as HTML, PDF, Word, or Latex.

-This Notebook as HTML is preferably open with Google Chrome.

-R-Code can be extracted as Rmd file under the button “Code” in the notebook.

-This Notebook using iterative development. It means the process starts with a simple implementation of a small set of idea requirements and iteratively enhances the evolving versions until the complete version is implemented and perfect.

What is Logistic Regression ?

• Regression analysis when the dependent variable is dichotomous (binary) such as Gender (male or female)
• Like all regression analyses, the logistic regression is a predictive analysis.
• The relationship between one dependent binary variable and one or more nominal, ordinal, interval or ratio-level independent variables.

Note: We don’t use Linear Regression for binary classification because its linear function results in probabilities outside [0,1] interval, thereby making them invalid predictions.

Type of questions that a binary logistic regression can examine:

• How does the probability of getting lung cancer (yes vs. no) change for every additional pound a person is overweight and for every pack of cigarettes smoked per day?
• Do body weight, calorie intake, fat intake, and age have an influence on the probability of having a heart attack (yes vs. no)?

Types of Logistic Regression

1. Multinomial Logistic Regression:

Let’s say our target variable has K = 4 classes. This technique handles the multi-class problem by fitting K-1 independent binary logistic classifier model. For doing this, it randomly chooses one target class as the reference class and fits K-1 regression models that compare each of the remaining classes to the reference class.

Due to its restrictive nature, it isn’t used widely because it does not scale very well in the presence of a large number of target classes. In addition, since it builds K - 1 models, we would require a much larger data set to achieve reasonable accuracy.

2. Ordinal Logistic Regression:

This technique is used when the target variable is ordinal in nature. Let’s say, we want to predict years of work experience (1,2,3,4,5, etc). So, there exists an order in the value, i.e., 5>4>3>2>1. Unlike a multinomial model, when we train K -1 models, Ordinal Logistic Regression builds a single model with multiple threshold values.

If we have K classes, the model will require K -1 threshold or cutoff points. Also, it makes an imperative assumption of proportional odds. The assumption says that on a logit (S shape) scale, all of the thresholds lie on a straight line.

Note: Logistic Regression is not a great choice to solve multi-class problems. But, it’s good to be aware of its types. In this tutorial we’ll focus on Logistic Regression for binary classification task.

How does it work

Logistic Regression assumes that the dependent (or response) variable follows a binomial distribution with the following characteristics:

• There must be a fixed number of trials denoted by n, i.e. in the data set, there must be a fixed number of rows.
• Each trial can have only two outcomes; i.e., the response variable can have only two unique categories.
• The outcome of each trial must be independent of each other; i.e., the unique levels of the response variable must be independent of each other.
• The probability of success (p) and failure (q) should be the same for each trial.

Major assumptions

• The dependent variable should be dichotomous in nature (e.g., presence vs. absent)
• There should be no outliers in the data
• There should be no high correlations (multicollinearity) among the predictors
• Mathematically, logistic regression estimates a multiple linear regression function defined as

Overfitting. When selecting the model for the logistic regression analysis, another important consideration is the model fit. Adding independent variables to a logistic regression model will always increase the amount of variance explained in the log odds (typically expressed as R²). However, adding more and more variables to the model can result in overfitting, which reduces the generalizability of the model beyond the data on which the model is fit.

Reporting the R2. Numerous pseudo-R2 values have been developed for binary logistic regression. These should be interpreted with extreme caution as they have many computational issues which cause them to be artificially high or low. A better approach is to present any of the goodness of fit tests available; Hosmer-Lemeshow is a commonly used measure of goodness of fit based on the Chi-square test.

Evaluate model and accuracy

1. Akaike Information Criteria (AIC) (or BIC)
2. Null Deviance and Residual Deviance

Use case:

Predict if an individual will earn more than $50K using logistic regression based on demographic variables data from https://github.com/itsmecevi/adult-csv

1. Import the data

inputData <- read.csv("adult.csv")
head(inputData)

A list of all variable:

names(inputData)

 [1] "AGE"           "WORKCLASS"     "FNLWGT"        "EDUCATION"     "EDUCATIONNUM" 
 [6] "MARITALSTATUS" "OCCUPATION"    "RELATIONSHIP"  "RACE"          "SEX"          
[11] "CAPITALGAIN"   "CAPITALLOSS"   "HOURSPERWEEK"  "NATIVECOUNTRY" "ABOVE50K"     

2. Check for class bias

Ideally, the proportion of events (Yes) and non-events (No) in the Y variable should approximately be the same. Lets check the proportion of classes in the dependent variable ABOVE50K.

table(inputData$ABOVE50K)

    0     1 
24720  7841 

Clearly, there is a class bias (the proportion of events is much smaller than proportion of non-events). We need sample the observations in approximately equal proportions to get better models.

3. Create training and test samples

# Create Training Data
input_ones <- inputData[which(inputData$ABOVE50K == 1), ]  # all 1's
input_zeros <- inputData[which(inputData$ABOVE50K == 0), ]  # all 0's
set.seed(100)  # for repeatability of samples
input_ones_training_rows <- sample(1:nrow(input_ones), 0.7*nrow(input_ones))  # 1's for training
input_zeros_training_rows <- sample(1:nrow(input_zeros), 0.7*nrow(input_ones))  # 0's for training. Pick as many 0's as 1's
training_ones <- input_ones[input_ones_training_rows, ]  
training_zeros <- input_zeros[input_zeros_training_rows, ]
trainingData <- rbind(training_ones, training_zeros)  # row bind the 1's and 0's 
# 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 the 1's and 0's 

4. Compute information value to find out important variables

The smbinning::smbinning function converts a continuous variable into a categorical variable using recursive partitioning. We will first convert them to categorical variables and then, capture the information values for all variables in iv_df

#install.packages("simbinning")
#library(smbinning)
# segregate continuous and factor variables
factor_vars <- c ("WORKCLASS", "EDUCATION", "MARITALSTATUS", "OCCUPATION", "RELATIONSHIP", "RACE", "SEX", "NATIVECOUNTRY")
continuous_vars <- c("AGE", "FNLWGT","EDUCATIONNUM", "HOURSPERWEEK", "CAPITALGAIN", "CAPITALLOSS")
iv_df <- data.frame(VARS=c(factor_vars, continuous_vars), IV=numeric(14))  # init for IV results
# compute IV for categoricals
for(factor_var in factor_vars){
  smb <- smbinning.factor(trainingData, y="ABOVE50K", x=factor_var)  # WOE table
  if(class(smb) != "character"){ # heck if some error occured
    iv_df[iv_df$VARS == factor_var, "IV"] <- smb$iv
  }
}

NAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercion

# compute IV for continuous vars
for(continuous_var in continuous_vars){
  smb <- smbinning(trainingData, y="ABOVE50K", x=continuous_var)  # WOE table
  if(class(smb) != "character"){  # any error while calculating scores.
    iv_df[iv_df$VARS == continuous_var, "IV"] <- smb$iv
  }
}

NAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercion

iv_df <- iv_df[order(-iv_df$IV), ]  # sort
iv_df

For more information: https://www.r-bloggers.com/woe-and-iv-variable-screening-with-information-in-r/

5. Build logit models and predict on test data

logitMod <- glm(ABOVE50K ~ RELATIONSHIP + AGE, data=trainingData, family=binomial(link="logit"))
predicted <- plogis(predict(logitMod, testData))  # predicted scores
# or
predicted <- predict(logitMod, testData, type="response")  # predicted scores

Why we used only RELATIONSHIP + AGE + CAPITALGAIN + OCCUPATION + EDUCATIONNUM.

Because of decide on optimal prediction probability cutoff for the model.

The default cutoff prediction probability score is 0.5 or the ratio of 1’s and 0’s in the training data. But sometimes, tuning the probability cutoff can improve the accuracy in both the development and validation samples. The InformationValue::optimalCutoff function provides ways to find the optimal cutoff to improve the prediction of 1’s, 0’s, both 1’s and 0’s and o reduce the misclassification error. Lets compute the optimal score that minimizes the misclassification error for the above model.

#install.packages("InformationValue")
#library(InformationValue)
optCutOff <- optimalCutoff(testData$ABOVE50K, predicted)[1] 
optCutOff

[1] 0.8771106

6. Do model diagnostics

The summary(logitMod) gives the beta coefficients, Standard error, z Value and p Value. If your model had categorical variables with multiple levels, you will find a row-entry for each category of that variable. That is because, each individual category is considered as an independent binary variable by the glm(). In this case it is ok if few of the categories in a multi-category variable don’t turn out to be significant in the model (i.e. p Value turns out greater than significance level of 0.5).

summary(logitMod)

Call:
glm(formula = ABOVE50K ~ RELATIONSHIP + AGE, family = binomial(link = "logit"), 
    data = trainingData)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0477  -0.7334   0.1133   0.8200   2.5855  

Coefficients:
                              Estimate Std. Error z value            Pr(>|z|)    
(Intercept)                  0.0003401  0.0899831   0.004               0.997    
RELATIONSHIP Not-in-family  -1.8881838  0.0569038 -33.182 <0.0000000000000002 ***
RELATIONSHIP Other-relative -3.1274740  0.2370408 -13.194 <0.0000000000000002 ***
RELATIONSHIP Own-child      -3.7651374  0.1567731 -24.016 <0.0000000000000002 ***
RELATIONSHIP Unmarried      -2.4644369  0.0936906 -26.304 <0.0000000000000002 ***
RELATIONSHIP Wife            0.1215044  0.0915656   1.327               0.185    
AGE                          0.0218334  0.0019637  11.119 <0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 15216  on 10975  degrees of freedom
Residual deviance: 11417  on 10969  degrees of freedom
AIC: 11431

Number of Fisher Scoring iterations: 5

VIF:-XXX

Like in case of linear regression, we should check for multicollinearity in the model. As seen below, all X variables in the model have VIF well below 4.

For more informarion: https://www.r-bloggers.com/collinearity-and-stepwise-vif-selection/

#install.packages("VIF")
#library(VIF)
vif(logitMod)

argument is not numeric or logical: returning NAError in as.vector(y) - mean(y) : non-numeric argument to binary operator

Misclassification Error:

Misclassification error is the percentage mismatch of predcited vs actuals, irrespective of 1’s or 0’s. The lower the misclassification error, the better is your model.

misClassError(testData$ABOVE50K, predicted, threshold = optCutOff)

[1] 0.0892

ROC (Receiver Operating Characteristics):

Receiver Operating Characteristics Curve traces the percentage of true positives accurately predicted by a given logit model as the prediction probability cutoff is lowered from 1 to 0. For a good model, as the cutoff is lowered, it should mark more of actual 1’s as positives and lesser of actual 0’s as 1’s. So for a good model, the curve should rise steeply, indicating that the TPR (Y-Axis) increases faster than the FPR (X-Axis) as the cutoff score decreases. Greater the area under the ROC curve, better the predictive ability of the model.

plotROC(testData$ABOVE50K, predicted)

The above model has area under ROC curve 88.78%, which is pretty good.

Concordance:

In simpler words, of all combinations of 1-0 pairs (actuals), Concordance is the percentage of pairs, whose scores of actual positive’s are greater than the scores of actual negative’s. For a perfect model, this will be 100%. So, the higher the concordance, the better is the quality of model.

Concordance(testData$ABOVE50K, predicted)

$`Concordance`
[1] 0.7969542

$Discordance
[1] 0.2030458

$Tied
[1] 0.00000000000000002775558

$Pairs
[1] 45252896

The above model with a concordance of 79.6% is indeed a good quality model.

Specificity and Sensitivity:-XXX

Sensitivity (or True Positive Rate) is the percentage of 1’s (actuals) correctly predicted by the model, while, specificity is the percentage of 0’s (actuals) correctly predicted. Specificity can also be calculated as 1???False Positive Rate.

#http://r-statistics.co/Logistic-Regression-With-R.html

sensitivity(testData$ABOVE50K, predicted, threshold = optCutOff)

Error in sensitivity.default(testData$ABOVE50K, predicted, threshold = optCutOff) : 
  inputs must be factors

Confusion Matrix:-XXX

confusionMatrix(testData$ABOVE50K, predicted, threshold = optCutOff)

Error: `data` and `reference` should be factors with the same levels.

Change log update

• 02.02.2019

References

• http://r-statistics.co/Logistic-Regression-With-R.html
• https://www.datacamp.com/community/tutorials/logistic-regression-R
• https://www.hackerearth.com/practice/machine-learning/machine-learning-algorithms/logistic-regression-analysis-r/tutorial/
• R Programming/Graphics
• Tutorials for learning R
• RDocumentation
• Statistical tools for high-throughput data analysis
• statmethods
• cognitiveclass.ai-Data Visualization with R
• Comprehensive Guide to Data Visualization in R
• datascienceplus
• Hands-On Programming with R
• R for Data Science
• R Markdown: The Definitive Guide
• R bookdown: Authoring Books and Technical Documents with R Markdown
• blogdown: Creating Websites with R Markdown
• r-statistics

License

MIT