AS4-1 投票行為的影響因素 (temp.)

主要議題：研究不同的提示訊息對投票行為的影響

學習重點：

• 用實驗設計確認應果關係
• 用回歸樹模型預測兩類別的機率
• 虛擬變數(dummy varable)
• 決策樹模型和邏輯式回歸模型的比較
• 交互作用項的使用

rm(list=ls(all=T))
options(digits=4, scipen=12)
library(dplyr)
library(rpart)
library(rpart.plot)
library(caret)
library(randomForest)
library(caTools)
source('DPP.R')

---

1. 資料探索與邏輯式回歸

1.1 讀進資料

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?

g = read.csv("data/gerber.csv")
mean(g$voting)                   # 31.6%

[1] 0.3159

從344084的人當中，找出會投票的人，用平均數mean()來找

1.2 Examine Data

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

sapply(g[,4:8], function(x) mean(g$voting[x == 1])) %>% sort #neighbors

  control civicduty hawthorne      self neighbors 
   0.2966    0.3145    0.3224    0.3452    0.3779 

sapply()是將字串轉換成向量

tapply(g$voting, g$civicduty, mean)

     0      1 
0.3161 0.3145 

tapply(g$voting, g$hawthorne, mean)

     0      1 
0.3151 0.3224 

tapply(g$voting, g$self, mean)

     0      1 
0.3122 0.3452 

tapply(g$voting, g$neighbors, mean)

     0      1 
0.3082 0.3779 

#neighbors

1.3 Logistic Model

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.

glm1 = glm(voting ~   hawthorne + civicduty + self + neighbors, 
           data=g, family="binomial")
summary(glm1)   #hawthorne、civicduty、neighbors、self

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

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-0.974  -0.869  -0.839   1.459   1.559  

Coefficients:
            Estimate Std. Error z value        Pr(>|z|)    
(Intercept) -0.86336    0.00501 -172.46         < 2e-16 ***
hawthorne    0.12048    0.01204   10.01         < 2e-16 ***
civicduty    0.08437    0.01210    6.97 0.0000000000031 ***
self         0.22294    0.01187   18.79         < 2e-16 ***
neighbors    0.36509    0.01168   31.26         < 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

1.4 Model Accuracy (threshold = 0.3)

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

pred = predict(glm1, type="response")
table(g$voting, pred > 0.3) %>% {sum(diag(.))/sum(.)} # 0.5420

[1] 0.542

pred=predict(glm1, type="response")
table(g$voting, pred>0.3)

   
     FALSE   TRUE
  0 134513 100875
  1  56730  51966

(134513+51966)/(134513+100875+56730+51966) #0.5420

[1] 0.542

1.5 Model Accuracy (threshold = 0.5)

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

table(g$voting, pred > 0.5) %>% {sum(diag(.))/sum(.)}  # ACC.glm1 = 0.6841

[1] 0.6841

pred=predict(glm1, type="response")
table(g$voting, pred>0.5)

   
     FALSE
  0 235388
  1 108696

(235388+0)/(235388+108696)   #0.6841

[1] 0.6841

1.6 AUC

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?

colAUC(pred, g$voting)  # AUC.glm1 = 0.5308 

          [,1]
0 vs. 1 0.5308

library(ROCR)
ROCRpred=prediction(pred, g$voting)
as.numeric(performance(ROCRpred, "auc")@y.values) #0.5308

[1] 0.5308

---

2. 決策樹模型

2.1 Default Regression Tree

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:

m1 = rpart(voting ~ civicduty + hawthorne + self + neighbors, data = g)
rpart.plot(m1)

# No variables are used (the tree is only a root node) - none of the variables make a big enough effect to be split on

2.2, 2.3 Tree with Zero Cost

Now build the tree using the command:

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

m2 = rpart(voting ~ civicduty+ hawthorne+self+neighbors, data=g, cp=0.0)
prp(m2)

# Neighbor is the first split, civic duty is the last.

2.4 Tree with Zero Cost

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?

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

m3 = rpart(voting ~ civicduty+ hawthorne+self+neighbors+sex, data=g, cp=0.0)
prp(m3)

# Neighbor is the first split, civic duty is the last.

---

3. 加入交互作用項

3.1

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.

par(mfcol=c(2,1))
m4 = rpart(voting ~ control, data=g, cp=0.0)
prp(m4, digits=6)
m5 = rpart(voting ~ control + sex, data=g, cp=0.0)
prp(m5, digits=6)

0.34 - 0.296638

[1] 0.04336

當control=1時，預測值為0.296638；當control=0，預測值為0.34，兩者相減為0.04336

3.2 Predictor’s Conditional Effect

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):

options(digits=6)
m5

n= 344084 

node), split, n, deviance, yval
      * denotes terminal node

1) root 344084 74359.0 0.315900  
  2) control>=0.5 191243 39901.7 0.296638  
    4) sex>=0.5 95419 19665.0 0.290456 *
    5) sex< 0.5 95824 20229.4 0.302795 *
  3) control< 0.5 152841 34297.5 0.340000  
    6) sex>=0.5 76376 16993.8 0.334176 *
    7) sex< 0.5 76465 17298.5 0.345818 *

0.334176 - 0.290456   # female: 0.04372

[1] 0.04372

0.345818 - 0.302795   # male:   0.043023

[1] 0.043023

(0.334176 - 0.290456) - (0.345818 - 0.302795)  # diff: 0.000697

[1] 0.000697

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

3.3 Calculate Effect from Logistic Reg. Coef.

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

options(digits=6)
m6 = glm(voting ~ control + sex, g, family="binomial")
summary(m6)

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

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-0.922  -0.901  -0.829   1.456   1.572  

Coefficients:
            Estimate Std. Error z value         Pr(>|z|)    
(Intercept) -0.63554    0.00651   -97.6          < 2e-16 ***
control     -0.20014    0.00736   -27.2          < 2e-16 ***
sex         -0.05579    0.00734    -7.6 0.00000000000003 ***
---
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

# Coefficient is negative, reflecting that women are less likely to vote

3.4 Comparing Estimations of Effect

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):

new = data.frame(sex=c(0,0,1,1),control=c(0,1,0,1)); new

predict(m6, newdata=new, type="response")

       1        2        3        4 
0.346256 0.302446 0.333738 0.290806 

# (Woman, Control)---logistic regression預測數0.290806
# (Woman, Control)---CART預測數0.296638

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.

abs(0.290806 - 0.290456) # 0.00035

[1] 0.00035

3.5 glm() with Interaction Term

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:

m7 = glm(voting ~ sex * control, g, family="binomial")
summary(m7)

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

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-0.921  -0.902  -0.828   1.457   1.573  

Coefficients:
            Estimate Std. Error z value  Pr(>|z|)    
(Intercept) -0.63747    0.00760  -83.84   < 2e-16 ***
sex         -0.05189    0.01080   -4.80 0.0000016 ***
control     -0.19655    0.01036  -18.98   < 2e-16 ***
sex:control -0.00726    0.01473   -0.49      0.62    
---
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

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

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

3.6 Compare the Effect Again

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

predict(m7, newdata=new, type="response")

       1        2        3        4 
0.345818 0.302795 0.334176 0.290456 

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.

(0.290456 - 0.290456)

[1] 0

(Woman, Control)logistic regression預測值0.2904558

3.7 Should we always include 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?

# No

---