빅분기실기연습(5)
R프로그램_로지스틱 회귀분석
합격 또는 불합격, 사망 또는 생존 등 연속된 수가 아닌 범주화된 종속변수가 있을 때 사용하는 방법이 로지스틱 회귀분석이다. 발생 확률이 0~1 사이에 있을 때 아래와 같이 X 값에 따라 Y 값은 0과 1 사이의 값을 가지게 된다. 이 경우, 확률 합수는 아래와 같이 로지스틱 함수의 형태가 된다.
이런 함수를 알기 위해서는 오즈(odds) 개념을 먼저 알아야 함. 성공할 확률이 p라고 하면 성공 확률/실패 확률이 오즈이다. odds=p/(1-p), 여기에 ln를 취하면 아래와 같다
ln(odds) = ln(p/(1-p))
우변을 테일러 정리하면
ln(p/(1-p)) = bo+ b1(x) + b2(x^2)+b3(x^3) ...
x^2 항 이후부터는 그 값이 작으므로 무시
ln(odds)=ln(p/1-p)) = bo+b1(x)
양변을 exp 취함
p/(1-p) = exp(bo+b1(x))
p에 대해 정리하면
p=1/[1+exp(-bo-b1x)]
즉, 로지스틱 회귀분석을 한다고 하면 회귀 함수의 계수인 bo와 b1을 찾는 것이다. linear 회귀 분석에서 y절편과 기울기를 찾는 것처럼.
# 예측 변수 x
x <- seq(-5, 5, length.out = 201)
# 임시로 회귀계수를 지정
beta0 <- 0
beta1 <- 1
# 확률 px
px <- 1 / (1 + exp(-(beta0 + beta1 * x)))
# 플롯
plot(x, px, type = "l", col = "red", ylim = c(0, 1))
# 예측 변수 x = 0 참조선 추가
lines(c(0, 0), c(0, 1), lty = "dashed", col = "blue")
text(0, 0.8, "x = 0")
# 확률 px = 0.5 참조선 추가
lines(c(-5, 5), c(0.5, 0.5), lty = "dashed", col = "blue")
text(-3, 0.55, "px = 0.5")
그럼 본격적으로 로지스틱 회귀 분석을 해 보자.
data는 GRE, GPA 출신학교 고등학교 Rank(1~4등급)에 따라 대학 합격 여부를 판단하는 것이다.
mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
> 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
admit의 0은 불합격, 1은 합격
data의 형태를 보자. admit과 rank는 연속형 변수가 아니므로 요인 값으로 바꾸어야 한다.
mydata$admit<-factor(mydata$admit)
mydata$rank<-factor(mydata$rank)
str(mydata)
> str(mydata)
'data.frame': 400 obs. of 4 variables:
$ admit: Factor w/ 2 levels "0","1": 1 2 2 2 1 2 2 1 2 1 ...
$ 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 ...
admit과 rank가 factor로 잘 바뀌었다.
그런 다음에 이제 로지스틱 회귀분석 glm을 사용하자. 합격/불합격이므로 family=binomial 을 사용하자
> model <- glm(admit ~ gre + gpa + rank, data =mydata, family = "binomial")
> summary(model)
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
출력 값을 보면 위에서 말한 bo에 해당하는 값이 -3.989979 이고, b1 에 해당하는 것이 gre~rank4 까지 각각 0.002264~-1.551464 로 나온다.
그리고 p value가 모두 0.05보다 작으므로 모든 변수가 유의한 변수로 사용할 수 있다.
그럼 각각 인자 1씩 변할 때 오즈비(odds = 성공확률/실패확률)을 구해보자. 그것은 각 인자에 exp를 취하면 된다.
exp(coef(model))
(Intercept) gre gpa rank2 rank3 rank4
0.0185001 1.0022670 2.2345448 0.5089310 0.2617923 0.2119375
즉 gre 1점이 증가하면 성공확률/실패 확률이 대략 1.002배이 증가한다는 것이다.
rank의 의미는 rank1 (좋은 학교)에서 rank2 로 바뀌면 오즈비가 0.5배로 떨어진다는 것이다.
만약 gre와 gpa는 평균값을 가지는 학생이 rank1~rank4 각각 학교에 다닐 때 붙을 수 있는 확률(admit)을 예측해보면 어떨까? newdata1을 아래와 같이 만들어보자
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
그런 다음 model을 이용하여 predict 함수를 써서 admit 확률 rankP를 구해보자
newdata1$rankP <- predict(model, newdata = newdata1, type = "response")
> 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
rank1은 붙을 확률이 50% 정도 되지만 동일 점수를 가지고 rank4에 다니면 붙는 확률이 18% 정도밖에 되지 않는다. 뭐지? 대학 입학에 고등학교별 점수가 매겨져 있는 것인가?
한 단계 더 나아가서 gre 점수가 200~800 점 사이에서 gpa는 평균값을 가지고 있고, 고등학교는 rank1 ~rank4 각각 다니는 학생들이 합격할 수 있는 범위는 얼마나 되는지 그래프로 표시해보자.
먼저 전체 gre 점수 200~800 사이에 100 등분을 해보자
seq(from = 200, to = 800, length.out = 100)
그 후에 rank 1~rank4에 각각 이 점수를 할당하려면 총 4번을 반복해야 함
rep(seq(from = 200, to = 800, length.out = 100),4)
rank는 1이 100번, 2가 100번, 3이 100번, 4가 100번 나오도록 함.
rank = factor(rep(1:4, each = 100)
gpa는 평균값을 가짐
gpa = mean(gpa)
이 모든 것을 newdata2에 넣는다.
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는 아래와 같이 400개 행으로 만들어진다.
> head(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
다음은 가장 중요한 것이다. 위에 만든 로지스틱 회귀분석 model을 이용하여 이 400개의 data를 넣어서 admit을 예측하는 것이다.
여기서 predict 함수 형식을 다시 한번 보자
predict(model, newdata, type = c("link", "response", "terms"))
여기서 type은 3종이 된다.
type = "link" 일 경우 log-odds 값이 출력되며, type = "response"의 경우 확률 $p$ 값이 출력됨.
따라서 우리는 log-odds를 출력하므로 link를 사용한다.
newdata3 <- cbind(newdata2, predict(model, newdata = newdata2, type = "link", se = TRUE))
다음으로 상한선(UL), 하한선(LL)을 넣어보자
newdata3 <- within(newdata3, {
PredictedProb <- plogis(fit)
LL <- plogis(fit - (1.96 * se.fit))
UL <- plogis(fit + (1.96 * se.fit))
})
plogis 함수를 link scale의 예측값에서 확률 추정 값 p를 구해준다. 즉 log(p/1-p) = Xb의 방정식을 풀어서 p를 구한다. 예를 들어 첫 번째 라인의 경우, log(p/(1-p)) = -0.811의 방정식을 양변에 exponential을 취해 풀면, 대략 p=0.308이 나온다.
> 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
마지막으로 ggplot2를 이용하여 그래프로 그려보자.
install.packages("ggplot2")
library(ggplot2)
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://3months.tistory.com/235 [Deep Play]rank1 ~rank4 학교에 다니는 학생들이 gpa는 평균점수를 가지면서 gre 점수가 200~800까지 분포할 때, 합격할 수 있는 확률을 잘 표현해주었다.
