R Study 5 - GLM

survey 에서 성별을 추측해보는 모델로 작업을 해봅니다.

library(MASS)
survey_tmp = survey
survey_tmp$Sex_Female = (survey_tmp$Sex == "Female")*1
input_ones = survey_tmp[which(survey_tmp$Sex == "Female"),]
input_zeros = survey_tmp[which(survey_tmp$Sex == "Male"),]

저렇게 하면 ones 에 여성이, zeros 에 남성이 들어갑니다

set.seed(100)
input_ones_training_rows = sample(1:nrow(input_ones), 0.7*nrow(input_ones))
input_zeros_training_rows = sample(1:nrow(input_zeros), 0.7*nrow(input_ones))
training_ones = input_ones[input_ones_training_rows, ]
training_zeros = input_zeros[input_zeros_training_rows, ]
trainingData = rbind(training_ones, training_zeros)
test_ones = input_ones[-input_ones_training_rows, ]
test_zeros = input_zeros[-input_zeros_training_rows, ]
testData = rbind(test_ones, test_zeros)

신기하네

이렇게 만들어서 테스트 데이터 는 test_ones 와 test_zeros 를 섞어서 넣고 training_ones 와 zeros 를 만들었다

library(smbinning)
factor_vars = c("W.Hnd", "Fold", "Clap", "Exer", "Smoke", "M.I")
continuous_vars = c("Wr.Hnd", "NW.Hnd", "Pulse", "Height", "Age")
iv_df = data.frame(VARS=c(factor_vars, continuous_vars), IV=numeric(11))
for(factor_var in factor_vars) {
  smb = smbinning.factor(trainingData, y="Sex_Female", x=factor_var)
  if ( class(smb) != "character") {
    iv_df[iv_df$VARS == factor_var, "IV"] = smb$iv
  }
}

강제형변환에 의해 생성된 NA 입니다강제형변환에 의해 생성된 NA 입니다

for(continuous_var in continuous_vars ) {
  smb = smbinning(trainingData, y="Sex_Female", x=continuous_var) 
  if ( class(smb) != "character") { 
    iv_df[iv_df$VARS == continuous_var, "IV"] = smb$iv
  }
}
iv_df = iv_df[order(-iv_df$IV), ]
iv_df

몇 가지 정리를 하면..

1. smbinning 을 돌렸을 때 리턴되는 smb 는, 딱히 매치가 안 된다면 여기서 “No significant split” 이라는 메시지가 나오고, 이것은 class() 를 넣으면 “character” 가 나온다. 말하자면 return -1 같은 느낌으로..

2. factor 랑 continuous 값들을 나눠서 각각을 따로 분류해서 IV값을 측정하는데, 그건 smbinning.factor() 나 smbinning() 으로 작업하면 각각을 만들어낼 수 있음. 만들어진 iv 값을 기준으로 어떤 값들이 더 가치있는 평가기준인지 확인할 수 있음.

3. 원래대로라면 얘들이 서로서로 관련이 없다는 게 증명되는 부분이 추가되어야 하지 않나?

table(trainingData$Sex_Female)

 0  1 
82 82 

균형은 잘 맞고

plot(density(survey_tmp$Sex_Female, na.rm=TRUE))

그래프를 보면 binomial 인게 확실하니까, 이건 glm 에서 binomial 로 + link 는 logit 으로 파악하면 된다. IV 테이블 상위 5개는 아래와 같으니 얘들을 그대로 쓰기로 하고

10 Height 2.3670
5 Smoke 0.0858
2 Fold 0.0244
3 Clap 0.0236
4 Exer 0.0094

model = glm( Sex_Female ~ Height + Smoke + Fold + Clap + Exer, data=trainingData, family=binomial(link="logit"))
testData$predicted = predict(model, testData, type="response")
summary(model)

Call:
glm(formula = Sex_Female ~ Height + Smoke + Fold + Clap + Exer, 
    family = binomial(link = "logit"), data = trainingData)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.19285  -0.36062   0.03434   0.52549   1.93776  

Coefficients:
            Estimate Std. Error z value     Pr(>|z|)    
(Intercept) 57.28964   10.42630   5.495 0.0000000391 ***
Height      -0.32045    0.05791  -5.534 0.0000000313 ***
SmokeNever  -0.75935    1.08505  -0.700       0.4840    
SmokeOccas  -1.76592    1.36317  -1.295       0.1952    
SmokeRegul  -3.78333    3.48797  -1.085       0.2781    
FoldNeither -2.28074    1.11378  -2.048       0.0406 *  
FoldR on L   0.22817    0.54579   0.418       0.6759    
ClapNeither -0.93497    0.86461  -1.081       0.2795    
ClapRight   -1.02665    0.79276  -1.295       0.1953    
ExerNone    -2.24143    0.87394  -2.565       0.0103 *  
ExerSome    -0.65517    0.64503  -1.016       0.3098    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 201.006  on 144  degrees of freedom
Residual deviance:  99.102  on 134  degrees of freedom
  (19 observations deleted due to missingness)
AIC: 121.1

Number of Fisher Scoring iterations: 7

library(car)
sqrt(vif(model)) > 2 

        GVIF    Df GVIF^(1/(2*Df))
Height FALSE FALSE           FALSE
Smoke  FALSE FALSE           FALSE
Fold   FALSE FALSE           FALSE
Clap   FALSE FALSE           FALSE
Exer   FALSE FALSE           FALSE

아까 3번에서 관계있는 거 가려내야 한다 했는데, 보니 쌍방 딱히 선형관계는 없는 모양이다.

library(InformationValue)
not_na_testData = testData[ complete.cases(testData['predicted']) , ]
threshold = optimalCutoff(not_na_testData$Sex_Female, not_na_testData$predicted)
misClassError(not_na_testData$Sex_Female, not_na_testData$predicted, threshold = threshold)

[1] 0.129

plotROC(testData$Sex_Female, testData$predicted)

커브가 대충 저렇게 나온다 뭐 이런 결론이 나옵니다.

AUROC 는 ROC 커브 아랫쪽 면적의 계산 (그러니까, 82%정도 맞추네요)

그럼 저기서 결론은

Height, SmokeNever, SmokeOccas, SmokeRegul, FoldNeither ExerNone ExerSome NW.Hnd

이렇게가 의미가 있다는 얘긴 것 같고

model = glm( Sex_Female ~ Height + Smoke + Fold + Exer + NW.Hnd, data=trainingData, family=binomial(link="logit"))
testData$predicted = predict(model, testData, type="response")
library(InformationValue)
plotROC(testData$Sex_Female, testData$predicted)

오 진짜네 조금 올랐네

참고로 기존에는 0.8271 이었습니다 8271 > 841

역시 사람보다 기계가 해야 (…)