多項ロジスティック回帰分析 – 2 – Google Geminiに聞きながらnnetで解析

Google Geminiを使って、多項ロジスティック回帰分析をさらに探索してみましょう。

Rのパッケージnnetを用いて、医療のデータではありませんがサンプルデータで実際に分析を実行し、説明変数の値から推定するところまで、Gemini 2.5 Flashに聞きながら、やってみます。

まず、試しに使えるデータセットが無いか聞いてみます。

Rパッケージのnnetに付属しているデータセットは？

付属するデータセットはないという回答でした。 その代わりに R に標準で付属しているデータセットがいくつか リストアップされました。 irisは、アヤメのデータセットで、分類問題の例としてよく使われますとのことです。

irisのデータセットを変数datに格納するRのスクリプトは？

dat <- iris

そこで次のスクリプトを実行し、どのようなデータか確認しました。コンソールへの出力です。

ここでは、サンプル1～4までのデータを示します。実際にはアヤメ150個体のデータです。連続変数の説明変数が4種類、従属変数はSpeciesで、アヤメの種類の様です。アヤメの種類はいくつあるのか調べてみましょう。

データセットdatのSpeciesの属性はいくつあるか調べるには？

Speciesは文字列の名義変数で、nlevels(dat$Species)で何種類含まれているかを知ることができるようです。

３種類であることがわかりました。

名義変数のそれぞれの値はlevelsと表現するようなので、それぞれの名義変数の値を知るにはどうするかを聞いてみます。

Speciesのレベルを知るには？

levels(dat$Species)でわかるようです。

setosa, versicolor, virginicaの3種類であることが確認できました。データフレームdatのカラム名を出力させたいと思います。

datのカラム名を出力するには？

colnames(dat)で、あるいはnames(dat)でカラム名を出力できるそうです。

さて、以上からこのデータセットirisは”Sepal.Length”、”Sepal.Width”、 “Petal.Length” “Petal.Width”の4つの説明変数から”Species”を推定するために、多項ロジスティック回帰分析を適用できそうなことが分かりました。nnetを使って多項ロジスティック回帰分析を実行してみようと思います。

Rのパケージnnetを使って、説明変数”Sepal.Length”、”Sepal.Width”、 “Petal.Length” “Petal.Width”から従属変数”Species”を推定するための多項ロジスティック回帰分析を実行するスクリプトは？変数datにirisのデータが格納されています。

nnetパッケージのインストール、パッケージの読み込み、データの準備の解説が提示されます。

# nnetパッケージを読み込む 

library(nnet)

これらに続いて、多項ロジスティック回帰分析の実行の解説が提示され、multinom関数を使って分析するスクリプトが提示されます。その部分を以下に示します。

# 多項ロジスティック回帰モデルの構築

# trace = FALSE は、最適化の過程を表示しないようにするためのものです

model <- multinom(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, data = dat, trace = FALSE)

# または、すべての説明変数を使用する場合

# model <- multinom(Species ~ ., data = dat, trace = FALSE)

そして、結果の確認のスクリプトは：

# モデルの要約を表示 summary(model)

これらをRのエディタにコピー・貼り付けして実行してみます。

このような結果が出力されました。実は、Geminiが出力例として提示した内容とは異なっていました。この結果を見ると、setosaを参照基準として、versicolor、virginicaに対する係数が計算されています。Geminiの出力例では、参照基準がvirginicaになっており、setosa、versicolorに対する係数が提示されていました。実際のデータがこちらです。

それでは、本当にそうなのか、multinom関数でreferenceをvirginicaに変えるにはどうしたらいいかGeminiに聞いてみましょう。

multinom関数でレファレンスをvirginicaにするにはどうしますか？

デフォルトでは、従属変数の並び順で最初の変数を参照基準referenceにするという説明がされ、もしvirginiaに変更するのであれば、次のスクリプトを実行するという説明でした。従って、上記の出力例はデフォルトでの解析結果ではないということになります。

参照基準の変更ができるかやってみます。

dat$Species <- relevel(dat$Species, ref = “virginica”)

# レベルが正しく変更されたか確認 

levels(dat$Species)

さて、今度は参照基準virginiaに対する、setosa、versioclorに対する係数の値が出力されました。しかし、これも上記の出力例とは違う値になっています。つまり、Geminiが出力例（一部）として提示したデータの由来は分かりませんが、どうも間違いのようです。

さて、元のデータのまま、再度同じ解析を実行しておきます。そして、Geminiからの回答にあった、係数と標準誤差、p値（両側検定）の値を出力してみます。

# 係数と標準誤差を取得 

z <- summary(model)$coefficients / summary(model)$standard.errors 

# p値を計算 (両側検定) 

p_values <- (1 – pnorm(abs(z), 0, 1)) * 2 print(p_values)

続いて、説明変数の値から、Speciesの推定predictionをするにはどうするかGeminiに聞いてみましょう。まず、分析の元のデータ、つまり訓練データに対して適用してみます。

datの4つの説明変数の値からSpeciesを推定するにはどうしますか？

predict関数を使用することと、訓練データの場合の推定結果、正解率の計算、推定確率の計算のためのスクリプトが提示されました。Species、つまりアヤメの名称をクラスと呼んでいます。

# 訓練データに対する推定されたクラス

predicted_class <- predict(model, newdata = dat, type = “class”)

print(head(predicted_class))

# 実際のSpeciesと比較して、どの程度合っているか確認する（混同行列）

confusion_matrix <- table(Actual = dat$Species, Predicted = predicted_class)

print(confusion_matrix)

# 正解率の計算

accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)

cat(“モデルの正解率:”, round(accuracy * 100, 2), “%\n”)

# 訓練データに対する各クラスの推定確率

predicted_probs <- predict(model, newdata = dat, type = “probs”)

print(head(predicted_probs))

irisのデータは、nrow(dat)で150、つまり150個体のアイリスのデータを含んでいます。ここでの出力はhead()関数で最初の6個だけになっています。confusion_matrixつまり混同行列の結果をみると、setosa50個体のはすべてsetosaとして推定され、versicolorは50個体の内49個体が正しくversicolorとして推定され、virginica1は50個体の内49個体が正しくvirginicaとして分類されていることが分かります。正解率の計算では98.67％でした。非常に高い精度です。

最後の、訓練データに対する各クラスの推定確率を見ると個体1はsetosaの確率が1.0でversicolor、verginicaの確率はほぼ0ということが分かります。以降個体6までの結果が出力されていますが、print(predicted_probs)で、150個体すべての結果が出力されます。

さて、Geminiの回答には、新しいデータポイントについてSpeciesを推定する方法、スクリプトについても説明があります。以下がその部分です。臨床の場合であれば、個々の対象者のデータから、診断、あるいは病態を推定することに相当します。ここでは、３つの個体のデータを作成して、推定を行っています。

——–
2. 新しいデータポイントについてSpeciesを推定する

これは、モデルを未知のデータに適用する場合の一般的な使用方法です。新しいデータは、モデル構築時に使用した説明変数と同じカラム名とデータ型を持つデータフレームである必要があります。

例として、新しいアヤメの測定値を使ってSpeciesを推定してみましょう。

R

# 新しいデータポイントを作成

# カラム名はモデル構築時に使用した説明変数と完全に一致させる必要があります

new_iris_data <- data.frame(

  Sepal.Length = c(5.1, 6.0, 7.0),

  Sepal.Width = c(3.5, 2.7, 3.2),

  Petal.Length = c(1.4, 4.5, 6.0),

  Petal.Width = c(0.2, 1.5, 2.0)

)

print(new_iris_data)

# 新しいデータポイントに対する推定されたクラス

predicted_new_class <- predict(model, newdata = new_iris_data, type = “class”)

cat(“\n新しいデータに対する推定されたSpecies:\n”)

print(predicted_new_class)

# 新しいデータポイントに対する各Speciesカテゴリの推定確率

predicted_new_probs <- predict(model, newdata = new_iris_data, type = “probs”)

cat(“\n新しいデータに対する各Speciesの推定確率:\n”)

print(predicted_new_probs)

このように、predict()関数を使うことで、学習済みモデルを用いて説明変数の値から従属変数Speciesを推定することができます。

—————

これらのスクリプトをRで実行してみます。

従って、個体1はsetosa、個体2はversicolor、個体3はvirginicaと推定されています。

さてirisのデータは説明変数が連続変数のみでした。説明変数に名義変数がある場合、どのように扱うのか、特に参照基準の設定について聞いてみましょう。

説明変数に名義変数がある場合、例えば、x1: A, B, C, D；x2: X, Y, Zのような場合、multinom関数を使う場合に、参照基準の設定は？

Rのmultinom関数で、説明変数に名義変数（カテゴリカル変数）がある場合、その説明変数の参照基準（レファレンスレベル）を設定する方法は、従属変数の場合と全く同じ考え方です。

Rの統計モデリング関数（lm, glm, multinomなど）は、因子（factor）型変数を自動的にダミー変数に変換して回帰分析に組み込みます。この際、デフォルトではその因子の最初のレベルが参照基準となり、そのレベルに対応するダミー変数はモデルから除外（係数が0として扱われる）されます。他のレベルの係数は、この参照基準レベルと比較してどのくらい異なるかを示します。

したがって、説明変数x1やx2の参照基準を設定したい場合は、それぞれの変数をfactor型であることを確認し、relevel()関数を使って参照基準にしたいレベルを最初のレベルに設定します。

以下に具体的なスクリプトを示します。

1. サンプルデータの作成

まず、ご提示のx1 (A, B, C, D) と x2 (X, Y, Z) のような名義変数を含むダミーのデータセットを作成します。

R

# データをシミュレート

set.seed(123)

n <- 100 # サンプルサイズ

# 説明変数 x1 (A, B, C, D)

x1_levels <- c(“A”, “B”, “C”, “D”)

x1 <- sample(x1_levels, n, replace = TRUE, prob = c(0.25, 0.25, 0.25, 0.25))

# 説明変数 x2 (X, Y, Z)

x2_levels <- c(“X”, “Y”, “Z”)

x2 <- sample(x2_levels, n, replace = TRUE, prob = c(0.33, 0.33, 0.34))

# 従属変数 Species (irisデータと同様に3カテゴリを想定)

species_levels <- c(“setosa”, “versicolor”, “virginica”)

species <- sample(species_levels, n, replace = TRUE) # ランダムなSpeciesを生成

# データフレームを作成

my_dat <- data.frame(

  Species = as.factor(species), # 従属変数はfactor型

  x1 = as.factor(x1),            # 説明変数x1はfactor型

  x2 = as.factor(x2)             # 説明変数x2はfactor型

)

# データの構造を確認

str(my_dat)

2. 説明変数の参照基準の設定

relevel()関数を使って、各説明変数の参照基準を変更します。

x1の参照基準をBにする場合:
R

my_dat$x1 <- relevel(my_dat$x1, ref = “B”)

levels(my_dat$x1) # 確認: “B” が最初になっているはず

x2の参照基準をZにする場合:
R

my_dat$x2 <- relevel(my_dat$x2, ref = “Z”)

levels(my_dat$x2) # 確認: “Z” が最初になっているはず

3. multinom関数の実行

参照基準を設定したデータフレームmy_datを使ってmultinom関数を実行します。

R

library(nnet)

# (オプション) 従属変数Speciesの参照基準も設定する例 (もし必要なら)

# my_dat$Species <- relevel(my_dat$Species, ref = “virginica”)

# 多項ロジスティック回帰モデルの構築

# trace = FALSE は、最適化の過程を表示しないようにするためのものです

model_nominal <- multinom(Species ~ x1 + x2, data = my_dat, trace = FALSE)

# モデルの要約を表示

summary(model_nominal)

出力の解釈

summary(model_nominal)のCoefficientsセクションを見ると、x1とx2の係数がどのように表示されるかを確認できます。

• x1について:
  x1Bのような係数は表示されず、x1A、x1C、x1Dのような係数が表示されます。これは、それぞれのレベル（A, C, D）が参照基準Bと比較してどうであるかを示します。
  
• x2について:
  x2Zのような係数は表示されず、x2X、x2Yのような係数が表示されます。これは、それぞれのレベル（X, Y）が参照基準Zと比較してどうであるかを示します。
  

重要な注意点:

• factor型であることの確認: 説明変数がfactor型であることを常に確認してください。文字列型のままではRは自動的にダミー変数化しません。
• デフォルトの動作: relevel()を使用しない場合、Rはアルファベット順（または数値順）で最初のレベルを参照基準とします。
• 他の参照基準設定方法: C()関数とcontr.treatments()を組み合わせて、モデルの式内で一時的に参照基準を設定することも可能ですが、relevel()でデータフレーム自体を変更する方が、コードが読みやすくなることが多いです。 例: multinom(Species ~ C(x1, contr.treatment(nlevels(x1), base=2)) + x2, data = my_dat) (この例ではx1の2番目のレベルを基準にする)

このように、説明変数に名義変数がある場合でも、relevel()関数を使用して参照基準を明示的に設定することができます。
—————

この解説で分かるように、説明変数が名義変数の場合は、参照基準はrlevelで設定できること、結果で得られる係数はそれぞれの値に対する値であることが分かります。

さらに、名義変数の場合、Multi-dichotomous variableにあらかじめ手作業で変換する必要があるか聞いたところ、その必要はないという回答でした。

説明変数が名義変数の場合、multi-dichootmous variableにあらかじめ変換する必要はありますか？

いいえ、Rのmultinom関数（および多くの統計モデリング関数、例：lm, glm）を使用する場合、説明変数である名義変数（カテゴリカル変数）を事前に手動でmulti-dichotomous variable（多値二項変数、つまりダミー変数）に変換する必要はありません。

Rは、データフレーム内の変数がfactor型として定義されていれば、multinom関数に渡されたときに、内部的に自動でダミー変数に変換してくれます。これがRのモデリング関数の便利な点の1つです。

Rが自動的に行うこと:

1. factor型変数の検出: モデル式（例: Species ~ x1 + x2）の中にfactor型の変数があると、Rはそれを自動的に認識します。
2. ダミー変数の生成: 例えば、x1がA, B, C, Dの4つのレベルを持つfactor型の場合、Rは自動的に3つのダミー変数を生成します。
   • x1B (1 if x1 is B, 0 otherwise)
   • x1C (1 if x1 is C, 0 otherwise)
   • x1D (1 if x1 is D, 0 otherwise) （これはデフォルトの設定で、最初のレベルAが参照基準となり、そのダミー変数は生成されません。）
3. モデルへの組み込み: これらの生成されたダミー変数が、回帰モデルの独立した予測変数として組み込まれます。

参照基準の設定の重要性:

Rは自動的にダミー変数を生成してくれますが、どのカテゴリを参照基準（ベースライン）とするかは、結果の解釈に影響します。デフォルトでは、factor型の最初のレベルが参照基準になります。

もし、特定のカテゴリを参照基準にしたい場合は、前回の回答で説明したようにrelevel()関数を使ってそのカテゴリを因子の最初のレベルに設定する必要があります。

例（再掲）：

R

# データをシミュレート

my_dat <- data.frame(

  Species = as.factor(sample(c(“setosa”, “versicolor”, “virginica”), 100, replace = TRUE)),

  x1 = as.factor(sample(c(“A”, “B”, “C”, “D”), 100, replace = TRUE)),

  x2 = as.factor(sample(c(“X”, “Y”, “Z”), 100, replace = TRUE))

)

# x1の参照基準を ‘B’ に設定

my_dat$x1 <- relevel(my_dat$x1, ref = “B”)

# x2の参照基準を ‘Z’ に設定

my_dat$x2 <- relevel(my_dat$x2, ref = “Z”)

library(nnet)

model_nominal <- multinom(Species ~ x1 + x2, data = my_dat, trace = FALSE)

summary(model_nominal)

このコードでは、x1とx2がfactor型として定義されていれば、Rが自動的に内部で必要なダミー変数変換を行い、multinom関数はそれらを使ってモデルを構築します。あなたは手動でx1_A, x1_C, x1_Dなどの新しいカラムを作成する必要はありません。

結論:

名義変数をfactor型として適切に定義し、必要であればrelevel()で参照基準を設定するだけで、Rのmultinom関数は適切に処理してくれます。手動でのダミー変数への変換は不要です。
——————

さてこれで、nnetを使った多項ロジスティック回帰分析のやり方は、分かったと思いますがいかがでしたか。Geminiに分からないことを聞きながら、そのまま実行可能なスクリプトが提示されるので、素晴らしいと思います。