R - 로지스틱 회귀분석


데이터 탐색

GRE, GPA, RANK이 입학(admission)에 어떤 영향을 주는지 로지스틱 회귀분석을 통해 분석한다.

library(aod) library(ggplot2) mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv") # view the first few rows of the data head(mydata) # 데이터의 대략적인 분포 확인 summary(mydata) # 데이터 구조 확인 str(mydata) # 변수별 표준편차 확인 sapply(mydata, sd) # contingency table : xtabs(~ x + y, data) xtabs(~admit+rank, data=mydata)

> head(mydata)
  admit gre  gpa rank
1     0 380 3.61    3
2     1 660 3.67    3
3     1 800 4.00    1
4     1 640 3.19    4
5     0 520 2.93    4
6     1 760 3.00    2

> summary(mydata) admit gre gpa rank Min. :0.0000 Min. :220.0 Min. :2.260 1: 61 1st Qu.:0.0000 1st Qu.:520.0 1st Qu.:3.130 2:151 Median :0.0000 Median :580.0 Median :3.395 3:121 Mean :0.3175 Mean :587.7 Mean :3.390 4: 67 3rd Qu.:1.0000 3rd Qu.:660.0 3rd Qu.:3.670 Max. :1.0000 Max. :800.0 Max. :4.000

> str(mydata) 'data.frame': 400 obs. of 4 variables: $ admit: int 0 1 1 1 0 1 1 0 1 0 ... $ gre : int 380 660 800 640 520 760 560 400 540 700 ... $ gpa : num 3.61 3.67 4 3.19 2.93 3 2.98 3.08 3.39 3.92 ... $ rank : Factor w/ 4 levels "1","2","3","4": 3 3 1 4 4 2 1 2 3 2 ...

> sapply(mydata, sd) admit gre gpa rank 0.4660867 115.5165364 0.3805668 0.9444602

> xtabs(~admit+rank, data=mydata) rank admit 1 2 3 4 0 28 97 93 55 1 33 54 28 12



GLM을 통한 로지스틱 회귀모형 구축


rank 변수를 factor 타입으로 변경시킨다.

glm 을통해 로지스틱회귀모형을 fitting시킨다. family='binomial' 인자를 통해, glm으로 로지스틱 회귀모형을 쓸 수 있다. (link function이 logit function이라는 의미)

mydata$rank <- factor(mydata$rank)
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")

summary(mylogit)


> summary(mylogit)

Call:
glm(formula = admit ~ gre + gpa + rank, family = "binomial", 
    data = mydata)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6268  -0.8662  -0.6388   1.1490   2.0790  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.989979   1.139951  -3.500 0.000465 ***
gre          0.002264   0.001094   2.070 0.038465 *  
gpa          0.804038   0.331819   2.423 0.015388 *  
rank2       -0.675443   0.316490  -2.134 0.032829 *  
rank3       -1.340204   0.345306  -3.881 0.000104 ***
rank4       -1.551464   0.417832  -3.713 0.000205 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 499.98  on 399  degrees of freedom
Residual deviance: 458.52  on 394  degrees of freedom
AIC: 470.52

Number of Fisher Scoring iterations: 4


로지스틱 회귀모형 결과 해석


Call : 구축한 모형에 대해 다시 상기시켜준다.

Deviance Residuals : Deviance residual에 대한 정보를 알려주는데, model fitting이 잘 되었는지에 대한 measure이다. 이를 통해 모델이 잘 적합됐는지를 평가할 수 있다.

Coefficient : 회귀계수와 그것들의 표준편차, z-statistics(wals's z-statistics), p-value를 나타낸다. 위 결과에서는 모든 변수가 유의한 것을 알 수 있다. 로지스틱 회귀모형에서는 회귀계수가 변수가 한 단위 증가했을 때 log(odds)의 증가량으로 해석할 수 있다.


example


gre가 1증가할 때, admission의 log odds(non-admission에 대한)가 0.0022 증가한다.

gpu의 경우도 마찬가지로 해석한다.

더미변수인 rank의 경우 해석이 약간 다른데, 예를 들어 rank2의 회귀계수 -0.67은 rank1에서 rank2로 바뀌었을 때, log(odds)의 변화량이다. 즉, rank1에 비해 rank2가 admission에 안좋은 영향을 준다는 것을 알 수 있다.



회귀계수들의 신뢰구간 얻기


Log-likelihood를 통해 구하는 법

confint(mylogit)
> confint(mylogit)
Waiting for profiling to be done...
                    2.5 %       97.5 %
(Intercept) -6.2716202334 -1.792547080
gre          0.0001375921  0.004435874
gpa          0.1602959439  1.464142727
rank2       -1.3008888002 -0.056745722
rank3       -2.0276713127 -0.670372346
rank4       -2.4000265384 -0.753542605


주어진 회귀계수의 표준편차를 이용해 신뢰구간을 구하는 법

## CIs using standard errors
confint.default(mylogit)
> confint.default(mylogit)
                    2.5 %       97.5 %
(Intercept) -6.2242418514 -1.755716295
gre          0.0001202298  0.004408622
gpa          0.1536836760  1.454391423
rank2       -1.2957512650 -0.055134591
rank3       -2.0169920597 -0.663415773
rank4       -2.3703986294 -0.732528724

이 신뢰구간은 이렇게도 계산 가능하다.예를 들어, gre의 95% 신뢰구간은 아래와 같이 구한다.

> 0.002264-1.96*0.001094
[1] 0.00011976
> 0.002264+1.96*0.001094
[1] 0.00440824


Wald Test를 통해 범주형 변수의 overall effect 파악하기


rank 변수를 dummy 변수로 만들어, reference인 rank1과 비교해 어떤 effect가 있는지 분석하였다. 하지만 aod 패키지의 walt.test를 이용하면, rank의 overall effect를 파악할 수 있다. [b:회귀계수, Sigma:error term의 공분산행렬, Terms:rank변수가 있는 열] wald.test는 순서를 통해 해당 범주형 변수의 위치를 파악하기 때문에 순서를 잘 신경 써야 한다.

wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)
> wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)
Wald test:
----------

Chi-squared test:
X2 = 20.9, df = 3, P(> X2) = 0.00011

카이스퀘어 테스트의 p-value가 0.00011이므로 rank 변수는 유의하다.



변수간의 회귀계수가 같은지 검정하기


wald test를 통해 rank가 admission에 유의한 효과가 있다는 것을 파악하였고, rank1이라는 reference에 비해 rank2와 rank2가 유의한 영향을 준다는 것도 확인하였다. 한가지 의문은 rank2, rank3, rank4 의 회귀계수가 동일한지에 대한 것이다. 회귀 계수가 동일하다면, rank1일 때에 비해 rank2, rank3, rank4의 효과가 동일하다는 것을 의미한다.


이 구문은 rank2의 효과와 rank3의 효과가 동일한지를 검정한다.

l <- cbind(0, 0, 0, 1, -1, 0)
wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)
> l <- cbind(0, 0, 0, 1, -1, 0)
> wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)
Wald test:
----------

Chi-squared test:
X2 = 5.5, df = 1, P(> X2) = 0.019


검정 결과, rank2, rank3의 효과는 유의하다 다른 것으라 나타났다.



OR과 회귀계수의 관계


앞서 로지스틱 회귀분석에서의 회귀계수는 log odds의 증가량이라고 언급하였다. 그러면 회귀계수에 exponential을 취하면, OR의 증분이 된다.

exp(coef(mylogit))
> exp(coef(mylogit))
(Intercept)         gre         gpa       rank2       rank3       rank4 
  0.0185001   1.0022670   2.2345448   0.5089310   0.2617923   0.2119375 


해석 ex) gpa가 1증가하면 admission의 non-admission에 대한 OR이 2.23배 증가한다.



예측하기


1. 다른 변수들을 고정하고 rank가 변할 때 예측값의 변화를 보기

newdata1 <- with(mydata, data.frame(gre = mean(gre), gpa = mean(gpa), rank = factor(1:4)))

> newdata1
    gre    gpa rank
1 587.7 3.3899    1
2 587.7 3.3899    2
3 587.7 3.3899    3
4 587.7 3.3899    4

newdata1$rankP <- predict(mylogit, newdata = newdata1, type = "response") # probability 점추정값 newdata1

> newdata1
    gre    gpa rank     rankP
1 587.7 3.3899    1 0.5166016
2 587.7 3.3899    2 0.3522846
3 587.7 3.3899    3 0.2186120
4 587.7 3.3899    4 0.1846684

2. rank, gpa를 고정한 후, gre의 효과 보기

newdata2 <- with(mydata, data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100),
                                              4), gpa = mean(gpa), rank = factor(rep(1:4, each = 100))))
newdata2
> newdata2
         gre    gpa rank
1   200.0000 3.3899    1
2   206.0606 3.3899    1
3   212.1212 3.3899    1
4   218.1818 3.3899    1
5   224.2424 3.3899    1
6   230.3030 3.3899    1
7   236.3636 3.3899    1
8   242.4242 3.3899    1
9   248.4848 3.3899    1
10  254.5455 3.3899    1
newdata3 <- cbind(newdata2, predict(mylogit, newdata = newdata2, type = "link", 
                                    se = TRUE))
newdata3 <- within(newdata3, {
  PredictedProb <- plogis(fit)
  LL <- plogis(fit - (1.96 * se.fit))
  UL <- plogis(fit + (1.96 * se.fit))
})
> head(newdata3)
       gre    gpa rank        fit    se.fit residual.scale        UL        LL PredictedProb
1 200.0000 3.3899    1 -0.8114870 0.5147714              1 0.5492064 0.1393812     0.3075737
2 206.0606 3.3899    1 -0.7977632 0.5090986              1 0.5498513 0.1423880     0.3105042
3 212.1212 3.3899    1 -0.7840394 0.5034491              1 0.5505074 0.1454429     0.3134499
4 218.1818 3.3899    1 -0.7703156 0.4978239              1 0.5511750 0.1485460     0.3164108
5 224.2424 3.3899    1 -0.7565919 0.4922237              1 0.5518545 0.1516973     0.3193867
6 230.3030 3.3899    1 -0.7428681 0.4866494              1 0.5525464 0.1548966     0.3223773


이 때 type="link" 옵션을 주면, link scale로 예측값을 내준다. 이는 logit 변환을 하기전의 예측값이다. 이를 통해 예측한 확률의 신뢰구간을 구할 수 있다. 왜냐하면 이 link scale의 예측값 (fit 변수)의 표준편차를 쉽게 알 수 있기 때문이다. 따라서 link scale의 예측값의 신뢰구간을 구한후 이를 다시 logit 변환하여 예측 확률의 신뢰구간을 얻을 수 있다. 이를 구현한 것이 바로 위의 코드이다.


plogis 함수를 link scale의 예측값에서 확률 추정값 p를 구해준다. 즉 log(p/1-p) = Xb 의 방정식을 풀어서 p를 구한다. 예를 들어 첫번째 라인의 경우, log(p/(1-p)) = -0.811 의 방정식을 양변에 exponential을 취해 풀면, 대략 p=0.308이 나온다.



newdata3의 예측값 시각화하기


rank를 고정시켜놓고 gre가 admission 확률에 어떤 영향을 주는지를 시각화

ggplot(newdata3, aes(x = gre, y = PredictedProb)) + geom_ribbon(aes(ymin = LL,
                                                                    ymax = UL, fill = rank), alpha = 0.2) + geom_line(aes(colour = rank),
                                                                                                                      size = 1)





출처

https://stats.idre.ucla.edu/r/dae/logit-regression/

  • 2018.04.25 23:26

    비밀댓글입니다

  • 2018.08.08 09:51

    정리 잘 해주셨습니다. 감사해요.

  • 2018.10.31 10:05

    비밀댓글입니다

    • Deepplay 2018.10.31 23:37 신고

      안녕하세요. 말씀하신 경우처럼 각각의 지역이 reference에 비해 유의한지 보는 게 아니라, "지역" 이라는 변수가 유의한지를 보고 싶은 경우, anova라고 하는 F 검정을 이용하시면 될 것 같습니다. F 검정을 하면 지역에 의한 변동량이 총 변동량 대비 유의한지를 검정하게 됩니다. SAS의 경우 proc glm을 이용하면 각 변수에 대한 F 검정 표가 나오게 됩니다. R에서는 aov 함수를 사용하시면 될 것 같습니다.