Census Income Project

Loading Libraries

library(tidyverse)
library(caret)
library(randomForest)
library(skimr)
library(ggplot2)
library(gridExtra )
library(caTools)
library(corrplot)
library(ggcorrplot)
library(kableExtra)
library(caret)
library(tree)
library(stringr)

Reading the Data

cen<-read.csv("./census.csv")

1. Data Preprocessing

Exploring the data

str(cen)

'data.frame':   32561 obs. of  15 variables:
 $ age           : int  39 50 38 53 28 37 49 52 31 42 ...
 $ workclass     : chr  " State-gov" " Self-emp-not-inc" " Private" " Private" ...
 $ fnlwgt        : int  77516 83311 215646 234721 338409 284582 160187 209642 45781 159449 ...
 $ education     : chr  " Bachelors" " Bachelors" " HS-grad" " 11th" ...
 $ education.num : int  13 13 9 7 13 14 5 9 14 13 ...
 $ marital.status: chr  " Never-married" " Married-civ-spouse" " Divorced" " Married-civ-spouse" ...
 $ occupation    : chr  " Adm-clerical" " Exec-managerial" " Handlers-cleaners" " Handlers-cleaners" ...
 $ relationship  : chr  " Not-in-family" " Husband" " Not-in-family" " Husband" ...
 $ race          : chr  " White" " White" " White" " Black" ...
 $ sex           : chr  " Male" " Male" " Male" " Male" ...
 $ capital.gain  : int  2174 0 0 0 0 0 0 0 14084 5178 ...
 $ capital.loss  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ hours.per.week: int  40 13 40 40 40 40 16 45 50 40 ...
 $ native.country: chr  " United-States" " United-States" " United-States" " United-States" ...
 $ X             : chr  " <=50K" " <=50K" " <=50K" " <=50K" ...

summary(cen)

      age         workclass             fnlwgt         education         education.num  
 Min.   :17.00   Length:32561       Min.   :  12285   Length:32561       Min.   : 1.00  
 1st Qu.:28.00   Class :character   1st Qu.: 117827   Class :character   1st Qu.: 9.00  
 Median :37.00   Mode  :character   Median : 178356   Mode  :character   Median :10.00  
 Mean   :38.58                      Mean   : 189778                      Mean   :10.08  
 3rd Qu.:48.00                      3rd Qu.: 237051                      3rd Qu.:12.00  
 Max.   :90.00                      Max.   :1484705                      Max.   :16.00  
 marital.status      occupation        relationship           race               sex           
 Length:32561       Length:32561       Length:32561       Length:32561       Length:32561      
 Class :character   Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                                               
                                                                                               
                                                                                               
  capital.gain    capital.loss    hours.per.week  native.country          X            
 Min.   :    0   Min.   :   0.0   Min.   : 1.00   Length:32561       Length:32561      
 1st Qu.:    0   1st Qu.:   0.0   1st Qu.:40.00   Class :character   Class :character  
 Median :    0   Median :   0.0   Median :40.00   Mode  :character   Mode  :character  
 Mean   : 1078   Mean   :  87.3   Mean   :40.44                                        
 3rd Qu.:    0   3rd Qu.:   0.0   3rd Qu.:45.00                                        
 Max.   :99999   Max.   :4356.0   Max.   :99.00                                        

Counting NAs

colSums(is.na(cen))

           age      workclass         fnlwgt      education  education.num marital.status 
             0              0              0              0              0              0 
    occupation   relationship           race            sex   capital.gain   capital.loss 
             0              0              0              0              0              0 
hours.per.week native.country              X 
             0              0              0 

Insights

- There is a “?” in workclass, occupation and native country

- No NA values at this point

Trimming all the white spaces

cen %>% 
  mutate(across(where(is.character), str_trim))->cen

Replacing blank/missing values with NA

cen[cen==" "]<-NA

Let us count NAs once again

colSums(is.na(cen))

           age      workclass         fnlwgt      education  education.num marital.status 
             0              0              0              0              0              0 
    occupation   relationship           race            sex   capital.gain   capital.loss 
             0              0              0              0              0              0 
hours.per.week native.country              X 
             0              0              0 

Since we have few “?” in the data. Let us replace that with NA

cen[cen=='?']<-NA

Let us count NAs again

colSums(is.na(cen))

           age      workclass         fnlwgt      education  education.num marital.status 
             0           1836              0              0              0              0 
    occupation   relationship           race            sex   capital.gain   capital.loss 
          1843              0              0              0              0              0 
hours.per.week native.country              X 
             0            583              0 

Now we have NAs in occupation, workclass and native.country

Let us drop the rows with NAs

cen<-drop_na(cen)

Let us count NAs again

colSums(is.na(cen))

           age      workclass         fnlwgt      education  education.num marital.status 
             0              0              0              0              0              0 
    occupation   relationship           race            sex   capital.gain   capital.loss 
             0              0              0              0              0              0 
hours.per.week native.country              X 
             0              0              0 

All rows with NAs are now removed

2. Data Manipulation

a. Extracting education column

census_ed<-cen %>% 
  select(education)

b. Extracting columns from age to relationship

census_seq<-cen %>% 
  select(1:8)

c. Extracting columns 5,8, and 11

census_col<-cen %>% 
  select(5,8,11)

d. Extracting male employees who work in state_gov.

male_gov<- cen %>% 
  filter( sex=="Male", workclass=="State-gov")

e. Extracting all 39 year olds who either have a bachelor degree or is a US citizen

census_us<- cen %>% 
  filter( age == 39 ) %>% 
  filter(education=="Bachelors" | native.country=="United-States")

f. Extracting 200 random rows

set.seed(123)
census_200<-cen %>% 
  sample_n(200)

g. Count of different levels in workclass column

table(cen$workclass)

     Federal-gov        Local-gov          Private     Self-emp-inc Self-emp-not-inc 
             943             2067            22286             1074             2499 
       State-gov      Without-pay 
            1279               14 

h. Calculating the mean of capital.gain column grouped by workclass

cen %>% 
  group_by(workclass) %>% 
  summarise(Mean=mean(capital.gain)) %>% 
  kable()

workclass	Mean
Federal-gov	832.3213
Local-gov	829.2303
Private	879.8582
Self-emp-inc	4810.7467
Self-emp-not-inc	1913.1345
State-gov	684.3065
Without-pay	487.8571

i. Create a separate dataframe with the details of males and females from the census data that has income more than 50,000.

census_male_50k<-cen %>% 
  filter(sex == "Male", X==">50K")
census_female_50K<-cen %>% 
  filter(sex == "Female", X==">50K")

j. Calculate the percentage of people from the United States who are private employees and earn less than 50,000 annually.

Number<-cen %>% 
  filter(native.country == "United-States", workclass == "Private", X =="<=50K") %>% 
  summarise(Count = n())

Number_us<-cen %>% 
  filter(native.country== "United-States") %>% 
  summarise(Count_us=n())

Result<-(Number/Number_us)*100

as.character(Result)

[1] "56.6972076788831"

k. Calculate the percentage of married people in the census data

total<-cen %>% 
  nrow()

married<-cen %>% 
  filter(marital.status=="Married-AF-spouse"|marital.status==
           "Married-civ-spouse"|marital.status =="Married-spouse-absent") %>% 
  nrow()
  

Result1<-(married/total)*100
Result1

[1] 47.92786

l. Calculate the percentage of high school graduates earning more than

50,000 annually.

hs<-cen %>% 
  filter(education=="HS-grad") %>% 
  nrow()


grtr<-cen %>% 
  filter(education=="HS-grad" & X==">50K") %>% 
  nrow()

Result2<-(grtr/hs)*100
Result2

[1] 16.43293

3. Data Visualization

Relationship vs Race vs Sex

bar1<- cen %>% 
  ggplot(aes(x=relationship, fill=race))+
  geom_bar(position = "dodge") + xlab("Categories of Relationship")+ylab("Count of Categories")+
  ggtitle("Distribution of Relationships by Race")+
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

bar2<- cen %>% 
  ggplot(aes(x=relationship, fill=sex))+
  geom_bar(position = "dodge") + xlab("Categories of Relationship")+ylab("Count of Categories")+
  ggtitle("Distribution of Relationships by Sex")+
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))

grid.arrange(bar1,bar2, ncol=1)

Distribution of age

cen %>% 
  ggplot(aes(x=age, fill=X))+geom_histogram(bins=50)+
  guides(fill=guide_legend("Yearly Income"))+
  ggtitle("Distribution of Age")+
  theme_bw()

Scatter Plot

cen %>% 
  ggplot(aes(x=capital.gain, y=hours.per.week, color=X))+
  geom_point(size= 2, alpha=.40)+
  ggtitle("Capital Gain vs Hours per week by income")+
  xlab("Capital Gain")+ylab("Hours per week")+
  guides(color=guide_legend("Yearly Income"))

Box plot

cen %>% 
  ggplot(aes(x=education, y=age, fill=sex))+ geom_boxplot()+
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  ggtitle("Box plot of Age by Education and Sex")

4. Logistic Regression

a) Build a simple logistic regression model as follows:

Divide the dataset into training and test sets in 65:35 ratio.

sample.split(cen$X,SplitRatio = .65)->split_tag1
train1<-subset(cen,split_tag1==T)
test1<-subset(cen,split_tag1==F)

Build a logistic regression model where the dependent variable is “X”(yearly income) and the independent variable is “occupation”.

train1$X<-as.factor(train1$X)
train1$occupation<-as.factor(train1$occupation)
test1$X<-as.factor(test1$X)
test1$occupation<-as.factor(test1$occupation)

slm<-glm(X~occupation,data = train1, family = "binomial")
summary(slm)

Call:
glm(formula = X ~ occupation, family = "binomial", data = train1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1661  -0.7923  -0.5266  -0.1404   3.0414  

Coefficients:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                  -1.90552    0.06048 -31.504  < 2e-16 ***
occupationArmed-Forces      -10.66054  122.74159  -0.087  0.93079    
occupationCraft-repair        0.66926    0.07658   8.739  < 2e-16 ***
occupationExec-managerial     1.87875    0.07203  26.083  < 2e-16 ***
occupationFarming-fishing    -0.06819    0.13654  -0.499  0.61749    
occupationHandlers-cleaners  -0.90199    0.15636  -5.769 8.00e-09 ***
occupationMachine-op-inspct  -0.02626    0.10336  -0.254  0.79941    
occupationOther-service      -1.21740    0.12513  -9.729  < 2e-16 ***
occupationPriv-house-serv    -2.70960    1.00676  -2.691  0.00712 ** 
occupationProf-specialty      1.65054    0.07212  22.885  < 2e-16 ***
occupationProtective-serv     1.20172    0.11973  10.037  < 2e-16 ***
occupationSales               0.90792    0.07641  11.883  < 2e-16 ***
occupationTech-support        1.11657    0.10753  10.384  < 2e-16 ***
occupationTransport-moving    0.57402    0.09823   5.844 5.11e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 22002  on 19604  degrees of freedom
Residual deviance: 19459  on 19591  degrees of freedom
AIC: 19487

Number of Fisher Scoring iterations: 11

Predict the values on the test set.

pred<-predict(slm, newdata = test1, type = 'response')

Build a confusion matrix and find the accuracy

#Let us consider a threshold of 30%
table(test1$X,pred>.3)

       
        FALSE TRUE
  <=50K  6076 1853
  >50K   1156 1472

acc1<-(6076+1472)/(6076+1472+1156+1853)
acc1

[1] 0.7149758

b)Build a multiple logistic regression model as follows:

Divide the dataset into training and test sets in 80:20 ratio.

sample.split(cen$X,SplitRatio = .8)->split_tag2
train2<-subset(cen,split_tag2==T)
test2<-subset(cen,split_tag2==F)

Build a logistic regression model where the dependent variable is “X”(yearly income) and independent variables are “age”, “workclass”, and“education”.

train2$X<-as.factor(train2$X)
train2$occupation<-as.factor(train2$occupation)
test2$X<-as.factor(test2$X)
test2$occupation<-as.factor(test2$occupation)

mlm<-glm(X~age+workclass+education,data = train2, family = "binomial")
summary(mlm)

Call:
glm(formula = X ~ age + workclass + education, family = "binomial", 
    data = train2)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.45453  -0.71218  -0.47802  -0.00173   2.89906  

Coefficients:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)                -4.298522   0.189794 -22.648  < 2e-16 ***
age                         0.044976   0.001361  33.042  < 2e-16 ***
workclassLocal-gov         -0.505908   0.101988  -4.960 7.03e-07 ***
workclassPrivate           -0.284586   0.084748  -3.358 0.000785 ***
workclassSelf-emp-inc       0.700984   0.113648   6.168 6.91e-10 ***
workclassSelf-emp-not-inc  -0.434137   0.099636  -4.357 1.32e-05 ***
workclassState-gov         -0.709368   0.114761  -6.181 6.36e-10 ***
workclassWithout-pay      -13.026821 169.196138  -0.077 0.938630    
education11th               0.139019   0.216794   0.641 0.521360    
education12th               0.285689   0.283850   1.006 0.314186    
education1st-4th           -0.787269   0.487168  -1.616 0.106092    
education5th-6th           -0.372229   0.338497  -1.100 0.271482    
education7th-8th           -0.426406   0.252444  -1.689 0.091199 .  
education9th               -0.326394   0.280973  -1.162 0.245375    
educationAssoc-acdm         1.794219   0.176564  10.162  < 2e-16 ***
educationAssoc-voc          1.766029   0.172192  10.256  < 2e-16 ***
educationBachelors          2.509785   0.159602  15.725  < 2e-16 ***
educationDoctorate          3.699814   0.210276  17.595  < 2e-16 ***
educationHS-grad            1.091805   0.158773   6.877 6.13e-12 ***
educationMasters            2.939926   0.166614  17.645  < 2e-16 ***
educationPreschool        -10.995295  88.674904  -0.124 0.901319    
educationProf-school        3.736719   0.194735  19.189  < 2e-16 ***
educationSome-college       1.468588   0.159919   9.183  < 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: 27079  on 24128  degrees of freedom
Residual deviance: 22458  on 24106  degrees of freedom
AIC: 22504

Number of Fisher Scoring iterations: 12

Predict the values on the test set.

pred1<-predict(mlm, newdata = test2, type = "response")

Build a confusion matrix and find the accuracy.

#Lets us consider a threshold of 30%
table(test2$X,pred1>.3)

       
        FALSE TRUE
  <=50K  3586  945
  >50K    616  886

acc2<-(3527+891)/(3527+891+1004+611)
acc2

[1] 0.7323057

5. Decision Tree

a) Build a decision tree model as follows:

Divide the dataset into training and test sets in 70:30 ratio.

sample.split(cen$X, SplitRatio = .7)->split_tag3
train3<-subset(cen,split_tag3==T)
test3<-subset(cen, split_tag3==F)

Build a decision tree model where the dependent variable is “X”(Yearly Income) and the rest of the variables as independent variables.

# We need to first convert the independent variables to factors

train3$workclass<-as.factor(train3$workclass)
train3$sex<-as.factor(train3$sex)
train3$occupation<-as.factor(train3$occupation)
train3$education<-as.factor(train3$education)
train3$X<-as.factor(train3$X)
train3$marital.status<-as.factor(train3$marital.status)

# Building model

treemod<-tree(X~sex+occupation+education+marital.status+workclass, data = train3)
plot(treemod)
text(treemod)

Predict the values on the test set.

# We need to first convert the independent variables to factors

test3$workclass<-as.factor(test3$workclass)
test3$sex<-as.factor(test3$sex)
test3$occupation<-as.factor(test3$occupation)
test3$education<-as.factor(test3$education)
test3$X<-as.factor(test3$X)
test3$marital.status<-as.factor(test3$marital.status)

# Predicting 
pred3<-predict(treemod, newdata = test3, type = "class")

Build a confusion matrix and calculate the accuracy

table(test3$X,pred3)

       pred3
        <=50K >50K
  <=50K  6470  326
  >50K   1343  909

acc3<-(6470+926)/(6470+926+1326+926)
acc3

[1] 0.7665837

6. Random Forest:

a) Build a random forest model as follows:

Divide the dataset into training and test sets in 80:20 ratio.

sample.split(cen$X, SplitRatio = .8)->split_tag4
train4<-subset(cen,split_tag4==T)
test4<-subset(cen, split_tag4==F)

Build a random forest model where the dependent variable is “X”(Yearly Income)and the rest of the variables as independent variables and number of trees as 300.

# We need to first convert the independent variables to factors

train4$workclass<-as.factor(train4$workclass)
train4$sex<-as.factor(train4$sex)
train4$occupation<-as.factor(train4$occupation)
train4$education<-as.factor(train4$education)
train4$X<-as.factor(train4$X)
train4$marital.status<-as.factor(train4$marital.status)

#Building Model
randomForest(X~workclass+sex+occupation+education+marital.status, 
             data=train4, mtry =3, ntree =500)->RFmod

# Finding the importance of independent variables
importance(RFmod)

               MeanDecreaseGini
workclass              297.0319
sex                    132.4960
occupation             730.2710
education             1098.5017
marital.status        1780.9738

#Plotting
varImpPlot(RFmod)

Predict values on the test set

# We need to first convert the independent variables to factors

test4$workclass<-as.factor(test4$workclass)
test4$sex<-as.factor(test4$sex)
test4$occupation<-as.factor(test4$occupation)
test4$education<-as.factor(test4$education)
test4$X<-as.factor(test4$X)
test4$marital.status<-as.factor(test4$marital.status)

#Predicting

levels(test4$workclass) <- levels(train4$workclass)
pred4<-predict(RFmod,newdata = test4, type = "class")

Build a confusion matrix and calculate the accuracy

table(test4$X,pred4)

       pred4
        <=50K >50K
  <=50K  4136  395
  >50K    679  823

acc4<-(4136+823)/(4136+823+679+395)
acc4

[1] 0.8219791