順序ロジスティック回帰モデルとは？

2023.03.07

要約

　順序ロジスティック回帰モデルは、目的変数Yが順位のつけられる離散的な変数を取る場合の回帰モデルの一種です。

　目的変数Yが順位のつけられる離散的な変数の例には、「とても満足＝4」「満足＝3」「不満＝2」「とても不満=1」があります。

全体像

（1）問題の構造

「順序ロジスティック回帰モデルとは何か？」という問いを設定します。この問いを

①アイデアは何か

①定義は何か

②定義はどう解釈できるか

③どう結果の予測に使えるか

④どう原因の説明に使えるか

に分解します。

（2）前提の選択

　前提として①を選択します。

（3）論点の選択

　論点として②③④を選択します。

（4）付録一覧

　冗長さを避けて、可読性を上げるために、以下の内容は付録に回します。

「Rコード」

前提

①順序ロジスティック回帰モデルのアイデアは何か

　順序ロジスティック回帰モデルは、次のアイデアで成り立っています。

　モデル構築の動機は、「優＝4」「良＝3」「可＝2」「不可=1」というような順序尺度の決定要因の理解です。そのために、ある閾値を超えるか超えないかで、順序尺度が決定されると考えるのです。例えば、50点未満は不可、80点以上は優というな閾値と順序尺度の関係を考えます。

　潜在変数モデルで考えると、順序ロジスティック回帰モデルはわかりやすいです。実際には観測されない潜在変数Y*を考え、潜在変数Y*が誤差項uがロジスティック分布に従う次のモデル

$$Y^*=\beta_0+\beta_1X+u$$

$$u|X　〜　ロジスティック分布$$

に従うとします。そして、いくつかの閾値を設定すれば、順序ロジスティック回帰モデルとなります。

　なお、ロジスティック分布の確率密度関数は、以下の通りです。標準正規分布より裾が重い分布です。

図4：確率密度関数

　具体例を挙げましょう。次の潜在変数モデル

$$Y^*=0.1+0.2 X +u$$

閾値は-3、-1、-2

をプロットすると、下図になります。

図2：ロジスティック回帰モデルの潜在変数Y*を縦軸にとった場合（例）

　この順序ロジスティック回帰モデルでは、潜在変数Y*と目的変数Yとの関係が

$$2≦Y^*ならば　Y=4$$

$$-1≦Y^*<2ならば　Y=3$$

$$-3≦Y^*＜-1ならば　Y=2$$

$$Y^*<-3ならば　Y=1$$

で決まります。実際に観測されるデータは目的変数なので、下図になります。

図3：図2の潜在変数Y*をもちいて、目的変数Yを生成した図

②順序ロジスティック回帰モデルの数式上の定義は何か

　順序ロジスティック回帰モデルは、目的変数Yが「順位のつけられる離散的な変数」しか取らないモデルです。ここでは「Y=1，2，3，4」の場合を考えます。

$$Yが4になる確率=p4$$

$$Yが3になる確率=p3$$

$$Yが2になる確率=p2$$

$$Yが1になる確率=p1$$

としましょう。順序ロジスティック回帰モデルでは、閾値μ1、μ2、μ3を用いて

$$p4=\frac{e^{\beta_0+\beta_1 X-\mu_3}}{1+e^{\beta_0+\beta_1 X-\mu_3}}$$

$$p3=\frac{e^{\beta_0+\beta_1 X-\mu_2}}{1+e^{\beta_0+\beta_1 X-\mu_2}}-p4$$

$$p2=\frac{e^{\beta_0+\beta_1 X-\mu_1}}{1+e^{\beta_0+\beta_1 X-\mu_1}}-p4-p3$$

$$p1=1-p4-p3-p2$$

というようにp4、p3、p2、p1が決定されます。

　例えば、パラメーターが

$$\beta_0=0.1,　\beta_1=0.2$$

$$\mu_1=-3,　\mu_2=-1,　\mu_3=2$$

の場合、p1〜p4は下図で表せます。「説明変数Xが-20のとき、Y=1になる確率は約0.7」「X=-10のとき、Y=2になる確率が約0.5で最も高い」「X=10のとき、Yは3か4で半々」といったことがわかります。

図1：確率p1、p2、p3、p4がXによって変化する様子

方法

結果

（１）Yの予測

　説明変数Xが決まったときの目的変数Yの期待値は

$$E(Y|X)=4\times p4+3\times p3+2\times p2+ p1$$

です。p1、p2・・・にそれぞれ代入してやると、計算式はだいぶ複雑になります。

　例えば、パラメーターが

$$\beta_0=0.1,　\beta_1=0.2,　\mu_1=-3,　\mu_2=-1,　\mu_3=2$$

のとき、Yの期待値は下図になります。

図5：期待値

　また、閾値μもとても大事です。潜在変数Y*が閾値を超えるか超えないかで、目的変数Yが決定されるからです。

（２）X→Yの効果

　Xが微小1単位大きいと、Yが平均的に大きい量を、限界効果といいます。ただ、単純に限界効果を定義するわけにはいきません。なぜなら、順序ロジスティック回帰では、Yは順序のみを意味し、3は 1の3倍というように量を意味しないからです。閾値の幅が等しくないことからもわかるでしょう。

　そこで、限界効果をY=1,2,3,4それぞれについて、つぎのように定義します。

$$Yが1のときの限界効果＝\frac{d P(Y=1|X)}{d X}$$

$$Yが2のときの限界効果＝\frac{d P(Y=2|X)}{d X}$$

$$Yが3のときの限界効果＝\frac{d P(Y=3|X)}{d X}$$

$$Yが4のときの限界効果＝\frac{d P(Y=4|X)}{d X}$$

　これの性質は、どれかの確率が増えれば、どれかの確率が減るので、限界効果の合計は0になることです。

　なお、この限界効果を因果関係と解釈するには、統計的因果推論の考え方が必要で留保すべきです。

考察

（1）結論

（2）妥当性評価

前提評価

方法評価

結論評価

（3）意義

付録：Rコード

（1）分析コード

#MASSパッケージの呼び出し　
library(MASS)
　
#順序ロジスティック回帰
polr(as.factor(y) ~ x1 + x2 + x3, data=data0,  method = c("logistic"))

（2）作図

　作図にはggplot2を用います。

#パッケージの呼び出し
library(ggplot2)

　図1：確率p1、p2、p3、p4がXによって変化する様子↓

#関数
p4 <- function(x) (exp(1)^(0.1+0.2*x-2))/(1+(exp(1)^(+0.1+0.2*x-2)))
p3 <- function(x) (exp(1)^(0.1+0.2*x+1))/(1+(exp(1)^(0.1+0.2*x+1)))-(exp(1)^(0.1+0.2*x-2))/(1+(exp(1)^(0.1+0.2*x-2)))
p2 <- function(x) (exp(1)^(+0.1+0.2*x+3))/(1+(exp(1)^(0.1+0.2*x+3)))-(exp(1)^(0.1+0.2*x+1))/(1+(exp(1)^(0.1+0.2*x+1)))
p1 <- function(x) 1-(exp(1)^(0.1+0.2*x+3))/(1+exp(1)^(0.1+0.2*x+3))
　
#描画
ggplot()+
  stat_function(fun=p1,color = "black")+
  stat_function(fun=p2,color = "red")+
  stat_function(fun=p3,color = "blue")+
  stat_function(fun=p4,color = "violet")+
  scale_x_continuous(limits = c(-30, 30))

　潜在変数と目的変数の作成↓

#サンプル数を500
sample_number <- 500
　
#sample_number個の正規分布にしたがうxを生成
x <- rnorm(n=sample_number, mean =-0.5, sd=10)
  
#sample_number個のロジスティック分布に従う乱数を発生
u <- rlogis(sample_number)
  
# yy = 0.1 + 0.2*x+ u　に従う潜在変数yyを作成
yy <- 0.1 + 0.2*x+ u
  
# yy≧2ならy=4、-1≦yy<2ならy=3、-3≦yy<-1ならy=2、yy<-3ならy=1と
y <- ifelse(yy>=2,4,
            ifelse(yy>=-1,3,
                   ifelse(yy>=-3,2,1)))
　
#データフレーム化
data0 <- data.frame(x,y,yy)

　図2：ロジスティック回帰モデルの潜在変数Y*を縦軸にとった場合（例）↓

#関数の作成
m3 <- function(x) 2
m2 <- function(x) -1
m1 <- function(x) -3
　
#描画
ggplot(data0,aes(x =x, y=yy))+
  geom_point()+
  stat_function(fun=m1)+
  stat_function(fun=m2)+
  stat_function(fun=m3)+
  scale_x_continuous(limits = c(-30, 30))+
  scale_y_continuous(limits = c(-7, 7))

　図3：図2の潜在変数Y*をもちいて、目的変数Yを生成した図↓

ggplot(data0,aes(x =x, y=y))+
  geom_point()+
  scale_x_continuous(limits = c(-30, 30))+
  scale_y_continuous(limits = c(0, 5))

　図4：確率密度関数↓

#確率密度関数を出力。横軸は-10から10
ggplot(data.frame(x = c(-10, 10)), aes(x))+
  stat_function(fun=dlogis)+ #ロジスティック分布
  stat_function(fun=dlogis)　#標準正規分布

図5

E <- function(x) (4*exp(1)^(0.1+0.2*x-2))/(1+(exp(1)^(+0.1+0.2*x-2)))+
    (3*exp(1)^(0.1+0.2*x+1))/(1+(exp(1)^(0.1+0.2*x+1)))-(3*exp(1)^(0.1+0.2*x-2))/(1+(exp(1)^(0.1+0.2*x-2)))+
    (2*exp(1)^(+0.1+0.2*x+3))/(1+(exp(1)^(0.1+0.2*x+3)))-(2*exp(1)^(0.1+0.2*x+1))/(1+(exp(1)^(0.1+0.2*x+1)))+
    1-(exp(1)^(0.1+0.2*x+3))/(1+exp(1)^(0.1+0.2*x+3))

#描画
ggplot()+
  stat_function(fun=E,color = "black")+
  scale_y_continuous(limits = c(0, 5))+
  scale_x_continuous(limits = c(-30, 30))