Ordered Response Models

이전 장에서 불러올 데이터

library("nnet")
data("BankWages", package = "AER")

BankWages$edcat <- factor(BankWages$education)
levels(BankWages$edcat)[3:10] <- rep(c("14-15", "16-18", "19-21"), c(2, 3, 3))

bw_tab <- xtabs(~ gender + minority + edcat + job, data = BankWages)

bwmale <- subset(BankWages, gender == "male")
bw_mnl <- multinom(job ~ education + minority,
                   data = bwmale, trace = FALSE)

6. Ordered Response Models

6.2. Standard Ordered Response Models

데이터 탐색을 위해 분할표를 만든다.

bw_tab <- xtabs(~minority + edcat + job, data = bwmale)
ftable(bw_tab, col.vars = "edcat")

                   edcat  8 12 14-15 16-18 19-21
minority job                                    
no       custodial        6  8     0     0     0
         admin            6 30    66     7     1
         manage           0  0     4    38    28
yes      custodial        7  5     1     0     0
         admin            4 18    18     7     0
         manage           0  1     0     2     1

모든 변수의 공동 분포를 사용하여 데이터를 시각화해본다.

par(mfrow = c(2,2))
mosaicplot(bw_tab, off = c(5, 0, 0),
           dir = c("v", "v", "h"), col = rev(gray.colors(3)), main = "")

mosaicplot(bw_tab[1,,], off = 0, col = rev(gray.colors(3)), main = "minority = no")

mosaicplot(bw_tab[2,,], off = 0, col = rev(gray.colors(3)), main = "minority = yes")

MASS 패키지의 polr() 함수는 logit 또는 probit link가 있는 모델을 제공한다. clm()은 동일한 모델뿐만 아니라 추가 링크 함수 및 구조화된 임계값과 같은 더 많은 유연성을 제공한다. 또한 clmm()을 사용하면 무작위 효과를 포함할 수 있다.

library(MASS)

bw_olm <- polr(job ~ education + minority, data = bwmale, Hess = TRUE)
summary(bw_olm)

Call:
polr(formula = job ~ education + minority, data = bwmale, Hess = TRUE)

Coefficients:
             Value Std. Error t value
education    0.870    0.09307   9.348
minorityyes -1.056    0.41199  -2.564

Intercepts:
                Value   Std. Error t value
custodial|admin  7.9514  1.0769     7.3833
admin|manage    14.1721  1.4744     9.6124

Residual Deviance: 260.6396 
AIC: 268.6396 

education 및 minority 변수로 설명되는 job category 모델

결과에는 절편이 없고 education과 minority의 계수만 포함되며, 반드시 순서가 지정된 custodial | admin 및 admin | management 사이의 임계값만 포함된다. job category의 잠재 척도에서 education이 긍정적, minority가 부정적 영향을 미친다는 것을 알 수 있다.

계수 테스트 및 우도비 테스트:

coefficients(bw_olm)

  education minorityyes 
  0.8699976  -1.0564379 

library(lmtest)
lrtest(bw_olm)

Likelihood ratio test

Model 1: job ~ education + minority
Model 2: job ~ 1
  #Df  LogLik Df  Chisq Pr(>Chisq)    
1   4 -130.32                         
2   2 -231.34 -2 202.05  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ordered probit 모델은 다음 코드를 통해 생성할 수 있다.

bw_opm <- polr(job ~ education + minority, data = bwmale,
               Hess = TRUE, method = "probit")

# function() 함수를 사용하는 방법도 있다.
# bw_opm <- update(bw_old, method = "probit")

결과적으로 로짓과 프로빗 모델은 유사한 결과를 보이고 있다.

library(modelsummary)

modelsummary(list("Probit model" = bw_opm,
                  "Logit model" = bw_olm), fmt = 3)

	Probit model	Logit model
custodial|admin	4.443	7.951
	(0.557)	(1.077)
admin|manage	7.844	14.172
	(0.744)	(1.474)
education	0.479	0.870
	(0.047)	(0.093)
minorityyes	-0.509	-1.056
	(0.214)	(0.412)
Num.Obs.	258	258
AIC	270.4	268.6
BIC	284.6	282.9
RMSE	1.97	1.97

modelsummary()를 통한 정리

모델의 적합도를 평가하기 위해 혼동 행렬을 계산해볼 수 있다.

table(true = bwmale$job, pred = predict(bw_olm, type = "class"))

           pred
true        custodial admin manage
  custodial        13    14      0
  admin            10   144      3
  manage            0    31     43

predict(bw_olm, newdata = bwmale[25, ], type = "probs")

 custodial      admin     manage 
0.01725637 0.88105484 0.10168879 

predict(bw_olm, newdata = bwmale[25, ], type = "class")

[1] admin
Levels: custodial admin manage

6.2.3. Parallel regression assumption

multinomial model과 ordered logit/probit models를 비교하기
data: bank wages(앞의 예시와 동일)

library(nnet)
bw_mnl <- multinom(job ~ education + minority, data = bwmale, trace = FALSE)
AIC(bw_mnl, bw_olm, bw_opm)

AIC(bw_mnl, bw_olm, bw_opm, k = log(nrow(bwmale)))

multinomial model의 AIC값이 가장 낮고, ordered logit/probit 모델들은 비슷하다.
또한 BIC를 계산하기 위해 로그에 대한 패널티를 증가시켜도 multinomial 모델이 가장 잘 수행된다.

library(effects)

eff_mnl <- allEffects(bw_mnl, xlevels = list(education = 50))
eff_olm <- allEffects(bw_olm, xlevels = list(education = 50))

plot(eff_olm, multiline = TRUE)

education 및 minority 효과 플롯 - Ordered Model

plot(eff_mnl, multiline =TRUE)

education 및 minority 효과 플롯 - Multinomial Model

두 모델 모두 admin과 management 사이의 전환이 날카롭다(기울기가 크다?).
즉, 일정 수준 이상의 교육을 거치면 management로 일하게 될 가능성이 급격히 높아진다. 이 두 범주 사이에는 상대적으로 명확한 education의 효과가 있지만, custodial과 admin 사이에는 동일하게 말할 수 없다. custodial 근로자는 더 많은 교육을 받은 후에 admin 범주로 승진하는 것이 아니라 단순히 다른 궤적일 뿐이다. 그러나 proportional odds model(전자) 에서는 custodial과 admin 두 곡선이 평행하고 서로 다른 임계값에서 서로 다른 효과를 가질 가능성이 없으므로 이 경우 multinomial model(후자)이 더 적합하다.

정확도를 평가하기 위한 분할표: multinomial 및 ordered model 예측

table(true = bwmale$job, pred = predict(bw_mnl, type = "class"))

           pred
true        custodial admin manage
  custodial        13    14      0
  admin            10   138      9
  manage            0     7     67

table(true = bwmale$job, pred = predict(bw_olm, type = "class"))

           pred
true        custodial admin manage
  custodial        13    14      0
  admin            10   144      3
  manage            0    31     43

m-1 = 2 인 binary logit model과 ordered model을 비교하기

m-1 = 2 이므로 하나의 임계값에 대해 하나씩 binary logit model을 두어 추정할 것이다.
각 범주에 대해 서로 다른 절편을 얻고 각 임계값에 대해 서로 다른 기울기 매개변수도 얻는다.
모델은 custodial vs “not custodial” 응답에 대한 하나의 이진 모델과, manage vs “not manage” 응답에 대한 하나의 이진 모델을 설정했다.

bw_logit1 <- glm(I(job != "custodial") ~ education + minority,
                 family = binomial, data = bwmale)
bw_logit2 <- glm(I(job == "manage") ~ education + minority,
                 family = binomial, data = bwmale)

modelsummary(list("Custodial/Not custodial" = bw_logit1, "Manage/Not manage"  = bw_logit2),
             fmt = 3, estimate = "{estimate}{stars}")

	Custodial/Not custodial	Manage/Not manage
(Intercept)	-5.062***	-26.215***
	(1.148)	(4.312)
education	0.586***	1.645***
	(0.095)	(0.277)
minorityyes	-0.481	-2.120**
	(0.509)	(0.794)
Num.Obs.	258	258
AIC
BIC
Log.Lik.
F	20.052	18.155
RMSE	0.26	0.26

결론적으로, 이 경우 병렬 회귀 가정(Parallel regression assumption)이 위반된 것으로 보인다. 즉, POLR 및 ordered probit 모델이 이 데이터 세트에 적합하지 않다는 것이다. 앞서 설명한 바와 같이 “custodial” vs “admin”, “admin” vs “manage”를 각각 담당하는 잠재 프로세스가 상이하기 때문에 이러한 결과는 그럴듯하다.

6.3. Extensions

Ordered logit and probit 모델을 사용하는 데는 두 가지 문제가 있다.

1. 단일 지수와 parellel regression을 가정함
   

2. 역치값 매개변수 \(\alpha_j\)는 모든 individuals에 일정하다. 즉, 공변량에 무관하다.

따라서 이러한 문제를 극복하는 일반화가 고려되어야 한다.

6.3.1. Generalized Threshold Models

일반화된 임계값 모델에서 임계값 매개변수는 공변량의 선형함수이다.

\[\alpha_{ij} = \tilde{\alpha}_j + {x_i}^\tau r_j\]

확률 모델:

\(\pi_{ij} = H(\tilde{\alpha}_j + {x_i}^\tau \gamma_j - {x_i}^\tau \beta) - H(\tilde{\alpha}_{j-1} + {x_i}^\tau \gamma_{j-1} - {x_i}^\tau \beta)\)
\(= H(\tilde{\alpha}_j - {x_i}^\tau \beta_j) - H(\tilde{\alpha}_{j-1} - {x_i}^\tau \beta_{j-1})\)

위 식에서 \(\beta_j = \beta - \gamma_j\)는 \(\beta 와 \gamma_j\)로 따로 구분할 수 없다.

6.3.2. Sequential models

일부 순서 변수는 연속적으로만 도달할 수 있는 범주를 아우르기도 한다.
예를 들면, 대학 학위 데이터에서 학사는 석사 학위에 필요한 전제 조건이다.

그러면 전이확률 \(y_i = j\)는 적어도 이 범주가 \(y_i \le j\)에 도달했다는 점에서 모델링이 가능하다.

확률 모델:

\[P(y_i = j | y_i \ge j, x_i) = H(\alpha_j - {x_i}^\tau \beta)\]

다시, 범주별 매개변수 \(\beta_j\)를 사용하여 모델을 일반화할 수 있다.

★ There is a problem in that when the current code is compiled, the formula does not appear, and when the formula is compiled, the code does not appear. ★

Reference

https://discdown.org/microeconometrics/ordered-response-models.html