p <- seq(0, 1, by = .0001) # 0から1の確率ヴェクトルを作成
odds <- p/(1 - p) # オッズ
log_odds <- log(odds) # 対数オッズ『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』
第12章 ロジスティック回帰分析
0 全体の章構成
1 一般化線型モデル(Generalized Linear Model)
一般化線型モデル(GLM)の中で社会学では最もポピュラーだと思われる二項ロジスティック回帰分析について説明する。
- リンク関数(link function)と線型予測子(linear predictor)
\[g(\hat{y})=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{j}x_{j}+...+\beta_{p}x_{p}\]
1-1 ロジスティック関数
以下のリンク関数をロジット関数という。
\[g(\hat{y})=log\left(\frac{p(\hat{y}=1)}{1-p(\hat{y}=1)}\right)\]
線形予測子を連結すると以下の通りである。
\[g(\hat{y})=log\left(\frac{p(\hat{y}=1)}{1-p(\hat{y}=1)}\right)=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{j}x_{j}+...+\beta_{p}x_{p}\] 以下の様にオッズの式にすると，オッズを独立変数による効果の積として予測するものであることがわかる。
\[\frac{p(\hat{y}=1)}{1-p(\hat{y}=1)}=e^{\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{j}x_{j}+...+\beta_{p}x_{p}}=e^{\beta_{0}}e^{\beta_{1}x_{1}}e^{\beta_{2}x_{2}}...e^{\beta_{j}x_{j}}...e^{\beta_{p}x_{p}}\] また，確率の式にすると以下の様になるが，このロジット変換の逆関数を，ロジスティック関数という。
\[p(\hat{y}=1)=\frac{1}{1+e^{-(\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{j}x_{j}+...+\beta_{p}x_{p})}}\]
オッズ比と対数オッズ比
確率と，オッズ，対数オッズの関係を図示する。
作図スクリプト
plot(p, odds,
ylim = c(-15, 15),
xlim = c(0, 1),
type = "l",
lwd = 3,
main = "オッズ（赤の実線）と対数オッズ（青の破線）",
xlab = "確率p",
ylab = "オッズと対数オッズの値",
col = "red",
bty = "n")
abline(h = 0,
lty = "solid",
col = "gray")
par(new=T)
plot(p, log_odds,
ylim = c(-15, 15),
xlim = c(0, 1),
type = "l",
lty = "dashed",
lwd = 3,
ylab = "",
xlab = "",
col = "blue",
bty = "n")
線形予測子の取り得る値の範囲は -∞ から +∞ であるが，確率は0~1の値しかとらない。
オッズも，正の値しかとらない。
しかし対数オッズは， -∞ から +∞ の値を取り得る。
よって，この中で線形予測子と連結するのにふさわしいのは，対数オッズであるといえる。
1-2 最尤法(Maximum Liklihood Estimation)
男性10人，女性8人の標本で，内閣支持者がそれぞれ5人と2人になった場合の，母集団における支持率を最尤法で推定してみよう。
n1 <- 10 # 男性の標本サイズ
n2 <- 8 # 女性の標本サイズ
w1 <- 5 # 男性標本における支持者数
w2 <- 2 # 女性標本における支持者数
pie <- seq(0, 1, by=.01) # 母比率の仮定値のヴェクトル
# 母集団の仮定値に対応する尤度のヴェクトル
likelihood <- choose(n1, w1) * pie ^ w1 * (1 - pie) ^ (n1 - w1) *
choose(n2, w2) * pie ^ w2 * (1 - pie) ^ (n2 - w2)様々な仮の母比率に対応する尤度の作図
plot(pie, likelihood,
bty = "n",
type = "l",
lty = "dashed",
xlab = "母比率の仮定値",
ylab = "尤度",
main = "最尤法による母数の推定")
abline(v = 7/18,
col = "#00000080",
lty = "dotted")
尤度（母数がある値だと仮定した場合に，このような標本が出現する確率）が最大になる母比率は，7/18 ≒ .39（男女合わせて18人中7人）である。グラフではグレーの点線で垂直線を書き込んである。
男女で母比率が異なると仮定
ここで仮に，これまでの調査研究の結果から，母集団における内閣支持率は，女性よりも男性の方が15ポイント高いと推測されるとしよう1。その場合にこの標本の尤度が最大になる（男性の）母比率は幾つになるだろうか。同様にグラフを描きつつ調べてみよう。
m_pie <- seq(.15, 1, by = .000001) # 男性の母比率の仮定値ヴェクトル
likelihood <- choose(n1, w1) * m_pie ^ w1 * (1 - m_pie) ^ (n1 - w1) *
choose(n2, w2) * (m_pie - .15) ^ w2 * (1 - (m_pie - .15)) ^ (n2 - w2)
plot(m_pie, likelihood,
type = "l",
lty = "dashed",
xlab = "男性の母比率の仮定値",
ylab = "尤度",
main = "最尤法による母数の推定")
尤度が最大になる母比率の仮説値は次の様に求める事が出来る。
m_pie[likelihood == max(likelihood)][1] 0.451528男性の標本比率は 5/10 = .50 であるのに対して，推定された母比率は 約.45，女性の標本比率は 2/8 = .25 であるのに対して，推定された母比率は 約.30 と云う事になる。
尤度と対数尤度
尤度(0～1)と対数尤度(-∞～0)の関係をグラフで示す。
自然対数は，1に対して0となり，0以上1未満の値に対してはマイナスの値になる。
よって対数尤度は基本的にマイナスであり，最大で0になる。
但し，尤度が大きければ対数尤度も大きい（マイナスだが絶対値が小さくなる）と云う事は確かである。
L <- seq(0, 1, by = .0001) # 0から1の尤度
log_L <- log(L) # 対数尤度尤度と対数尤度の作図
plot(L, log_L,
ylim = c(-10, 0),
xlim = c(0, 1),
type = "l",
main = "尤度と対数尤度の関係",
xlab = "尤度",
ylab = "対数尤度の値",
bty = "n")
abline(h = 0,
lty = 3,
col = "grey")
対数尤度と逸脱度
尤度が大きければ対数尤度は大きくなり（マイナスで絶対値が小さくなり）-2対数尤度（逸脱度）はプラスで値が小さくなる。
“より大きな”モデル（独立変数がより多いモデル）と“より小さな”モデルの二つのモデルがあった時，逸脱度は“より大きな”モデルの方が小さい。
より小さなモデルの逸脱度とより大きなモデルの逸脱度の差がカイ二乗分布に従う時，その差が有意であれば（偶然とは言えないくらい異なれば），逸脱度が小さい，“より大きな”モデルを採用する。
逸脱度の差が有意でなければ（偶然の範囲内である事が否定できなければ），同程度に逸脱しているがより簡潔なモデル，すなわち“より小さな”モデルを採用する。
対数尤度と逸脱度の作図
plot(L, log_L,
ylim = c(-10, 10),
xlim = c(0, 1),
type = "l",
main = "対数尤度(LL)と-2対数尤度(-2LL)の関係",
xlab = "尤度L",
ylab = "LLと-2LLの値",
bty = "n",
lty = "dotted")
par(new = T)
plot(L, -2*log_L,
ylim = c(-10, 10),
xlim = c(0, 1),
type = "l",
main = "",
xlab = "",
ylab = "",
bty = "n")
text(.1, -5, "対数尤度(LL)",
adj = 0)
text(.1, 8, "-2対数尤度(-2LL)",
adj = 0)
1-3 対数オッズ比(log odds ratio)と2項ロジットモデル(binomial logit model)
glm( ) の実行
glm( )関数を用いて，二項分布binomial，ロジットリンクlogitを指定。
性別，年齢，教育年数，個人年収の4つの独立変数を投入する2。
result <- glm(single ~ SEX + age + edu + income,
data = d02,
family = binomial(link = "logit"))
summary(result)
Call:
glm(formula = single ~ SEX + age + edu + income, family = binomial(link = "logit"),
data = d02)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.8769972 1.4952433 2.593 0.00952 **
SEX女性 -1.7164462 0.3807364 -4.508 6.54e-06 ***
age -0.0623619 0.0190020 -3.282 0.00103 **
edu -0.0569334 0.0730136 -0.780 0.43553
income -0.0023085 0.0007374 -3.130 0.00175 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 344.82 on 337 degrees of freedom
Residual deviance: 306.77 on 333 degrees of freedom
AIC: 316.77
Number of Fisher Scoring iterations: 5Coefficientsの欄は，重回帰分析の結果表とほぼ同じ様に理解出来る。
- 各独立変数について，Estimate（対数オッズ比）が正であれば従属変数が1となる確率を高め，負であれば1となる確率を低める。
- Estimate を Std. Error で割ったものが z value であり，zと名前がついている様に，これは近似的に標準正規分布に従うとされ，標準正規分布を用いたゼロ仮説有意性検定の結果の有意確率が Pr(>|z|) である。両側検定の確率なので，正負両方で絶対値がz value以上になる確率が表示される。
- Null devianceはゼロ逸脱度，Residual devianceは残差逸脱度である。
- AIC = -2LL + 2k = deviance + 2k で，kは独立変数の個数だけを数える流儀と，切片も含める場合がある。上記では切片を含めてk = 5となっている。
対数オッズ比をオッズ比に
計算結果は，各独立変数の係数，つまり対数オッズ比で示されるが，これをオッズ比に変換して示す事は多い。
exp(result$coefficients)(Intercept) SEX女性 age edu income
48.2790265 0.1797036 0.9395428 0.9446570 0.9976942 結果を並べて見たければ次のようにする。
1列目が対数オッズ比，2列目がオッズ比である。
読み取りやすさの為に，同時に四捨五入して表示する。
round(cbind(result$coefficients, exp(result$coefficients)), 3) [,1] [,2]
(Intercept) 3.877 48.279
SEX女性 -1.716 0.180
age -0.062 0.940
edu -0.057 0.945
income -0.002 0.998信頼区間の計算
係数（対数オッズ比）の信頼区間，オッズ比の信頼区間を四捨五入して表示。
round(confint(result), 3) 2.5 % 97.5 %
(Intercept) 0.981 6.864
SEX女性 -2.483 -0.987
age -0.101 -0.026
edu -0.200 0.087
income -0.004 -0.001round(exp(confint(result)), 3) 2.5 % 97.5 %
(Intercept) 2.667 957.501
SEX女性 0.083 0.373
age 0.904 0.975
edu 0.818 1.091
income 0.996 0.999平均限界効果(AME)
ロジスティック回帰分析の結果は，オッズが何倍になるかというオッズ比で示されるが，これは直感的には必ずしも分かり易くない。
限界効果(marginal effect)は，ある独立変数が1単位増加した時に，従属変数である「y=1となる確率」がどれだけ変化するかを意味し，独立変数が数量変数である場合には限界的な変化率（極限的な傾き，微分）となる。
AMEは，各ケースにおいて計算されたMEの平均である。
比較的新しいパッケイジ。
marginaleffects::avg_slopes(result)
Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 %
SEX 女性 - 男性 -0.270621 0.058429 -4.632 < 0.001 18.1 -0.385139
age dY/dX -0.009055 0.002641 -3.429 < 0.001 10.7 -0.014231
edu dY/dX -0.008267 0.010581 -0.781 0.43466 1.2 -0.029006
income dY/dX -0.000335 0.000103 -3.260 0.00112 9.8 -0.000537
97.5 %
-0.156103
-0.003879
0.012472
-0.000134
Type: response広く普及しているがやや古典的なパッケイジ。
AME <- margins::margins(result)
summary(AME)AMEを見ても女性だと未婚になりにくい効果が強いので，単純にデータで確認しておく。
これは他の変数を全く統制していない結果なので，平均限界効果などと直ちに一致するものではない。
t_Single <- with(d02, table(SEX, single))
addmargins(t_Single) single
SEX 0 1 Sum
男性 98 41 139
女性 170 29 199
Sum 268 70 338round(prop.table(t_Single, margin = 1)*100, 1) single
SEX 0 1
男性 70.5 29.5
女性 85.4 14.6女性の未婚者割合は15%弱なのに対し男性では30%近くでちょうど2倍程度である。
1-4 適合度検定(Goodness of Fit Test)と疑似決定係数(Pseudo R-squared)
Hosmer-Lemeshowの適合度検定
二項ロジスティック回帰分析の適合度検定であるHosmer-Lemeshow検定をRで実行するには，2026年2月現在，Web上で検索すると，以下の様なパッケイジと関数が見つかる 6。
- glmtoolbox::hltest( )
- ResourceSelection::hoslem.test( )
- performance::performance_hosmer( )
- generalhoslem::logitgof( )
ここでは，上記のうち，テクスト本文で使用している2と，この検定専用のパッケイジである4を試してみよう。
ResourceSelection::hoslem.test( )
ResourceSelection::hoslem.test(result$y, result$fitted.values, g = 10)
Hosmer and Lemeshow goodness of fit (GOF) test
data: result$y, result$fitted.values
X-squared = 13.46, df = 8, p-value = 0.09698g = で，いくつのグループに分割して比較するかを指定する。通常は10であるが，これを変えると検定結果も変わる。
この例では，グループを9分割，10分割，11分割と変えると，5%有意，10%有意，1%有意とかなり変動する。
ResourceSelection::hoslem.test(result$y, result$fitted.values, g = 9)
Hosmer and Lemeshow goodness of fit (GOF) test
data: result$y, result$fitted.values
X-squared = 17.437, df = 7, p-value = 0.01479ResourceSelection::hoslem.test(result$y, result$fitted.values, g = 11)
Hosmer and Lemeshow goodness of fit (GOF) test
data: result$y, result$fitted.values
X-squared = 22.277, df = 9, p-value = 0.008041著名な統計学者のPaul AllisonはHosmer-Lemeshow Testを信用しないと書いている［2026年2月26日閲覧］。
generalhoslem::logitgof( )
generalhoslem::logitgof(result$y, result$fitted.values, g = 10, ord = FALSE)
Hosmer and Lemeshow test (binary model)
data: result$y, result$fitted.values
X-squared = 13.46, df = 8, p-value = 0.09698結果も情報量も ResourceSelection::hoslem.test( ) と同じである。
擬似決定係数
最尤法を用いた一般線形モデルでは，最小二乗法を用いた線形モデルのように分散説明率を説明力の指標とはしない。
その代わりに，尤度を元にした擬似決定係数(Pseudo R-squared)を計算する事が多い。
但しこれはあくまで疑似（似て非なるもの）であって本来の決定係数（分散説明率）とは異なる7ので，工夫の仕方によって幾つもの種類が提案されている。
現在検索すると，擬似決定係数の計算にはDescTools パッケイジの PseudoR2( ) が使用可能と出る。
DescTools::PseudoR2(result, "all") McFadden McFaddenAdj CoxSnell Nagelkerke AldrichNelson
0.11034136 0.08134078 0.10646341 0.16648700 0.10117853
VeallZimmermann Efron McKelveyZavoina Tjur AIC
0.20035571 0.13461027 0.22674936 0.12519257 316.77268914
BIC logLik logLik0 G2
335.88791862 -153.38634457 -172.41033495 38.04798077 非常に多くの擬似決定係数が表示されているが，一般的に知られているのはMcFadden，CoxSnell，Nagelkerkeくらいで，最も目にする様に思うのは何故かNagelkerkeである。
疑似決定係数の定義
各種の疑似決定係数の定義式は次のサイトを見ると良い。
UCLA’s Office of Advanced Research Computing, ‘FAQ: What are pseudo R-squareds?’（2026年2月26日閲覧）
LogLiklihoodをLLと表記すると，
McFadden’s pseudo R-squared: \[R^2=1-\frac{LL_{model}}{LL_{null}}\]
Cox-Snell’s pseudo R-squared: \[R^2=1-\left(\frac{L_{null}}{L_{model}}\right)^{\frac{2}{n}}\]
Nagelkerke’s pseudo R-squared: \[R^2=\frac{1-\left(\frac{L_{null}}{L_{model}}\right)^{\frac{2}{n}}}{1-{L_{null}}^{\frac{2}{n}}}\]
追加パッケイジ無しで計算できる
主なPseudo R2を，追加パッケイジをインストールする事無く，自分で計算してみよう。
スクリプトの汎用性を高める為に， 最初にglm( )の結果オブジェクトを glm_result にcopyする。
glm_result <- result # この右辺のresultを適宜書き換える。
LLm <- glm_result$deviance/(-2)
LL0 <- glm_result$null.deviance/(-2)
McFadden <- 1 - LLm/LL0 # McFadden
adjMcFadden <- 1 - (LLm - glm_result$rank)/LL0 # adjusted McFadden
# Cox & Snell
n_model <- glm_result$rank + glm_result$df.residual; # n_model
CoxSnell <- 1-(exp((2/n_model)*(LL0 - LLm)))
# Nagelkerke
Nagelkerke <- (1-(exp((2/n_model)*(LL0 - LLm))))/(1-exp((2/n_model)*LL0)) McFadden[1] 0.1103414adjMcFadden[1] 0.08134078CoxSnell[1] 0.1064634Nagelkerke[1] 0.166487DescToolsのMcFaddenAdjとBaylorEdPsychのAdj.McFaddenは一致しないが，それぞれがどの様な計算なのかを探索してみる。
1 - (LLm - glm_result$rank )/LL0 # DescToolsのMcFaddenAdjに一致[1] 0.081340781 - (LLm - glm_result$rank - 1)/LL0 # BaylorEdPsychのAdj.McFaddenに一致[1] 0.075540662 ロジットモデル
2-1 多項ロジット(Multinomial Logit Model)
nnetパッケイジを使用して，婚姻上の地位（未婚，有配偶，離死別）を性別，年齢，教育年数，本人年収で説明する多項ロジット分析を行う。
従属変数の作成
d01$marital <- c(2, 2, 1, 3, 3)[d01$q2100] # 婚姻上の地位: 有配偶2，未婚1，離死別3
d01$Marital <- factor(d01$marital, levels = 1:3,
labels = c("unmarried", "married", "separated"))
d01$Single <- relevel(d01$Marital, ref = "married")完備ケースデータの準備
full <- with(d01, complete.cases(Single, SEX, age, edu, income))
d04 <- d01[full, c("Single", "SEX", "age", "edu", "income")]多項ロジット分析の実行(nnet)
nnetパッケイジの multinom( ) 関数を使用する。
mres <- nnet::multinom(Single ~ SEX + age + edu + income, d04)# weights: 18 (10 variable)
initial value 371.330954
iter 10 value 235.619548
iter 20 value 233.657740
iter 20 value 233.657740
iter 20 value 233.657740
final value 233.657740
convergedmsum <- summary(mres)
msumCall:
nnet::multinom(formula = Single ~ SEX + age + edu + income, data = d04)
Coefficients:
(Intercept) SEX女性 age edu income
unmarried 4.054518 -1.613429 -0.06040090 -0.0745727 -0.0022640068
separated -1.333564 1.622518 0.03039848 -0.2799723 0.0007676712
Std. Errors:
(Intercept) SEX女性 age edu income
unmarried 1.507523 0.3829439 0.01912610 0.07381118 0.0007384185
separated 2.222100 0.6526882 0.02861961 0.12324943 0.0007255267
Residual Deviance: 467.3155
AIC: 487.3155 AIC = Deviance + 2k，kは推定されている係数の数の10個になっている。
パラメタ推定値のz検定は以下の様にして行える。
z <- msum$coefficients/msum$standard.errors
p <- pnorm(abs(z), mean = 0, sd = 1, lower.tail = F)*2検定統計量を小数第3位までで表示
round(z, 3) (Intercept) SEX女性 age edu income
unmarried 2.69 -4.213 -3.158 -1.010 -3.066
separated -0.60 2.486 1.062 -2.272 1.058有意確率を小数第3位までで表示
round(p, 3) (Intercept) SEX女性 age edu income
unmarried 0.007 0.000 0.002 0.312 0.002
separated 0.548 0.013 0.288 0.023 0.290mlogitパッケイジによる多項ロジット
このnnetのmultinom( )関数では，係数のz検定の結果も，疑似決定係数も出力しないのでやや不親切にも感じる。
mlogitパッケイジのmlogit( )関数はもっと多くの情報を出力するので便利に思えるが，データの形式を変換してから分析を行わなければならないので却って面倒な部分もある。
しかも，z検定や疑似決定係数は，multinom( )の情報を使って自分で計算する事は可能である。どちらのパッケイジ，どちらの関数の方が自分にとって簡単であるかによって判断しよう。
以下では，まずmlogitパッケイジのmlogit( )関数による分析結果を示したのち，最初に紹介したmultinom( )関数の情報を使って同様の結果を得る方法を紹介する。
データの変換
d05 <- mlogit::mlogit.data(d04, shape = "wide", choice = "Single")多項ロジット分析
mres2 <- mlogit::mlogit(Single ~ 0|SEX + age + edu + income|0, d05)
summary(mres2)
Call:
mlogit::mlogit(formula = Single ~ 0 | SEX + age + edu + income |
0, data = d05, method = "nr")
Frequencies of alternatives:choice
married separated unmarried
0.713018 0.079882 0.207101
nr method
6 iterations, 0h:0m:0s
g'(-H)^-1g = 2.98E-08
gradient close to zero
Coefficients :
Estimate Std. Error z-value Pr(>|z|)
(Intercept):separated -1.33353834 2.22210282 -0.6001 0.548423
(Intercept):unmarried 4.05461164 1.50752321 2.6896 0.007154 **
SEX女性:separated 1.62247261 0.65268804 2.4858 0.012925 *
SEX女性:unmarried -1.61339867 0.38294221 -4.2132 2.518e-05 ***
age:separated 0.03040280 0.02861985 1.0623 0.288101
age:unmarried -0.06040214 0.01912612 -3.1581 0.001588 **
edu:separated -0.27998519 0.12324954 -2.2717 0.023105 *
edu:unmarried -0.07457828 0.07381115 -1.0104 0.312307
income:separated 0.00076758 0.00072555 1.0579 0.290088
income:unmarried -0.00226392 0.00073841 -3.0659 0.002170 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log-Likelihood: -233.66
McFadden R^2: 0.10122
Likelihood ratio test : chisq = 52.627 (p.value = 1.2734e-08)- 係数の推定値や標準誤差推定値はmultinom( )関数の結果と，誤差の範囲で一致していると言える。
- 検定統計量Zや有意確率もほぼ一致しているが，こちらは自分で計算しなくても見易い一覧表で示してくれる。
- McFaddenの疑似決定係数や尤度比検定の結果も，multinom( )関数の結果には表示されていなかったものである。
multinom( )の結果から疑似決定係数を計算
以下では，multinom( )関数の結果を用いて，疑似決定係数や尤度比検定を行う方法を示す。
これが出来れば，multinom( )関数の代わりに，データ形式の変換が必要なmlogit( )関数を使う必要は必ずしもないとも言える。
分析モデルとnullモデルの対数尤度を計算
LLm <- logLik(nnet::multinom(Single ~ SEX + age + edu + income, d04))
LL0 <- logLik(nnet::multinom(Single ~ 1, d04))LLm; LL0'log Lik.' -233.6577 (df=10)'log Lik.' -259.9712 (df=2)multinom( )の結果から疑似決定係数
# McFadden
1-LLm[1]/LL0[1] [1] 0.1012168# Cox & Snell
n <- length(d04$Single)
1-(exp(LL0[1])/exp(LLm[1]))^(2/n) [1] 0.1441849# Nagelkerke
(1-(exp(LL0[1])/exp(LLm[1]))^(2/n))/(1-exp(LL0[1])^(2/n)) [1] 0.183616モデル全体の尤度比検定
null <- summary(nnet::multinom(Single ~ 1, d04)) # 切片だけのモデル# weights: 6 (2 variable)
initial value 371.330954
iter 10 value 259.971189
iter 10 value 259.971189
iter 10 value 259.971189
final value 259.971189
convergednull$deviance - msum$deviance[1] 52.6269pchisq(null$deviance - msum$deviance, df = 2*(msum$rank - null$rank), lower.tail = F)[1] 1.273438e-08mlogit( )の結果と見比べて，対数尤度の値，McFaddenの疑似決定係数の値，尤度比検定の結果も全て一致している。
更に，Cox & Snellなど他の疑似決定係数の値も計算出来る。
2つのパッケイジの結果の比較
改めて両方の結果を並べて見よう。
summary(mres2)
Call:
mlogit::mlogit(formula = Single ~ 0 | SEX + age + edu + income |
0, data = d05, method = "nr")
Frequencies of alternatives:choice
married separated unmarried
0.713018 0.079882 0.207101
nr method
6 iterations, 0h:0m:0s
g'(-H)^-1g = 2.98E-08
gradient close to zero
Coefficients :
Estimate Std. Error z-value Pr(>|z|)
(Intercept):separated -1.33353834 2.22210282 -0.6001 0.548423
(Intercept):unmarried 4.05461164 1.50752321 2.6896 0.007154 **
SEX女性:separated 1.62247261 0.65268804 2.4858 0.012925 *
SEX女性:unmarried -1.61339867 0.38294221 -4.2132 2.518e-05 ***
age:separated 0.03040280 0.02861985 1.0623 0.288101
age:unmarried -0.06040214 0.01912612 -3.1581 0.001588 **
edu:separated -0.27998519 0.12324954 -2.2717 0.023105 *
edu:unmarried -0.07457828 0.07381115 -1.0104 0.312307
income:separated 0.00076758 0.00072555 1.0579 0.290088
income:unmarried -0.00226392 0.00073841 -3.0659 0.002170 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log-Likelihood: -233.66
McFadden R^2: 0.10122
Likelihood ratio test : chisq = 52.627 (p.value = 1.2734e-08)結果表示のスクリプト
R.sq <-
sprintf("McFadden R2: %.3f\nCox & Snell R2: %.3f\nNagelkerke R2: %.3f",
1-LLm[1]/LL0[1],
1-(exp(LL0[1])/exp(LLm[1]))^(2/n) ,
(1-(exp(LL0[1])/exp(LLm[1]))^(2/n))/(1-exp(LL0[1])^(2/n)) )
names(R.sq) <- c("pseudo R-squared")
LR.test <-
sprintf("chi-squared: %.3f, p.value =%.5f",
null$deviance - msum$deviance,
pchisq(null$deviance - msum$deviance, df = 2*(msum$rank - null$rank), lower.tail = F))
names(LR.test) <- c("Likelihood ratio test")
t(round(first_row, 4))
t(round(second_row, 4))
R.sq
LR.test stats
unmarried coefficient std.error z-value Pr(>|z|)
(Intercept) 4.0545 1.5075 2.6895 0.0072
SEX女性 -1.6134 0.3829 -4.2132 0.0000
age -0.0604 0.0191 -3.1580 0.0016
edu -0.0746 0.0738 -1.0103 0.3123
income -0.0023 0.0007 -3.0660 0.0022
stats
separated coefficient std.error z-value Pr(>|z|)
(Intercept) -1.3336 2.2221 -0.6001 0.5484
SEX女性 1.6225 0.6527 2.4859 0.0129
age 0.0304 0.0286 1.0622 0.2882
edu -0.2800 0.1232 -2.2716 0.0231
income 0.0008 0.0007 1.0581 0.2900
pseudo R-squared
"McFadden R2: 0.101\nCox & Snell R2: 0.144\nNagelkerke R2: 0.184"
Likelihood ratio test
"chi-squared: 52.627, p.value =0.00000" 2-2 順序ロジット(Ordinal Logit Model)
まずはMASSパッケイジのpolr( )関数，次にordinalパッケイジのclm( )関数を紹介する。
MASSパッケイジのpolr( )関数
ordinal <- MASS::polr(STRATA ~ SEX + age + edu + JOB + fincome,
data = d06,
method = "logistic")
summary(ordinal)
Re-fitting to get HessianCall:
MASS::polr(formula = STRATA ~ SEX + age + edu + JOB + fincome,
data = d06, method = "logistic")
Coefficients:
Value Std. Error t value
SEX女性 0.049224 0.2865274 0.1718
age -0.017447 0.0119039 -1.4657
edu -0.093406 0.0448983 -2.0804
JOB非正規 0.246047 0.3393981 0.7250
JOB自営 -0.377164 0.4410765 -0.8551
JOB無職・学生 -0.341299 0.3738658 -0.9129
fincome -0.002183 0.0003204 -6.8138
Intercepts:
Value Std. Error t value
上|中の上 -8.7693 0.0248 -353.9837
中の上|中の下 -4.3160 0.5063 -8.5251
中の下|下の上 -2.0492 0.5253 -3.9013
下の上|下の下 0.3243 0.6246 0.5191
Residual Deviance: 621.2733
AIC: 643.2733 有意確率は表示されないが，t統計量の絶対値から判断して，世帯収入は明確に階層帰属意識を高め8，教育年数がギリギリ有意になっているだけである。
ordinalパッケイジのclm( )
比較しやすい様に独立変数を減らした同一モデルの分析結果を比較。
polr01 <- MASS::polr(STRATA ~ SEX + JOB + fincome,
data = d06,
method = "logistic")
clm01 <- ordinal::clm(STRATA ~ SEX + JOB + fincome,
data = d06)summary(polr01)
Re-fitting to get HessianCall:
MASS::polr(formula = STRATA ~ SEX + JOB + fincome, data = d06,
method = "logistic")
Coefficients:
Value Std. Error t value
SEX女性 0.082948 0.2869116 0.2891
JOB非正規 0.236262 0.3328376 0.7098
JOB自営 -0.307219 0.4401518 -0.6980
JOB無職・学生 -0.338396 0.3681615 -0.9192
fincome -0.002325 0.0003138 -7.4100
Intercepts:
Value Std. Error t value
上|中の上 -6.7394 0.7056 -9.5514
中の上|中の下 -2.2853 0.3072 -7.4386
中の下|下の上 -0.0367 0.2813 -0.1304
下の上|下の下 2.3336 0.4279 5.4533
Residual Deviance: 624.3366
AIC: 642.3366 t分布に近似してt検定を行っている。
summary(clm01)formula: STRATA ~ SEX + JOB + fincome
data: d06
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 291 -312.17 642.34 8(0) 5.46e-12 3.2e+07
Coefficients:
Estimate Std. Error z value Pr(>|z|)
SEX女性 0.0826114 0.2868123 0.288 0.773
JOB非正規 0.2361513 0.3327151 0.710 0.478
JOB自営 -0.3075289 0.4399220 -0.699 0.485
JOB無職・学生 -0.3382828 0.3684626 -0.918 0.359
fincome -0.0023254 0.0003127 -7.438 1.03e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
上|中の上 -6.74054 0.70436 -9.570
中の上|中の下 -2.28532 0.30944 -7.385
中の下|下の上 -0.03703 0.28139 -0.132
下の上|下の下 2.33333 0.42776 5.455標準正規分布に近似してz検定を行っている。
僅かに数値は異なるが9，ほぼ同様の結果である。
clm( )の方が有意確率も添えて出力する。
発展1 ログリニアモデル(log-linear model; 対数線型モデル)
本文で紹介しているログリニア・モデル（対数線型モデル）を実行してみる。
最初に，標準のstats::loglin( ) で実演していくが，実際には末尾に示す MASS::loglm( ) を使った方が便利である。
対象となる2元分割表の作成
t01 <- with(d01, table(Marital, SEX))addmargins(t01) SEX
Marital 男性 女性 Sum
unmarried 47 32 79
married 111 161 272
separated 5 24 29
Sum 163 217 380round(prop.table(t01, margin = 2)*100, 1) SEX
Marital 男性 女性
unmarried 28.8 14.7
married 68.1 74.2
separated 3.1 11.1logF <- loglin(t01, margin = list(c(1,2)), fit = T, param = T)2 iterations: deviation 0 logM <- loglin(t01, margin = list(1,2), fit = T, param = T)2 iterations: deviation 0 3元分割表
次に，性，学歴，婚姻上の地位の3つの変数についてログリニアで分析してみる。
学歴変数のfactor化
d01$School <- factor(d01$edu2, levels = 1:3, labels = c("高校", "短大・高専", "四大以上")) この飽和モデルは必然的にデータに一致するものであり，これ自体は意味はない。
ここからどれだけモデルを簡潔にしても適合度が有意に低下しないかがログリニア・モデルの基本的考え方である。
まずは，3次の交互作用を除外してみる。取り敢えずモデルの逸脱度に関心があるので，param = Tとfit = Tは省略する。
# 2次交互作用を全て含んだモデル
log12.23.13 <- loglin(t3, margin = list(c(1, 2), c(2, 3), c(1, 3)))4 iterations: deviation 0.01693409 log12.23.13$lrt
[1] 4.941315
$pearson
[1] 6.262066
$df
[1] 4
$margin
$margin[[1]]
[1] "Marital" "School"
$margin[[2]]
[1] "School" "SEX"
$margin[[3]]
[1] "Marital" "SEX" pchisq(log12.23.13$lrt, log12.23.13$df, lower.tail = F)[1] 0.2933721pchisq(log12.23.13$pearson, log12.23.13$df, lower.tail = F)[1] 0.1804132逸脱度もカイ二乗統計量も有意にはならない。つまり，飽和モデルよりもズレの指標は多少悪化するが，その悪化の程度は統計的に有意ではない。飽和モデルよりも簡潔なこのモデルを採用しても，データには誤差の範囲で同じ程度に適合していると考える事が出来る。よって，必要以上に複雑な飽和モデルよりもこのモデルの方が好ましいとみなす。
ここから更にモデルを簡略化出来るかどうか検討する。
分析結果のうち，逸脱度，カイ二乗統計量と自由度のみ表示させ，それらの有意確率を計算する。
# 2次交互作用を1つ除外したモデル
log12.23 <- loglin(t3, margin = list(c(1, 2), c(2, 3)))2 iterations: deviation 8.881784e-16 log12.23$margin[[1]]
[1] "Marital" "School"
[[2]]
[1] "School" "SEX" sprintf("逸脱度%.2f, カイ二乗値%.2f，自由度%2d", log12.23$lrt, log12.23$pearson, log12.23$df)[1] "逸脱度20.06, カイ二乗値19.63，自由度 6"sprintf("尤度比検定 p=%.3f，カイ二乗検定 p=%.3f",
pchisq(log12.23$lrt, log12.23$df, lower.tail = F),
pchisq(log12.23$pearson, log12.23$df, lower.tail = F))[1] "尤度比検定 p=0.003，カイ二乗検定 p=0.003"1%水準で有意になる。つまり，逸脱度の悪化は偶然の範囲とは言えないので，このモデルは採用できない。
log12.13 <- loglin(t3, margin = list(c(1, 2), c(1, 3)))2 iterations: deviation 8.881784e-16 log12.13$margin[[1]]
[1] "Marital" "School"
[[2]]
[1] "Marital" "SEX" sprintf("逸脱度%.2f, カイ二乗値%.2f，自由度%2d",
log12.13$lrt,
log12.13$pearson,
log12.13$df)[1] "逸脱度35.29, カイ二乗値34.02，自由度 6"sprintf("尤度比検定 p=%.3f，カイ二乗検定 p=%.3f",
pchisq(log12.13$lrt, log12.13$df, lower.tail = F),
pchisq(log12.13$pearson, log12.13$df, lower.tail = F))[1] "尤度比検定 p=0.000，カイ二乗検定 p=0.000"これも高度に有意になる。
log23.13 <- loglin(t3, margin = list(c(2, 3), c(1, 3)))2 iterations: deviation 7.105427e-15 log23.13$margin[[1]]
[1] "School" "SEX"
[[2]]
[1] "Marital" "SEX" sprintf("逸脱度%.2f, カイ二乗値%.2f，自由度%2d",
log23.13$lrt,
log23.13$pearson,
log23.13$df)[1] "逸脱度11.85, カイ二乗値12.02，自由度 8"sprintf("尤度比検定 p=%.3f，カイ二乗検定 p=%.3f",
pchisq(log23.13$lrt, log23.13$df, lower.tail = F),
pchisq(log23.13$pearson, log23.13$df, lower.tail = F))[1] "尤度比検定 p=0.158，カイ二乗検定 p=0.150"このモデルは，飽和モデルからの逸脱度の悪化が有意ではないので，同等のズレでより簡潔な，好ましいモデルであると言える。
2次交互作用を全て含んだモデルと比べても，逸脱度は6.9083105悪化しているが，自由度は4増えており，有意確率は0.1408137で有意ではない。
除外したのは婚姻上の地位と学歴の関連で，学歴によって婚姻上の地位は変わらない，と云う事を意味する。
では，このモデルを基準として，更に2次の交互作用項を除外する事が出来るかどうか検討してみよう。
log23 <- loglin(t3, margin = list(c(2, 3)))2 iterations: deviation 7.105427e-15 log23$margin[[1]]
[1] "School" "SEX" sprintf("Δ逸脱度%.2f，Δ自由度%2d",
log23$lrt - log23.13$lrt,
log23$df - log23.13$df)[1] "Δ逸脱度273.42，Δ自由度 4"sprintf("尤度比検定 p=%.3f",
pchisq(log23$lrt - log23.13$lrt,
log23$df - log23.13$df,
lower.tail = F))[1] "尤度比検定 p=0.000"log13 <- loglin(t3, margin = list(c(1, 3)))2 iterations: deviation 0 log13$margin[[1]]
[1] "Marital" "SEX" sprintf("Δ逸脱度%.2f，Δ自由度%2d",
log13$lrt - log23.13$lrt,
log13$df - log23.13$df)[1] "Δ逸脱度31.34，Δ自由度 4"sprintf("尤度比検定 p=%.3f",
pchisq(log13$lrt - log23.13$lrt,
log13$df - log23.13$df,
lower.tail = F))[1] "尤度比検定 p=0.000" 性別を含んだ交互作用項を除外すると，いずれにしても有意に逸脱度（適合度）が悪化するので，これ以上モデルを単純化する事は出来ない。
“MASS”のloglm( )関数
次は，MASSパッケイジ中のloglm( )関数でログリニアモデルの分析をしてみよう。まずは上で最終的に採用されたモデルを分析してみる。
loglm2 <- loglm(~ Marital * SEX + School * SEX, data = t3)
loglm2Call:
loglm(formula = ~Marital * SEX + School * SEX, data = t3)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 11.84963 8 0.1580404
Pearson 12.01917 8 0.1503508extractAIC(loglm2)[1] 10.00000 31.84963loglm2$param # パラメタ推定値を表示$`(Intercept)`
[1] 2.486074
$Marital
unmarried married separated
0.004035079 1.255662974 -1.259698053
$School
高校 短大・高専 四大以上
0.07419316 -0.11532277 0.04112961
$SEX
男性 女性
-0.2788494 0.2788494
$Marital.SEX
SEX
Marital 男性 女性
unmarried 0.4304338 -0.4304338
married 0.0726238 -0.0726238
separated -0.5030576 0.5030576
$School.SEX
SEX
School 男性 女性
高校 0.04095381 -0.04095381
短大・高専 -0.39368457 0.39368457
四大以上 0.35273076 -0.35273076尤度比やピアソンカイ二乗値，自由度は，先のloglin( )関数の結果と一致しているが，有意確率やAICまで自動で表示してくれるので便利である。
パラメタ推定値もloglin( )関数の結果と同じである。
ここまではわざわざMASSパッケイジのloglm( )を使う理由はない。
しかし，loglin( )関数の結果のオブジェクトでは，変数選択の為のstep( )関数もstepAIC( )関数も使用できないが，MASSのloglm( )関数ではこれらが使えるので，変数を削除してモデルを簡潔に出来るかどうかを検討するにはとても便利である。
最初に飽和モデルを設定してstep( )にかけるところからやってみよう。
loglmF <- loglm(~ Marital * SEX * School, data = t3) # 飽和モデル
step(loglmF)Start: AIC=36
~Marital * SEX * School
Df AIC
- Marital:SEX:School 4 32.941
<none> 36.000
Step: AIC=32.94
~Marital + SEX + School + Marital:SEX + Marital:School + SEX:School
Df AIC
- Marital:School 4 31.850
<none> 32.941
- Marital:SEX 2 44.059
- SEX:School 2 59.290
Step: AIC=31.85
~Marital + SEX + School + Marital:SEX + SEX:School
Df AIC
<none> 31.850
- Marital:SEX 2 43.812
- SEX:School 2 59.043Call:
loglm(formula = ~Marital + SEX + School + Marital:SEX + SEX:School,
data = t3, evaluate = FALSE)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 11.84963 8 0.1580404
Pearson 12.01917 8 0.1503508 この結果を見ると，最初に3次の交互作用が除外され，次に婚姻上の地位と学歴の2次の交互作用が除外されている。それ以上は削除される交互作用項が存在しない事が最後に示され，最終的に採用されたモデルの尤度比検定結果が表示されている。AICなどの値は本文で表に纏めたものと一致している事を確認して欲しい。パラメタ推定値は既にloglm2の結果を表示してあるので省略する。
表記法もlm( )やglm( )などと共通性が高く，MASS::loglm( ) を使う方が便利であると言えよう。
Footnotes
実際にその様に推測される事は余りないが。↩︎
個人年収は関連があるとしてもその因果関係の解釈はやや難しいだろう。ここでは示さないが，性別と個人収入の交互作用項を入れると，個人年収が未婚確率を下げるのは男性だけで，女性には関係ない事が分かる。↩︎
殆どの場合は丁度.5になる事は無いので，.5以上としても.5を超えるとしても大抵結果は変わらない。↩︎
自由度の1は切片（平均値）に相当する。↩︎
多元配置の分散分析で使用したものである。↩︎
ダブルコロン(::)の前がパッケイジ名，その後が関数名である。↩︎
この点，一般に疑似相関と呼ばれるものとは異なる。疑似相関は，「相関と似ているが非なるもの」ではなく，紛れもなく相関である。単にそれが因果関係と解釈できないと云うだけである。ここから，疑似相関ではなく表層的相関と呼ぶのが正しいともいえる。英語ではそもそもpseudo correlationではなくてspurious correlationである。↩︎
上が1で下の下が5であるから，係数の符号がマイナスだと階層が高まる。↩︎
誤差の範囲と言えるだろうか。↩︎