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/
'Tools > R' 카테고리의 다른 글
| R - 종속변수가 순서 척도인 로지스틱 회귀분석 (4) | 2017.12.03 |
|---|---|
| R - Multinomial Logistic Regression(다중 로지스틱 회귀분석) (0) | 2017.11.22 |
| R - Reshape를 통한 데이터 재구성 (melt, reshape) (0) | 2017.11.15 |
| R - 랜덤포레스트 실습 (0) | 2017.11.11 |
| 우분투에 R과 Rstudio Server 설치 (0) | 2017.11.05 |
