41. 順序ロジスティック回帰分析
Masahiko Asano
2025-09-02
1. 順序ロジスティック回帰分析とは
Ordered Logistic Regressionのこと
- 順序尺度（例：1=反対、2=中立、3=賛成）のようにカテゴリーに自然な順序はあるが間隔は等間隔とは限らないデータに使う手法
- Rでは主に MASS パッケージの
polr()関数がよく使われる
1.1 データ
- ここでは自民党総裁選に関する架空のデータ (df) を作成する
# -------------------------
# 1. データのシミュレーション
# -------------------------
N <- 500
age <- rnorm(N, mean = 50, sd = 15) # 年齢
female <- rbinom(N, 1, 0.5) # 女性ダミー
# 潜在スコア（線形予測子）
# 年齢が高いほど賛成↑、女性はやや反対寄り↓
latent <- 0.05*age - 0.5*female + rnorm(N)
# カテゴリ閾値を設定
cuts <- c(-1, 1)
# 潜在スコアを3カテゴリに切る (反対=1, 中立=2, 賛成=3)
attitude <- cut(latent,
breaks = c(-Inf, cuts, Inf),
labels = c("反対","中立","賛成"),
ordered_result = TRUE)
df <- data.frame(attitude, age, female)- 作成したデータフレーム (
df) は 5 変数、72 行
| attitude | 自民党の総裁選を実施するかどうか（賛成、中立、反対）→ 従属変数 |
| age | 政治家の年齢 |
| femele | 政治家の性別（女性なら 1、男性なら 0） |
1.2 順序ロジスティック回帰モデル式
- ここで想定する順序ロジスティック回帰モデル（累積ロジットモデル）は次のとおり
- 順序ロジスティック回帰（cumulative logit model）
一般式（累積ロジット）
線形予測子
\[η = \beta_{age} age + \beta_{female}female\]
2 つの累積ロジット（cutpoint は \(\alpha_1, \alpha_2\)）
\[logit[{P(Y≤j|X)}]= \theta_j - η,\]
\[logit[P(attitude ≦ 反対|age, female)]= \alpha_1 -η, \]
\[logit[P(attitude ≦ 中立|age, female)]= \alpha_2 -η, \]
- 累積ロジットはカテゴリ数が J の場合、cutpoint数は
J-1個
- ここでは従属変数 (attitude) が 3 カテゴリ（反対 < 中立 <
賛成）あり
→ cutpoint数は 2 個
→ 累積ロジット（cumulative logit）は 2 つ
- それぞれ「あるカテゴリ以下である確率」に対してロジットをとる
1 本目のロジット（cutpoint α₁）
→「反対 vs （中立 or 賛成）」を区切る境界
→ 反対以下である確率
\[logit[P(Y ≦ 反対)] = \alpha_1 -\beta_{age}age - \beta_{female}female\]
2 本目のロジット（cutpoint α₂）
→「（反対 or 中立） vs 賛成」を区切る境界
→ 反対〜中立である確率
\[logit[P(Y ≦ 中立)] = \alpha_2 -\beta_{age}age - \beta_{female}female\]
βの係数にデフォルトでマイナスがついている理由
・分析結果の解釈を容易にするため
分析結果であるβの係数が
・プラスなら「上のカテゴリに行きやすい」
・マイナスなら「下のカテゴリに留まりやすい」
と解釈できる
\(Y\) = 態度（総裁選開催に反対=1, 中立=2, 賛成=3）
\(j\) = カテゴリの閾値（ここでは「反対|中立」「中立|賛成」の2つ）
\(\alpha_j\) = カテゴリごとの切片 (threshold)
\(\beta\) = 説明変数の係数：＋なら上のカテゴリ（＝「賛成寄り」）になりやすい
順序ロジットでは、カテゴリが順序を持っていることを前提にしている
もし従属変数 (attitude) の値が単なる factor（unordered）だと「賛成、中立、反対」が「赤・青・黄」と同じ扱いになる
→ 順序関係が無視されてしまう
→ ここでは反対 < 中立 < 賛成という順序を明示的に与えてみる.
→mutate()関数を使ってでattitudeを「順序付き因子（ordered factor）」に作り直す
female もファクターに変換
1.3 順序ロジスティック回帰分析
- 順序ロジスティック回帰分析では、同じデータでも分析のたびに係数の値が異なる
その理由：
- 順序ロジスティック回帰では、閾値（
cutpoint,threshold）と係数が同時に推定されるため
→ サンプルごとの分布の偏りによってcutpointの位置が変動する
→ それが係数推定にも影響する
→ 係数そのものより「限界効果」や「予測確率の変化」で解釈する
- ここでは
log-odds、オッズ比、限界効果と予測確率、それぞれの解釈方法を解説する
1.4 結果の解釈（log-odds）
推定結果を累積ロジット式と線形予測子に当てはめてみる
ここで得られた係数値は log-odds
- \(\beta_{age} = 0.09111\)
- \(\beta_{female} = -1.12218\)
カット点:
- \(\theta_{反対|中立} =
-2.5149\)
- \(\theta_{中立|賛成} = 1.8682\)
累積ロジット式
\[logit[P(Y ≦ 反対)] = \alpha_1 +\beta_{age}age - \beta_{female}female\]
\[= -2.5149 - 0.09111・age + 1.12218・female\\→「反対\ vs \ （中立＋賛成）」のオッズのロジットを表す\]
\[logit[P(Y ≦ 中立)] = \alpha_2 +\beta_{age}age - \beta_{female}female\]
\[ = 1.8682 - 0.09111・age + 1.12218・female\\→「(反対＋中立)\ vs\ 賛成」のオッズのロジットを表す\]
線形予測子: η（イータ）
- 「線形予測子」とは「累積ロジット式の右辺で、説明変数をまとめる部分」のこと
- 線形予測子 η 自体は「ただの線形結合」で、まだ確率ではない
\[η = \beta_{age} age + \beta_{female}female\]
\[ η = 0.09111・age -1.12218・female\]
- ここで得られた線形予測子の
ageとfemaleの回帰係数はlog(odds)=log-odds=logit
この結果からわかること
まとめ
・age の係数 = 0.09111 → 0 より大きい
→ 年齢が高いほど上のカテゴリ（=「賛成」）に移りやすい
・female の係数 = -1.12218 → 0 より小さい
→ 女性は男性に比べて下のカテゴリ（=「反対」）に移りやすい
・これ以上の解釈はできない
1.5 結果の解釈（オッズ比）
log(odds)をexp()してオッズ比に変換すると「どれくらい高いカテゴリに行きやすいか」を直感的に説明できる
- 例えば、age の係数 (0.09111) をオッズ比に変換してみる
[1] 1.095389broomパッケージのbroom()関数を使って、ageとfemale両方の変数を同時にオッズ比に変換する
各係数の解釈
age のオッズ比が 1.10
→ 年齢が1歳上がるごとに
・「反対 → 中立・賛成に入るオッズ」
・「反対・中立 → 賛成に入るオッズ」
の両方が約 1.1 倍（＝10%増加）する、という意味
female のオッズ比が 0.326
→ 女性だと
・「反対 → 中立・賛成に入るオッズ」
・「反対・中立 → 賛成に入るオッズ」
の両方が約 0.33 倍（＝67%減少）する、という意味
まとめ
・年齢が 1 歳増えると「賛成」寄りになるオッズが約 10 %増加
→ 高齢ほど賛成しやすい
・女性は男性に比べて「賛成」寄りになるオッズが約 67 %低下
→ 女性は反対・中立寄り
1.6 結果の解釈（予測確率と限界効果）
- オッズ比を使った解釈に加えて、予測確率を示すと理解しやすい
- 線形予測子 η を使って「閾値 \(\alpha_k\)
」と比較することで、各カテゴリの確率を計算する
- ここでは、男女別に自民党候補者が総裁選開催に関して反対、中立、賛成する予測確率を可視化してみる
予測確率の推移（性別ごと）
# 必要パッケージ
library(MASS)
library(ordinal)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(marginaleffects)
# ========= 前処理 & 推定（既にfit/df_cleanがあればスキップ可） =========
options(contrasts = c("contr.treatment","contr.poly"))
df_clean <- df %>%
transmute(
attitude = factor(attitude, levels = c("反対","中立","賛成"), ordered = TRUE),
age = as.numeric(age),
female = factor(female, levels = c(0,1), labels = c("男性","女性"))
) %>%
filter(complete.cases(across(everything()))) %>%
droplevels()
stopifnot(all(table(df_clean$attitude) > 0), all(table(df_clean$female) > 0))
cats <- levels(df_clean$attitude)
fit <- try(polr(attitude ~ age + female, data = df_clean, Hess = TRUE), silent = TRUE)
if (inherits(fit, "try-error")) {
fit <- clm(attitude ~ age + female, data = df_clean, link = "logit")
fit_engine <- "clm"
} else {
fit_engine <- "polr"
}
# ========= 男女別の予測確率（年齢=平均に固定） =========
mean_age <- mean(df_clean$age, na.rm = TRUE)
newdat <- data.frame(age = mean_age, female = c("男性","女性"))
# 予測確率＋95%CI
preds <- predictions(
fit,
newdata = newdat,
type = if (inherits(fit,"clm")) "prob" else "probs"
)
# カテゴリ列の自動検出
detect_cat_col <- function(d, cats){
cand <- intersect(names(d),
c("category","response","response.level","response_level",
"outcome","term","level","type","group"))
for(nm in cand){
v <- as.character(d[[nm]])
if (all(v %in% cats)) return(nm)
}
chr <- names(d)[sapply(d, function(x) is.character(x) || is.factor(x))]
for(nm in chr){
v <- as.character(d[[nm]])
if (all(v %in% cats)) return(nm)
}
NA_character_
}
cat_col <- detect_cat_col(preds, cats)
if (is.na(cat_col)) stop("カテゴリ列が特定できませんでした（predictions 出力）。")
preds_df <- preds %>%
rename(category = !!cat_col) %>%
mutate(
category = factor(as.character(category), levels = cats, ordered = TRUE),
female = factor(female, levels = c("男性","女性")),
label = percent(estimate, accuracy = 0.1)
)
# ========= 女性−男性の確率差（カテゴリ別）と p値 =========
comp <- comparisons(
fit,
variables = list(female = c("男性","女性")), # 上−下 = 女性−男性
newdata = newdat,
type = if (inherits(fit,"clm")) "prob" else "probs"
)
cat_col2 <- detect_cat_col(comp, cats)
if (is.na(cat_col2)) stop("カテゴリ列が特定できませんでした（comparisons 出力）。")
comp_out <- comp %>%
rename(category = !!cat_col2) %>%
mutate(category = factor(as.character(category), levels = cats, ordered = TRUE)) %>%
select(category, estimate, conf.low, conf.high, p.value)
# p値表示用データ（各カテゴリのバーの上に置く）
p_labels <- preds_df %>%
group_by(category) %>%
summarise(y = min(0.98, max(estimate) + 0.06), .groups = "drop") %>%
left_join(comp_out, by = "category") %>%
mutate(p_text = paste0("p = ", formatC(p.value, format = "f", digits = 3)))
# ========= グラフ：棒＋95%CI＋数値ラベル＋p値（カテゴリごと上部） =========
res_gender_ci_p <- ggplot(preds_df, aes(x = category,
y = estimate,
fill = female)) +
geom_col(position = position_dodge(width = 0.6), width = 0.55) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
position = position_dodge(width = 0.6), width = 0.28) +
geom_text(aes(label = label),
position = position_dodge(width = 0.6), vjust = -2, size = 4) +
# 各カテゴリの中央に p値を注記（男女差：女性−男性）
geom_text(data = p_labels,
aes(x = category, y = y, label = p_text),
inherit.aes = FALSE, vjust = -2.5, size = 4.2) +
scale_fill_manual(values = c("男性"="#4575b4","女性"="#d73027")) +
scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(0, 1.05)) +
labs(
title = sprintf("男女別の予測確率（年齢 = 平均値に固定）", fit_engine),
x = "カテゴリ", y = "予測確率", fill = "性別"
) +
theme_minimal(base_size = 14) +
theme_bw(base_family = "HiraKakuProN-W3")
res_gender_ci_p
結果の解釈（まとめ）
横軸：態度のカテゴリ（反対 → 中立 → 賛成）
縦軸：そのカテゴリを選ぶ確率（%）
1.4
結果の解釈（log-odds）:
- Note 1: 係数は
log-odds
- Note 2:
female1の意味は「female の値は女性が 1、男性が 0」
・female の log-odds = -1.122
→ 0 より小さい
→ 女性は男性に比べて「反対」寄りになりやすい
→ 統計的に有意
1.5 結果の解釈（オッズ比）:
・female の odds比 = 0.326
→ 女性は男性に比べて「賛成」寄りになるオッズが約 67 %低下
→ 女性は反対・中立寄り
1.6 結果の解釈（予測確率と限界効果）：
賛成:
・男性の予測確率 (93.1%)
・女性の予測確率 (86.0%)
→ 女性ダミーの限界効果 = 7.1%ポイント
→ 差は大きく、p値も有意 (p = 0.006)
中立:
・男性の予測確率 (13.7%)
・女性の予測確率 (6.7%)
→ 女性ダミーの限界効果 = 7.0%ポイント
→ 差は大きく、p値も有意 (p = 0.006)
反対:
・男性の予測確率 (0.1%)
・女性の予測確率 (0.3%)
→ 差は小さく、p値は有意ではない (p = 0.207)