AS4-1

Problem 1.1 - Exploration and Logistic Regression

We will first get familiar with the data. Load the CSV file gerber.csv into R. What proportion of people in this dataset voted in this election ?

gerber <- read.csv("gerber.csv")
temp <- table(gerber$voting)
temp[2] / sum(temp)

        1 
0.3158996 

Problem 1.2 - Exploration and Logistic Regression

Which of the four “treatment groups” had the largest percentage of people who actually voted (voting = 1)?

mean(gerber$voting[gerber$civicduty == 1])

[1] 0.3145377

mean(gerber$voting[gerber$hawthorne == 1])

[1] 0.3223746

mean(gerber$voting[gerber$self == 1])

[1] 0.3451515

mean(gerber$voting[gerber$neighbors== 1])

[1] 0.3779482

print("Neighbors" , quote = FALSE)

[1] Neighbors

# 利用平均值計算投票的比例

Problem 1.3 - Exploration and Logistic Regression

Build a logistic regression model for voting using the four treatment group variables as the independent variables (civicduty, hawthorne, self, and neighbors). Use all the data to build the model (DO NOT split the data into a training set and testing set). Which of the following coefficients are significant in the logistic regression model? Select all that apply.

model1 <- glm(voting ~ civicduty + hawthorne + self + neighbors, data = gerber, family = "binomial")
summary(model1)

Call:
glm(formula = voting ~ civicduty + hawthorne + self + neighbors, 
    family = "binomial", data = gerber)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9744  -0.8691  -0.8389   1.4586   1.5590  

Coefficients:
             Estimate Std. Error  z value Pr(>|z|)    
(Intercept) -0.863358   0.005006 -172.459  < 2e-16 ***
civicduty    0.084368   0.012100    6.972 3.12e-12 ***
hawthorne    0.120477   0.012037   10.009  < 2e-16 ***
self         0.222937   0.011867   18.786  < 2e-16 ***
neighbors    0.365092   0.011679   31.260  < 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: 429238  on 344083  degrees of freedom
Residual deviance: 428090  on 344079  degrees of freedom
AIC: 428100

Number of Fisher Scoring iterations: 4

print("Civicduty , Hawthorne , Self , Neighbors ", quote = FALSE)

[1] Civicduty , Hawthorne , Self , Neighbors 

Problem 1.4 - Exploration and Logistic Regression

Using a threshold of 0.3, what is the accuracy of the logistic regression model? (When making predictions, you don’t need to use the newdata argument since we didn’t split our data.)

pred1 <- predict(model1 , type = "response")
temp <- table(gerber$voting , pred1 > 0.3)
(temp[2,2] + temp[1,1]) / sum(temp)

[1] 0.5419578

# 將 confusion matrix 存到 temp 中，再計算 accuracy 的數值

Problem 1.5 - Exploration and Logistic Regression

Using a threshold of 0.5, what is the accuracy of the logistic regression model?

temp <- table(gerber$voting , pred1 > 0.5)
temp[1] / sum(temp)

[1] 0.6841004

# 同上，更改門檻值即可

Problem 1.6 - Exploration and Logistic Regression

Compare your previous two answers to the percentage of people who did not vote (the baseline accuracy) and compute the AUC of the model. What is happening here?

library(ROCR)
ROCRpred <- prediction(pred1, gerber$voting)
as.numeric(performance(ROCRpred, "auc")@y.values)

[1] 0.5308461

print("Even though all of the variables are significant, this is a weak predictive model." , quote = FALSE)

[1] Even though all of the variables are significant, this is a weak predictive model.

# AUC 僅高出 0.5 一點點點

Problem 2.1 - Trees

We will now try out trees. Build a CART tree for voting using all data and the same four treatment variables we used before. Don’t set the option method=“class” - we are actually going to create a regression tree here. We are interested in building a tree to explore the fraction of people who vote, or the probability of voting. We’d like CART to split our groups if they have different probabilities of voting. If we used method=‘class’, CART would only split if one of the groups had a probability of voting above 50% and the other had a probability of voting less than 50% (since the predicted outcomes would be different). However, with regression trees, CART will split even if both groups have probability less than 50%.

Leave all the parameters at their default values. You can use the following command in R to build the tree:

CARTmodel = rpart(voting ~ civicduty + hawthorne + self + neighbors, data=gerber) Plot the tree. What happens, and if relevant, why?

library(rpart)
library(rpart.plot)
modelCART <- rpart(voting ~ civicduty + hawthorne + self + neighbors, data = gerber)
prp(modelCART)

Problem 2.2 - Trees

Now build the tree using the command:

CARTmodel2 = rpart(voting ~ civicduty + hawthorne + self + neighbors, data=gerber, cp=0.0) to force the complete tree to be built. Then plot the tree. What do you observe about the order of the splits?

modelCART2 <- rpart(voting ~ civicduty + hawthorne + self + neighbors, data = gerber, cp=0.0)
prp(modelCART2)

print("Neighbor is the first split, civic duty is the last." , quote = FALSE)

[1] Neighbor is the first split, civic duty is the last.

# CP 表示為表示精度的參數

Problem 2.3 - Trees

Using only the CART tree plot, determine what fraction (a number between 0 and 1) of “Civic Duty” people voted:

0.31

[1] 0.31

# 從上圖可以看出來

Problem 2.4 - Trees

Make a new tree that includes the “sex” variable, again with cp = 0.0. Notice that sex appears as a split that is of secondary importance to the treatment group.

In the control group, which gender is more likely to vote?

modelCART3 <- rpart(voting ~ civicduty + hawthorne + self + neighbors + sex, data = gerber, cp=0.0)
prp(modelCART3)

print("Men" , quote = FALSE)

[1] Men

In the “Civic Duty” group, which gender is more likely to vote?

print("Men" , quote = FALSE)

[1] Men

Problem 3.1 - Interaction Terms

We know trees can handle “nonlinear” relationships, e.g. “in the ‘Civic Duty’ group and female”, but as we will see in the next few questions, it is possible to do the same for logistic regression. First, let’s explore what trees can tell us some more.

Let’s just focus on the “Control” treatment group. Create a regression tree using just the “control” variable, then create another tree with the “control” and “sex” variables, both with cp=0.0.

In the “control” only tree, what is the absolute value of the difference in the predicted probability of voting between being in the control group versus being in a different group? You can use the absolute value function to get answer, i.e. abs(Control Prediction - Non-Control Prediction). Add the argument “digits = 6” to the prp command to get a more accurate estimate.

modelCARTcontrol <- rpart(voting ~ control, data = gerber, cp = 0.0)
modelCARTcontrol_sex <- rpart(voting ~ control + sex, data = gerber, cp = 0.0)
prp(modelCARTcontrol, digits = 6)

abs(0.296638 - 0.34)

[1] 0.043362

Problem 3.2 - Interaction Terms

Now, using the second tree (with control and sex), determine who is affected more by NOT being in the control group (being in any of the four treatment groups):

prp(modelCARTcontrol_sex, digits = 6)

0.334176 - 0.290456

[1] 0.04372

0.345818 - 0.302795

[1] 0.043023

print("They are affected about the same (change in probability within 0.001 of each other)." , quote = FALSE)

[1] They are affected about the same (change in probability within 0.001 of each other).

Problem 3.3 - Interaction Terms

Going back to logistic regression now, create a model using “sex” and “control”. Interpret the coefficient for “sex”:

modelControl_Sex <- glm(voting ~ control + sex, data = gerber, family = "binomial")
summary(modelControl_Sex)

Call:
glm(formula = voting ~ control + sex, family = "binomial", data = gerber)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9220  -0.9012  -0.8290   1.4564   1.5717  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.635538   0.006511 -97.616  < 2e-16 ***
control     -0.200142   0.007364 -27.179  < 2e-16 ***
sex         -0.055791   0.007343  -7.597 3.02e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 429238  on 344083  degrees of freedom
Residual deviance: 428443  on 344081  degrees of freedom
AIC: 428449

Number of Fisher Scoring iterations: 4

print("Coefficient is negative, reflecting that women are less likely to vote " , quote = FALSE)

[1] Coefficient is negative, reflecting that women are less likely to vote 

# 由於 Sex = 1 代表女性，且係數是負的，代表說女性的投票率是較低的

Problem 3.4 - Interaction Terms

The regression tree calculated the percentage voting exactly for every one of the four possibilities (Man, Not Control), (Man, Control), (Woman, Not Control), (Woman, Control). Logistic regression has attempted to do the same, although it wasn’t able to do as well because it can’t consider exactly the joint possibility of being a women and in the control group.

We can quantify this precisely. Create the following dataframe (this contains all of the possible values of sex and control), and evaluate your logistic regression using the predict function (where “LogModelSex” is the name of your logistic regression model that uses both control and sex):

Possibilities = data.frame(sex=c(0,0,1,1),control=c(0,1,0,1)) predict(LogModelSex, newdata=Possibilities, type=“response”)

The four values in the results correspond to the four possibilities in the order they are stated above ( (Man, Not Control), (Man, Control), (Woman, Not Control), (Woman, Control) ). What is the absolute difference between the tree and the logistic regression for the (Woman, Control) case? Give an answer with five numbers after the decimal point.

Possibilities <- data.frame(sex = c(0,0,1,1) , control = c(0,1,0,1))
temp <- predict(modelControl_Sex, newdata = Possibilities, type = "response")
round(abs(0.2908065 - 0.290456),5)

[1] 0.00035

# 創建一個有這四種狀況的 data.frame

Problem 3.5 - Interaction Terms

So the difference is not too big for this dataset, but it is there. We’re going to add a new term to our logistic regression now, that is the combination of the “sex” and “control” variables - so if this new variable is 1, that means the person is a woman AND in the control group. We can do that with the following command:

LogModel2 = glm(voting ~ sex + control + sex:control, data=gerber, family=“binomial”)

How do you interpret the coefficient for the new variable in isolation? That is, how does it relate to the dependent variable?

model2 <- glm(voting ~ sex + control + sex:control, data = gerber, family = "binomial")
summary(model2)

Call:
glm(formula = voting ~ sex + control + sex:control, family = "binomial", 
    data = gerber)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9213  -0.9019  -0.8284   1.4573   1.5724  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.637471   0.007603 -83.843  < 2e-16 ***
sex         -0.051888   0.010801  -4.804 1.55e-06 ***
control     -0.196553   0.010356 -18.980  < 2e-16 ***
sex:control -0.007259   0.014729  -0.493    0.622    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 429238  on 344083  degrees of freedom
Residual deviance: 428442  on 344080  degrees of freedom
AIC: 428450

Number of Fisher Scoring iterations: 4

print("If a person is a woman and in the control group, the chance that she voted goes down." , quote = FALSE)

[1] If a person is a woman and in the control group, the chance that she voted goes down.

# 因為交叉相乘項的係數是負的，代表說是女性且為 control group 的投票率會下降

Problem 3.6 - Interaction Terms

Run the same code as before to calculate the average for each group:

predict(LogModel2, newdata=Possibilities, type=“response”)

Now what is the difference between the logistic regression model and the CART model for the (Woman, Control) case? Again, give your answer with five numbers after the decimal point.

predict(model2, newdata = Possibilities, type = "response")

        1         2         3         4 
0.3458183 0.3027947 0.3341757 0.2904558 

round(abs(0.2904558 - 0.290456),5)

[1] 0

Problem 3.7 - Interaction Terms

This example has shown that trees can capture nonlinear relationships that logistic regression can not, but that we can get around this sometimes by using variables that are the combination of two variables. Should we always include all possible interaction terms of the independent variables when building a logistic regression model?

print("No" , quote = FALSE)

[1] No

# 加入過多的變數容易導致模型過度適配、失去一般性的問題