R - 종속변수가 순서 척도인 로지스틱 회귀분석

R - 종속변수가 순서 척도인 로지스틱 회귀분석

다항 로지스틱 회귀분석에서 종속변수가 명목형이 아닌 순서형인 경우 분석하는 방법에 관한 포스팅입니다. 이 포스팅은 이 사이트의 튜토리얼을 보고 작성하였습니다.

라이브러리 로드

require(foreign)
require(ggplot2)
require(MASS)
require(Hmisc)
require(reshape2)

데이터 읽기

dat <- read.dta("https://stats.idre.ucla.edu/stat/data/ologit.dta")
head(dat)

> head(dat)
            apply pared public  gpa
1     very likely     0      0 3.26
2 somewhat likely     1      0 3.21
3        unlikely     1      1 3.94
4 somewhat likely     0      0 2.81
5 somewhat likely     0      0 2.53
6        unlikely     0      1 2.59

X : pared(부모중 한명이 대학 학위가 있는지 여부), public(졸업한 학교가 공립인지), gpa(학점)

Y : apply (“unlikely”, “somewhat likely”, “very likely”, 각각 1, 2, 3으로 코딩되어있음) => 대학원에 진학할 가능성

이 때 Y는 순서형 척도이다.

범주형 변수의 빈도 확인

lapply 함수를 이용해 범주형 변수의 빈도표를 확인

## one at a time, table apply, pared, and public
lapply(dat[, c("apply", "pared", "public")], table)

> lapply(dat[, c("apply", "pared", "public")], table)
$apply

       unlikely somewhat likely     very likely 
            220             140              40 

$pared

  0   1 
337  63 

$public

  0   1 
343  57 

범주형 독립변수와 종속변수 관계 파악 : 3-변수 빈도표

변수 빈도표를 통해 로지스틱 회귀분석을 하기에 앞서 예상되는 결과를 미리 확인한다.

## three way cross tabs (xtabs) and flatten the table
ftable(xtabs(~ public + apply + pared, data = dat))

> ftable(xtabs(~ public + apply + pared, data = dat))
                       pared   0   1
public apply                        
0      unlikely              175  14
       somewhat likely        98  26
       very likely            20  10
1      unlikely               25   6
       somewhat likely        12   4
       very likely             7   3

연속형 변수 분포 확인

summary(dat$gpa)

> summary(dat$gpa)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.900   2.720   2.990   2.999   3.270   4.000 

sd(dat$gpa)

> sd(dat$gpa)
[1] 0.3979409

연속형 독립 변수와 종속 변수와의 관계 파악

Apply에 따라 gpa가 어떻게 분포하여있는지 상자그림을 통해 확인한다. 이러한 경향을 수치화하는 것이 로지스틱 회귀분석의 목표이다.

ggplot(dat, aes(x = apply, y = gpa)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) +
  facet_grid(pared ~ public, margins = TRUE) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

지금까지는 데이터에 대해서 살펴보았고, 본격적으로 순서 척도가 종속변수인 로지스틱 회귀분석을 어떻게 하는지 알아본다.

이 때 가능한 분석방법은 꼭 순서형 척도만은 아닐 것이다. GPA만을 독립변수로하여 ANOVA를 할 수도 있고, y를 1,2,3으로 놓고 다중회귀분석을 할 수도 있다. 또 y를 순서가 없는 명목형 척도로 보고 명목형 다항 로지스틱 분석(multinomial logistic regression)을 할 수도 있을 것이다. 하지만 해당 데이터에 대해 위 분석 방법들은 가정에 맞지 않는다. 따라서 순서형 척도를 종속변수로하는 로지스틱 회귀분석(ordinal logistic regression)이 적합하다.

순서형 척도를 종속변수로하는 로지스틱 회귀분석의 중요한 가정은 회귀계수가 같다고 가정하는 것이다. y=1,2,3인 이  데이터의 경우 어떠한 변수의 변화가 1 vs 2,3의 OR에 미치는 영향과 1,2 vs 3의 OR에 미치는 영향이 같다고 가정한다. y의 범주가 3개인 경우 2개의 식이 나오면 이 두식은 회귀계수가 같고, 단지 절편만 다르다.

MASS 패키지의 polr 함수를(proportional odds logistic regession) 이용하여 순서형 변수를 종속변수로하는 다항 로지스틱 회귀분석을 해보자.

회귀식을 적어주고 summary를 통해 적합된 결과를 본다. Hess=TRUE는 optimization을 할 때 사용하는 Hessian matrix를 의미한다. (https://en.wikipedia.org/wiki/Hessian_matrix)

## fit ordered logit model and store results 'm'
m <- polr(apply ~ pared + public + gpa, data = dat, Hess=TRUE)

## view a summary of the model
summary(m)

> summary(m)
Call:
polr(formula = apply ~ pared + public + gpa, data = dat, Hess = TRUE)

Coefficients:
          Value Std. Error t value
pared   1.04769     0.2658  3.9418
public -0.05879     0.2979 -0.1974
gpa     0.61594     0.2606  2.3632

Intercepts:
                            Value   Std. Error t value
unlikely|somewhat likely     2.2039  0.7795     2.8272
somewhat likely|very likely  4.2994  0.8043     5.3453

Residual Deviance: 717.0249 
AIC: 727.0249 

결과 해석

회귀계수에 exponential을 취하면 OR(Odds ratio)이 되며 이를 통해 결과를 해석해본다.

## odds ratios
exp(coef(m))

> exp(coef(m))
    pared    public       gpa 
2.8510579 0.9429088 1.8513972 

순서척도 로지스틱 회귀분석에서 OR은 proportional OR이라고도 부른다.

• 다른 변수가 고정되어있을 때, pared=0일 때에 비하여 pared=1일 때, “very likely” vs “somewhat likely” or “unlikely”의 OR이 2.851이다.
  
• 다른 변수가 고정되어있을 때, pared=0일 때에 비하여 pared=1일 때, “very likely” or “somewhat likely” vs “unlikely”의 OR이 2.851이다.
• 이는 즉, pared=1일 때, 대학원 진학(apply)을 할 가능성이 높다는 것을 뜻한다.
• 같은 binary 변수인 public도 위와 마찬가지로 해석한다.
• 다른 변수가 고정되어있을 때, gpa가 1 증가할 때, “very likely” vs “somewhat likely” or “unlikely” Odds가 1.85배가 된다.
  
• 다른 변수가 고정되어있을 때, gpa가 1 증가할 때, “very likely” or “somewhat likely” vs “unlikely” Odds가 1.85배가 된다.

이를 식으로 표현하면

log Odds(Y=very) = 2.2039 + 1.04pared - 0.05879public + 0.61594gpa

log Odds(Y=very or somewhat) = 4.2994 + 1.04pared - 0.05879public + 0.61594gpa

• 이 때, Odds(Y=1) = P(P1/(1-P1)), Odds(Y=1,2) = P((P1+P2)/(1-P1-P2))
• log Odds 는 logit function 이라고도 부른다.
  

위 두 방정식을 풀면 P(Y=1), P(Y=2), P(Y=3) 3개의 식을 얻을 수 있다. 이 식은 X가 주어졌을 때, Y의 예측값을 구한다.