R - Multinomial Logistic Regression(다중 로지스틱 회귀분석)

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/