Rの関数:multinom {nnet}

Rの関数から multinom {nnet} を確認します。

`multinom`関数とは

multinom関数は、多項ロジスティック回帰モデルの当てはめに用いられます。このモデルは、目的変数が3つ以上のカテゴリを持つ名義尺度である場合に適しています。例えば、選挙における支持政党の選択、製品のブランド選択、あるいは疾患のタイプの分類など、多岐にわたる事象の予測に活用されます。

multinom関数の内部的な動作は、nnetパッケージの基幹関数であるnnet.defaultを、隠れ層のないニューラルネットワーク（size = 0）、そして出力層にソフトマックス関数（softmax = TRUE）を指定して呼び出すことにより、多項ロジスティック回帰と等価な計算を実行しています。

ソースコード内のnnet.default呼び出し

multinom関数の内部では、目的変数Yがどのような形式か（行列か、因子か）によって処理が分岐しますが、最終的にはどちらの分岐でもnnet.defaultという関数を呼び出しています。

分岐 A: 目的変数が因子（factor）の場合

else {
    if (length(offset) <= 1L) {
        mask <- c(FALSE, rep(TRUE, r))
        # ↓↓↓ ここで呼び出しています ↓↓↓
        fit <- nnet.default(X, Y, w, mask = mask, size = 0, 
            skip = TRUE, entropy = TRUE, rang = 0, ...)
    }
# ... (offsetがある場合も同様)

分岐 B: 目的変数が行列（matrix）の場合

if (is.matrix(Y)) {
# ... (中略) ...
    else {
        mask <- c(rep(FALSE, r + 1L), rep(c(FALSE, rep(TRUE, 
            r)), p - 1L))
        # ↓↓↓ ここで呼び出しています ↓↓↓
        fit <- nnet.default(X, Y, w, mask = mask, size = 0, 
            skip = TRUE, softmax = TRUE, censored = censored, 
            rang = 0, ...)
    }

キーとなる3つの引数

size = 0:
これはニューラルネットワークの隠れ層（中間層）のノード（ユニット）の数を0に設定することを意味します。隠れ層を持たないニューラルネットワークは、入力層から直接出力層へ信号が伝達される構造になります。これは、本質的に一般化線形モデル（Generalized Linear Model, GLM）と等価な構造です。
softmax = TRUE (または entropy = TRUE):
- softmax = TRUE: これは、出力層の活性化関数としてソフトマックス関数を使用することを指定します。ソフトマックス関数は、複数の出力の合計が1になるように正規化する関数であり、多クラス分類問題において各クラスに属する「確率」を計算するために用いられます。これは、まさしく多項ロジスティック回帰の確率を計算する式そのものです。
- entropy = TRUE: 目的変数が因子（クラスラベル）の場合、この引数がTRUEに設定されます。これは、モデルの当てはめに使用する損失関数として、交差エントロピー（cross-entropy）を使用することを意味します。多クラス分類において、ソフトマックス関数を出力層に持ち、交差エントロピーを損失関数として学習するモデルは、多項ロジスティック回帰と等価です。

`multinom`関数の引数

library(nnet)
args(multinom)

function (formula, data, weights, subset, na.action, contrasts = NULL, 
    Hess = FALSE, summ = 0, censored = FALSE, model = FALSE, 
    ...) 
NULL

formula:
モデルの構造を定義する式です。目的変数 ~ 説明変数1 + 説明変数2 のように記述します。
data:
モデルで使用するデータが含まれるデータフレームを指定します。
weights:
各観測値に対する重みを指定するベクトルです。デフォルトでは全ての観測値は等しく扱われます。
subset:
データフレームの一部を用いてモデルを構築する場合に、使用する観測値のインデックスを指定します。
na.action:
データに欠損値（NA）が含まれている場合の処理方法を指定します。デフォルトは na.omit で、欠損値を含む行を無視します。
contrasts:
モデルに含まれる因子（カテゴリカル変数）の対比（コントラスト）を指定します。
Hess:
論理値（TRUEまたはFALSE）で、モデルのヘッセ行列（パラメータの共分散行列の推定に用いられる）を計算し、結果に含めるかどうかを指定します。TRUEにすると、パラメータの標準誤差が計算されますが、計算時間は増加します。
summ:
データを要約するための内部的な手法を指定する整数値（0, 1, 2, 3）です。通常、ユーザーがこの値を変更する必要はありません。
censored:
目的変数が右側打ち切り（right-censored）データである場合、TRUEに設定します。これは生存時間解析のような特定の文脈で用いられます。
model:
論理値で、TRUEに設定すると、返されるオブジェクトにモデルフレームが含まれます。
...:
nnet.default関数に渡されるその他の引数です。例えば、最適化の最大反復回数を指定する maxit などがあります。

シミュレーション

multinom関数を確認するため、シミュレーションを実行します。ここでは、目的変数をZ（カテゴリZ1, Z2, Z3を持つ）、説明変数をA, B, C, Dとするサンプルデータを作成し、分析を行います。

なお、有意水準は5%とします。

サンプルデータの作成

始めに、分析の対象となるサンプルデータを作成します。データには、目的変数Zと4つの説明変数A, B, C, Dが含まれます。

# 再現性を確保するために乱数の種を設定
seed <- 20251118
set.seed(seed)

# サンプルサイズを設定
n <- 500

# 4つの説明変数を作成（標準正規分布に従う乱数）
A <- rnorm(n)
B <- rnorm(n)
C <- rnorm(n)
D <- rnorm(n)

# ロジットを計算するための係数を設定
# 全ての説明変数が有意な結果となるよう、比較的大きな値を設定
# 基準カテゴリは「Z1」とする
# Z2 vs Z1
b_z2 <- c(0.5, 1.5, -1.0, 0.8, -1.2) # (Intercept), A, B, C, D
# Z3 vs Z1
b_z3 <- c(-0.5, -0.9, 1.2, 1.5, 0.7) # (Intercept), A, B, C, D

# デザイン行列を作成（切片項を含む）
X <- cbind(1, A, B, C, D)

# 各選択肢の線形予測子（対数オッズ）を計算
logit_z2 <- X %*% b_z2
logit_z3 <- X %*% b_z3

# 各選択肢の確率を計算（ソフトマックス関数）
# 基準カテゴリ（Z1）の対数オッズは0とする
prob_z1 <- 1 / (1 + exp(logit_z2) + exp(logit_z3))
prob_z2 <- exp(logit_z2) / (1 + exp(logit_z2) + exp(logit_z3))
prob_z3 <- exp(logit_z3) / (1 + exp(logit_z2) + exp(logit_z3))

# 各確率を結合して行列にする
probs <- cbind(prob_z1, prob_z2, prob_z3)

# 各観測値のカテゴリを確率に基づいて決定
z_numeric <- apply(probs, 1, function(p) sample(1:3, 1, prob = p))
Z <- factor(z_numeric,
  levels = c(1, 2, 3),
  labels = c("Z1", "Z2", "Z3")
)

# データフレームを作成
sample_data <- data.frame(Z, A, B, C, D)

# 作成したデータの一部を確認
cat("--- 作成されたサンプルデータの一部を確認 ---\n")
print(str(sample_data))

# 各カテゴリの観測数を確認
cat("\n--- 各カテゴリの観測数 ---\n")
print(table(sample_data$Z))

--- 作成されたサンプルデータの一部を確認 ---
'data.frame':   500 obs. of  5 variables:
 $ Z: Factor w/ 3 levels "Z1","Z2","Z3": 2 3 2 2 1 2 3 3 2 2 ...
 $ A: num  -0.416 -1.759 -0.262 0.474 2.3 ...
 $ B: num  -0.866 0.646 -0.977 -1.067 0.921 ...
 $ C: num  0.1602 -0.0855 -0.8467 -1.1625 -0.6313 ...
 $ D: num  -0.1266 0.0343 0.123 -0.0577 0.3391 ...
NULL

--- 各カテゴリの観測数 ---

 Z1  Z2  Z3 
119 237 144

multinom関数によるモデル構築と結果の解釈

作成したサンプルデータを用いて、multinom関数で多項ロジットモデルを構築します。目的変数をZ、説明変数をA,B,C,Dとします。

# 多項ロジスティック回帰モデルを構築
# Z ~ . は、Z を目的変数とし、データフレーム内の他の全ての変数を説明変数とすることを意味します
model <- multinom(Z ~ ., data = sample_data, Hess = TRUE)

# モデルの要約情報を表示
cat("\n\n--- モデルの要約 ---\n")
summary_model <- summary(model)
print(summary_model)

# weights:  18 (10 variable)
initial  value 549.306144 
iter  10 value 343.095556
iter  20 value 295.763137
final  value 295.761435 
converged


--- モデルの要約 ---
Call:
multinom(formula = Z ~ ., data = sample_data, Hess = TRUE)

Coefficients:
   (Intercept)          A         B         C          D
Z2   0.5643645  1.4679517 -1.214588 0.9501563 -1.4304888
Z3  -0.5013325 -0.8378469  1.093368 1.6582902  0.4564617

Std. Errors:
   (Intercept)         A         B         C         D
Z2   0.1651730 0.1910905 0.2025989 0.1705136 0.1851975
Z3   0.2248707 0.2013078 0.2166420 0.1997023 0.1820109

Residual Deviance: 591.5229 
AIC: 611.5229

Call

Call:
multinom(formula = Z ~ ., data = sample_data, Hess = TRUE)

この行は、このモデルを生成するために実行されたRのコマンド（関数呼び出し）をそのまま表示しています。Z ~ .という式は、データフレームsample_data内のZを目的変数とし、それ以外の全ての変数（A, B, C, D）を説明変数として使用したことを示します。また、Hess = TRUEが指定されているため、標準誤差の計算に必要なヘッセ行列が算出されています。

Coefficients（係数）

Coefficients:
   (Intercept)          A         B         C          D
Z2   0.5643645  1.4679517 -1.214588 0.9501563 -1.4304888
Z3  -0.5013325 -0.8378469  1.093368 1.6582902  0.4564617

こちらのセクションは、推定された回帰係数を示しています。multinom関数は、目的変数のカテゴリの中から一つを基準（ベースライン）カテゴリとして扱います。今回のケースでは、Z1が基準カテゴリとして自動的に選択されています。

解釈の基本:
各係数は、基準カテゴリ（Z1）と比較した際の、他のカテゴリの選択されやすさを対数オッズ（log-odds）で表します。
行の意味:
- Z2: 「Z1が選ばれる場合と比較してZ2が選ばれる」対数オッズに関する係数です。
- Z3: 「Z1が選ばれる場合と比較してZ3が選ばれる」対数オッズに関する係数です。
具体的な係数の解釈:
- 例えば、Z2の行にあるAの係数 1.468 は、「他の変数が一定の時、説明変数Aが1単位増加すると、Z1よりもZ2が選択される対数オッズが1.468増加する」ことを意味します。対数オッズの増加は、その選択肢がより選ばれやすくなることを示します。
- 同様に、Z2の行にあるBの係数 -1.215 は、「Bが1単位増加すると、Z1よりもZ2が選択される対数オッズが1.215減少する」ことを示します。
オッズ比への変換:
対数オッズは、指数関数 exp() を用いることで、より直感的に理解しやすいオッズ比に変換できます。例えば、AがZ2に与える影響のオッズ比は exp(1.468) ≈ 4.34 となり、「Aが1単位増加すると、Z1に対するZ2の選択オッズが約4.34倍になる」と解釈できます。

Std. Errors（標準誤差）

Std. Errors:
   (Intercept)         A         B         C         D
Z2   0.1651730 0.1910905 0.2025989 0.1705136 0.1851975
Z3   0.2248707 0.2013078 0.2166420 0.1997023 0.1820109

こちらのセクションは、上で推定された各係数の標準誤差を示します。標準誤差は、係数の推定値がどの程度の精度を持つか、すなわち推定値のばらつきの大きさを示す指標です。この値が小さいほど、係数の推定がより正確であることを意味します。標準誤差は、係数の統計的有意性を検定するためのp値や、信頼区間を計算する際に利用されます。

Residual Deviance（残差逸脱度）と AIC（赤池情報量規準）

Residual Deviance: 591.5229 
AIC: 611.5229

これらの指標は、モデル全体の適合度（goodness-of-fit）、すなわちモデルが観測データをどれだけうまく説明できているかを評価するためのものです。

Residual Deviance:
モデルの当てはまりの悪さを示す指標で、逸脱度とも呼ばれます。値が小さいほど、モデルがデータによく適合していることを示唆します。この値は、異なるモデル同士を比較する際の尤度比検定に用いられます。
AIC:
モデルの適合度と、モデルの複雑さ（パラメータの数）の両方を考慮した指標です。モデルが複雑になる（説明変数を増やす）と残差逸脱度は必ず減少しますが、AICは不要な変数を加えることに対してペナルティを課します。それゆえ、複数のモデル候補の中から最適なモデルを選択する際に利用され、AICの値が最も小さいモデルが最良であると判断されます。今回のAIC 611.5229 は、残差逸脱度 591.5229 に、モデルのパラメータ数（10個）の2倍である20を加算して計算されています (591.5229 + 2 * 10 = 611.5229)。

# confint() 関数を multinom オブジェクトに適用
# level = 0.95 はデフォルトですが、明示的に指定しています
confidence_intervals <- confint(model, level = 0.95)
cat("--- confint()関数による95%信頼区間 ---\n")
print(confidence_intervals)

--- confint()関数による95%信頼区間 ---
, , Z2

                 2.5 %     97.5 %
(Intercept)  0.2406314  0.8880977
A            1.0934212  1.8424822
B           -1.6116746 -0.8175016
C            0.6159557  1.2843569
D           -1.7934692 -1.0675085

, , Z3

                  2.5 %      97.5 %
(Intercept) -0.94207089 -0.06059411
A           -1.23240291 -0.44329083
B            0.66875747  1.51797837
C            1.26688088  2.04969948
D            0.09972686  0.81319659

シミュレーションのデータ生成時に設定した真の係数値と、multinomモデルが算出した各係数の95%信頼区間を比較し、真の値が区間内に含まれているかを確認します。

カテゴリ	係数名	真の値	信頼区間下限 (2.5%)	信頼区間上限 (97.5%)	区間内に含まれるか?
Z2	(Intercept)	0.5	0.2406	0.8881	はい
	A	1.5	1.0934	1.8425	はい
	B	-1.0	-1.6117	-0.8175	はい
	C	0.8	0.6160	1.2844	はい
	D	-1.2	-1.7935	-1.0675	はい
Z3	(Intercept)	-0.5	-0.9421	-0.0606	はい
	A	-0.9	-1.2324	-0.4433	はい
	B	1.2	0.6688	1.5180	はい
	C	1.5	1.2669	2.0497	はい
	D	0.7	0.0997	0.8132	はい

# 係数のp値を計算（z検定）
# このp値が設定した有意水準を下回っていれば、その係数が統計的に有意であると判断されます
cat("--- 係数の有意性（p値）---\n")
z_stats <- summary_model$coefficients / summary_model$standard.errors
p_values <- (1 - pnorm(abs(z_stats), 0, 1)) * 2
print(p_values)

--- 係数の有意性（p値）---
    (Intercept)            A            B           C            D
Z2 0.0006335978 1.576517e-14 2.034373e-09 2.51372e-08 1.132427e-14
Z3 0.0257855698 3.154463e-05 4.490624e-07 0.00000e+00 1.214576e-02

全ての切片および係数のp値は設定した有意水準を下回っています。

以上です。