R - Multinomial Logistic Regression
Multinomial Logistic Regression이란 y의 범주가 3개 이상(multi)이며 명목형(nomial)일 때 사용하는 로지스틱 회귀분석이다. 본 포스팅은 이 사이트의 튜토리얼을 보고 실습하고 해석한 자료이다.
데이터 설명
해당 데이터셋은 hsbdemo 데이터셋으로 총 n은 200이다.
ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")
> head(ml)
id female ses schtyp prog read write math science socst honors awards cid
1 45 female low public vocation 34 35 41 29 26 not enrolled 0 1
2 108 male middle public general 34 33 41 36 36 not enrolled 0 1
3 15 male high public vocation 39 39 44 26 42 not enrolled 0 1
4 67 male low public vocation 37 37 42 33 32 not enrolled 0 1
5 153 male middle public vocation 39 31 40 39 51 not enrolled 0 1
6 51 female high public general 42 36 42 31 39 not enrolled 0 1
Y : prog (프로그램 타입)
X : ses(social economic status), write (writing score)
다중회귀분석을 통해 X가 Y에 어떤 영향을 미치는지 알아본다.
빈도 테이블과 그룹 평균을 통해 X와 Y의 관계 알아보기
with(ml, table(ses, prog))
# 또는 table(ml$ses, ml$prog)
# 또는 xtabs(~ses+prog, data=ml)
> with(ml, table(ses, prog))
prog
ses general academic vocation
low 16 19 12
middle 20 44 31
high 9 42 7
with(ml, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))
> with(ml, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))
M SD
general 51.33333 9.397775
academic 56.25714 7.943343
vocation 46.76000 9.318754
nnet 패키지의 multinom 함수를 통해 다중 로지스틱 회귀분석 실시
ml$prog2 <- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)
multinom을 실행하기 전에 relevel이라는 함수를 통해 reference를 지정해준다. 이 reference를 기준으로 결과를 해석할 수 있다. 따라서 이 reference는 어떤 baseline이 되어야한다. 이 경우에 academic이 baseline이며 이를 기준으로 하여 분석한다.
multinom 함수를 실행하면 모델이 fitting되어 계수가 결정되고 이 정보는 test에 담긴다.
summary(test)
> summary(test)
Call:
multinom(formula = prog2 ~ ses + write, data = ml)
Coefficients:
(Intercept) sesmiddle seshigh write
general 2.852198 -0.5332810 -1.1628226 -0.0579287
vocation 5.218260 0.2913859 -0.9826649 -0.1136037
Std. Errors:
(Intercept) sesmiddle seshigh write
general 1.166441 0.4437323 0.5142196 0.02141097
vocation 1.163552 0.4763739 0.5955665 0.02221996
Residual Deviance: 359.9635
AIC: 375.9635
이 예제의 경우 다중 로지스틱 회귀분석의 계수는 2세트가 나오게 된다.
1. academic vs general 의 log odds에 관한 계수
2. academic vs vocation의 log odds에 관한 계수
결과 해석
- write가 1단위 증가할 때, academic이 general이 될 log odds가 -0.058이 된다.
- write가 1단위 증가할 때, academic이 vocation이 될 log odds가 -0.1136037이 된다.
- ses:row에서 ses:high로 변할 때, academic이 general이 될 log odds가 -1.1628이 된다.
- ses:row에서 ses:middle로 변할 때, academic이 general이 될 log odds가 -0.5332가 된다.
- ses:row에서 ses:high로 변할 때, academic이 vocation이 될 log odds가 -0.9826이 된다.
- ses:row에서 ses:middle로 변할 때, academic이 vocation이 될 log odds가 +0.29가 된다.
변수의 유의성 판단
coefficient가 0라는 귀무가설 하에서 coef ~ N(0, sd^2)라는걸 이용하여 z-test 한다.
z <- summary(test)$coefficients/summary(test)$standard.errors
z
> z <- summary(test)$coefficients/summary(test)$standard.errors
> z
(Intercept) sesmiddle seshigh write
general 2.445214 -1.2018081 -2.261334 -2.705562
vocation 4.484769 0.6116747 -1.649967 -5.112689# 2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p
> p <- (1 - pnorm(abs(z), 0, 1)) * 2
> p
(Intercept) sesmiddle seshigh write
general 0.0144766100 0.2294379 0.02373856 6.818902e-03
vocation 0.0000072993 0.5407530 0.09894976 3.176045e-07
OR 추정값
exp(coef(test))
> exp(coef(test))
(Intercept) sesmiddle seshigh write
general 17.32582 0.5866769 0.3126026 0.9437172
vocation 184.61262 1.3382809 0.3743123 0.8926116
회귀 계수들에 exponential을 취하면, 이는 OR에 대한 추정값이 된다.
Write가 평균일 때, ses에 따라 예측값이 어떻게 달라지는가?
dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
dses
predict(test, newdata = dses, "probs")> dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
> dses
ses write
1 low 52.775
2 middle 52.775
3 high 52.775
> predict(test, newdata = dses, "probs")
academic general vocation
1 0.4396845 0.3581917 0.2021238
2 0.4777488 0.2283353 0.2939159
3 0.7009007 0.1784939 0.1206054
Write를 변화시켜가면서 예측값의 변화 관찰하기
dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), write = rep(c(30:70),
3))
## store the predicted probabilities for each value of ses and write
pp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))
## calculate the mean probabilities within each level of ses
by(pp.write[, 3:5], pp.write$ses, colMeans)
## melt data set to long for ggplot2
lpp <- melt(pp.write, id.vars = c("ses", "write"), value.name = "probability")
head(lpp) # view first few rows> head(lpp) # view first few rows
ses write variable probability
1 low 30 academic 0.09843588
2 low 31 academic 0.10716868
3 low 32 academic 0.11650390
4 low 33 academic 0.12645834
5 low 34 academic 0.13704576
6 low 35 academic 0.14827643
## plot predicted probabilities across write values for each level of ses
## facetted by program type
ggplot(lpp, aes(x = write, y = probability, colour = ses)) + geom_line() + facet_grid(variable ~
., scales = "free")
- 이 그래프는 식을 전개하여 y를 sigmoid 함수 형태로 만든 것이다. 2개의 적합식에서 3개의 함수를 만들어낼 수 있다.
Y~write의 결과랑 무엇이 다른가?
예측의 관점에서 보면 ses를 빼고 write만을 통해 Y를 예측하면 예측을 세분화하지 못하여 예측력이 떨어진다고 볼 수 있다. 이는 ses를 고려하지 않고, write와 Y의 관계만을 보기 때문이다. 따라서 ses가 추가됨으로써 write의 계수가 바뀐다는 것을 알 수 있다. 이 과정을 ses를 보정한다고도 말한다.
출처 - https://stats.idre.ucla.edu/r/dae/multinomial-logistic-regression/
sample.csv