Copilotに聞きながらメタアナリシスを実行：metaパッケージを使って

まずRのメタアナリシス用のパッケージに何があるか聞いてみます。

Copilotの回答は青いブロックの右端のアイコンをクリックすると表示したり隠したりできます。

メタアナリシスのためのRのパッケージを教えて。

Rにはメタアナリシスを行うための便利なパッケージがいくつかあります。代表的なものを紹介しますね。

  • metafor: メタアナリシスのための強力なパッケージで、固定効果モデルやランダム効果モデルの適用、フォレストプロットの作成などが可能です。
  • meta: メタアナリシスのための包括的なパッケージで、異質性の評価や出版バイアスの検討など、多くの機能を備えています。
  • rmeta: よりシンプルなメタアナリシスの実装を提供するパッケージで、基本的な統計解析を行うのに適しています。

これらのパッケージを使えば、メタアナリシスのデータ統合や異質性の評価、視覚化などがスムーズに行えます。どのパッケージを選ぶかは、分析の目的や必要な機能によります。試してみたいものはありますか？

Copilotの提案の中からmetaパッケージを使ってみようと思います。二値変数アウトカムの場合のメタアナリシスを実行するスクリプトを提案してもらいましょう。

Rのmetaパッケージで二値変数アウトカムの場合のメタアナリシスを
実行し、フォレストプロット、ファンネルプロットを作成します。
分散逆数法でランダム効果モデルを用います。対照群の症例数にnc、イベント数にrc、介入群の症例数にnt、イベント数にrtという変数名を使います。プロットは新規のウインドウにそれぞれ出力してください。

meta パッケージを使って、分散逆数法 に基づくランダム効果モデルを適用し、二値変数アウトカムのメタアナリシスを実行するRスクリプトを紹介します。フォレストプロットとファンネルプロットを 新規ウインドウ に出力します。

スクリプト

# 必要なパッケージの読み込み
library(meta)

# 仮のデータセット
dat <- data.frame(
  study = c("Study 1", "Study 2", "Study 3", "Study 4", "Study 5"),
  rt = c(10, 20, 15, 25, 30),  # 介入群のイベント数
  nt = c(100, 200, 150, 250, 300),  # 介入群の総数
  rc = c(5, 15, 10, 20, 25),  # 対照群のイベント数
  nc = c(100, 200, 150, 250, 300)   # 対照群の総数
)

# メタアナリシスの実行（分散逆数法を用いたランダム効果モデル）
meta_analysis <- metabin(rt, nt, rc, nc, data = dat, 
                         studlab = study, sm = "RR", method = "Inverse", 
                         fixed = FALSE, random = TRUE)
                         common = FALSE
# メタアナリシスの結果をコンソールに出力
print(meta_analysis)

# 新規ウインドウでフォレストプロットを作成
windows()
forest(meta_analysis, xlab = "Risk Ratio", ref = 1)

# 新規ウインドウでファンネルプロットを作成
windows()
funnel(meta_analysis)

ポイント

  • metabin() 関数を使用して、二値変数アウトカムのメタアナリシスを実行。
  • sm = "RR" でリスク比（Risk Ratio）を指定。
  • method = "Inverse" で 分散逆数法 を適用。
  • random = TRUE で ランダム効果モデル を使用。
  • windows() を使って 新規ウインドウ にプロットを出力。

このスクリプトを実行すれば、分散逆数法を用いたメタアナリシスを適切に実施できます。試してみて、調整したい点があれば教えてくださいね！詳しい解説は こちら で確認できます。

仮のデータセットの作成から、メタアナリシスの実行のスクリプトがコピーできる形で提案されたのでそのままRで実行したらこのような結果になりました。τ二乗値、I二乗値、統合値に対するZ値、p値も出力されています。

metabin()関数の実行でエラーメッセージが出ています。固定効果モデルの場合fixed = FALSEではなく、common = FALSEと記述するということです。実際にそのように書き換えて実行するとエラーメッセージは出なくなりました。上のスクリプトの赤字の部分です。

フォレストプロットに提示されているp値はQ統計値に対するp値で、統合値に対するp値はコンソールに出力されている、z p-valueの方です。この例では、0.0520でした。

さて、自分のデータセットでメタアナリシスをする場合、仮のデータセットの部分を書き換えて、以下同様にスクリプトを実行するのがひとつの方法です。

Excelでデータを用意する場合、study, nc, rc, nt, rtのラベルでそれぞれのデータを入力し、データ範囲を選択して、コピー操作をして、Rで次のスクリプトを実行します。

dat = read.table("clipboard", sep = "\t", header = TRUE)

例えば、このようなデータを用意して、左上のセルから右下のセルの範囲を選択して、コピー（Ctr+C)し、Rの方で、上記のスクリプトを実行すると変数datにデータが格納され、それをdatで確認できます。それ以下は同じスクリプトを実行します。次のような結果が得られました。

それでは、ファンネルプロットの非対称性のBeggの検定、Eggerの検定はどうするか聞いてみます。

metaパッケージでファンネルプロットの非対称性のBeggの検定、
Eggerの検定のためのスクリプトは？

その部分だけを示します。

# Begg検定（順位相関検定）
begg_test <- metabias(meta_analysis, method = "rank")
print(begg_test)

# Egger検定（回帰分析による検定）
egger_test <- metabias(meta_analysis, method = "linreg")
print(egger_test)

今回の例では、これらの検定は研究数が10個以上でないとできないという、エラーメッセージが出ました。

さて、統合値に対するPrediction Intervalも出力するにはどうするか聞いてみます。

metaパッケージで統合値に対するPrediction Intervalも出力し、
フォレストプロットにも追加するにはどうしますか？

関連したところだけを以下に示します。また、addpredとなっていた引数名をpredictionに訂正してあります。


# メタアナリシスの結果をコンソールに出力（予測区間を含む）
print(meta_analysis, prediction = TRUE)

# フォレストプロットの作成（予測区間を追加）
forest(meta_analysis, xlab = "Risk Ratio", ref = 1, prediction = TRUE)

これで出力されたフォレストプロットがこちらです。

さて、今回もCopilotの提案には間違いがありました。commonとすべきところがfixedに、predictionとすべきところがaddpredになっていました。それについて次の様に聞いてみたところ、正しいスクリプトが提案されました。

addpredではないのでは？

その通りです！meta パッケージのフォレストプロットで 予測区間 (Prediction Interval, PI) を追加するには、prediction = TRUE を指定します。

予測区間を出力し、フォレストプロットに追加する方法

# 必要なパッケージの読み込み
library(meta)

# 仮のデータセット
dat <- data.frame(
  study = c("Study 1", "Study 2", "Study 3", "Study 4", "Study 5"),
  rt = c(10, 20, 15, 25, 30),  # 介入群のイベント数
  nt = c(100, 200, 150, 250, 300),  # 介入群の総数
  rc = c(5, 15, 10, 20, 25),  # 対照群のイベント数
  nc = c(100, 200, 150, 250, 300)   # 対照群の総数
)

# メタアナリシスの実行（リスク比を使用）
meta_analysis <- metabin(rt, nt, rc, nc, data = dat, 
                         studlab = study, sm = "RR", method = "Inverse", 
                         fixed = FALSE, random = TRUE)

# メタアナリシスの結果をコンソールに出力（予測区間を含む）
print(meta_analysis, prediction = TRUE)

# フォレストプロットの作成（予測区間を追加）
forest(meta_analysis, xlab = "Risk Ratio", ref = 1, prediction = TRUE)

ポイント

  • print(meta_analysis, prediction = TRUE) を指定すると、統合値の 予測区間 がコンソールに出力されます。
  • forest(meta_analysis, prediction = TRUE) を指定すると、フォレストプロットに 予測区間 が追加されます。

この方法を試してみて、さらにカスタマイズしたい点があれば教えてくださいね！詳しい解説は こちら で確認できます。

metabin()のfixedも違うのでは？

その通りです！meta パッケージの metabin() 関数では、固定効果モデルとランダム効果モデルの指定方法が fixed ではなく、common（固定効果モデル）と random（ランダム効果モデル）を使います。

正しいメタアナリシスの実行方法

# 必要なパッケージの読み込み
library(meta)

# 仮のデータセット
dat <- data.frame(
  study = c("Study 1", "Study 2", "Study 3", "Study 4", "Study 5"),
  rt = c(10, 20, 15, 25, 30),  # 介入群のイベント数
  nt = c(100, 200, 150, 250, 300),  # 介入群の総数
  rc = c(5, 15, 10, 20, 25),  # 対照群のイベント数
  nc = c(100, 200, 150, 250, 300)   # 対照群の総数
)

# メタアナリシスの実行（リスク比を使用）
meta_analysis <- metabin(rt, nt, rc, nc, data = dat, 
                         studlab = study, sm = "RR", method = "Inverse", 
                         common = FALSE, random = TRUE)

# メタアナリシスの結果をコンソールに出力（予測区間を含む）
print(meta_analysis, prediction = TRUE)

# フォレストプロットの作成（予測区間を追加）
forest(meta_analysis, xlab = "Risk Ratio", ref = 1, prediction = TRUE)

ポイント

  • common = FALSE で 固定効果モデルを無効化。
  • random = TRUE で ランダム効果モデルを適用。
  • print(meta_analysis, prediction = TRUE) で 予測区間を出力。
  • forest(meta_analysis, prediction = TRUE) で フォレストプロットに予測区間を追加。

この方法を試してみて、さらにカスタマイズしたい点があれば教えてくださいね！詳しい解説は こちら で確認できます。

実際にRでスクリプトを実行するとエラーメッセージがでたり、思っているのと違ったことが実行されたりした場合、間違いと思われるところを指摘すると、訂正を提案してくれるようです。

最後に、forest()関数の引数refは何かを聞いてみました。

metaパッケージのforest()関数のref引数は何ですか？

meta パッケージの forest() 関数における ref 引数は、フォレストプロットの参照線の位置 を指定するためのものです。通常、効果量の基準値（例えば 1 や 0）を設定します。

ref 引数の使い方

forest(meta_analysis, xlab = "Risk Ratio", ref = 1)
  • ref = 1 → リスク比 (RR) の場合、1 を基準線として設定
  • ref = 0 → 標準化平均差 (SMD) の場合、0 を基準線として設定

ポイント

  • ref を指定すると、フォレストプロットに 垂直の基準線 が追加されます。
  • オッズ比やリスク比 の場合は ref = 1 が一般的。
  • 標準化平均差 (SMD) の場合は ref = 0 を使用。

詳しい解説は こちら で確認できます。試してみたい設定はありますか？

フォレストブロットでX軸の無効果の位置に垂直線を引く時に値を指定するのが一般的な使い方の様です。

今回は、Copilotを利用することで、スクリプトをタイプする必要が無いというところは大きな利点であることを再確認しましました。

オッズ比、ハザード比、平均値差、標準化平均値差などほかの効果指標の場合はどうなるか、Copilotに聞いたら、スクリプトを提案してくれるでしょう。

追加：metaパッケージのforest()関数ではlayoutにJAMAのスタイル、RevMan5のスタイルを指定できます。

windows()
forest(meta_analysis, xlab = "Risk Ratio", ref = 1, layout="JAMA")

windows()
forest(meta_analysis, xlab = "Risk Ratio", ref = 1, layout="RevMan5")

上がJAMAのスタイル、下がRevMan5のスタイルで、デフォルトとはデザインが違います。

CopilotでRのスクリプトを作成する：相関、散布図、直線回帰、重回帰

今回は、CopilotにRのスクリプトを提案してもらいましょう。

青のボックス内に書いてあるのがプロンプトです。右端のアイコンをクリックするとCopilotの回答を表示したり隠したりできます。その下にはRでのスクリプトの実行結果を示します。

まず、Rの標準データセットで相関に関するデータサンプルがあるか聞いてみます。

Rの標準データセットに相関に関するサンプルデータはありますか？

はい、Rには相関を調査するのに適した標準データセットがいくつか含まれています。以下のようなデータセットがあります。

  • mtcars：32種類の車に関する燃費や馬力、重量などのデータが含まれており、これらの変数間の相関を分析するのに便利です。
  • iris：アヤメの種類ごとの花弁の長さや幅などのデータがあり、数値変数間の関係性を調べるのに適しています。
  • airquality：大気汚染に関するデータで、温度や風速とオゾン濃度の相関を検討できます。
  • mtcars：自動車の仕様データセットで、燃費と車両重量の関係などを分析できます。

例えば、mtcarsデータセットの変数間の相関係数を求めるには、以下のように記述できます。

data(mtcars)  # データセットをロード
cor(mtcars)   # 相関係数を計算

特定の変数間の相関を調べたい場合は、例えばmpgとhpの関係を確認するには以下のようにします。

cor(mtcars$mpg, mtcars$hp)

他にも、cor.test()関数を使うことで、相関の統計的有意性を検定できます。

興味のあるデータセットがあれば、詳しく説明できますので、お知らせください！

特に指示はしていないのですが、データセットの解説と、相関係数の計算のスクリプトが提示されました。それをRのエディタにコピー・貼り付けて実行すると、以下の結果が、コンソールに出力されました。32種類の車の燃費、馬力などの値について、cor()関数での演算結果が、相関行列として得られました。

次に、散布図を描画してみようと思います。馬力hpと燃費mpgの関係を見たいと思います。グラフ作成にはggplot2というパッケージが良く使われているので、それを使うよう指示します。

mtcarsのhpをX軸に、mpgをY軸に散布図をggplot2で描くスクリプトは？

ggplot2を使ってmtcarsデータセットのhpをX軸、mpgをY軸にした散布図を描くには、以下のRスクリプトを使用できます。

# ggplot2パッケージを読み込む
library(ggplot2)

# mtcarsデータセットを使用して散布図を描画
ggplot(mtcars, aes(x = hp, y = mpg)) +
  geom_point() +                      # 散布図のポイントを追加
  labs(title = "Horsepower vs. Miles Per Gallon", 
       x = "Horsepower (hp)", 
       y = "Miles Per Gallon (mpg)") +
  theme_minimal()                      # シンプルなテーマを適用

このスクリプトでは、ggplot()を使ってデータを指定し、geom_point()で散布図を作成しています。また、labs()を使ってタイトルや軸ラベルを設定し、theme_minimal()でシンプルなデザインにしています。

もし回帰線を追加したい場合は、geom_smooth()を使うと便利です。

ggplot(mtcars, aes(x = hp, y = mpg)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "blue") +
  labs(title = "Horsepower vs. Miles Per Gallon with Regression Line",
       x = "Horsepower (hp)",
       y = "Miles Per Gallon (mpg)") +
  theme_minimal()

このコードでは、geom_smooth(method = "lm")を追加することで、線形回帰をプロットしています。

カスタマイズの方法や他の描画スタイルについても、お手伝いできますので、気軽に聞いてくださいね！

散布図描画のための、スクリプトが提示されました。まず、それをRのエディタにコピー・貼り付けて実行すると、以下のグラフが作成されました。

これも特にプロンプトでは指示していないのですが、回帰直線を描画に追加するスクリプトが提示されたので、それを実行すると以下のグラフが作成されました。この場合、前のグラフに上書きされます。

それでは、hpとmpgの散布図とhpとcylの散布図を重ねて描画するにはどうするかCopilotに聞いてみましょう。

hpとmpgの散布図の回帰直線のグラフに重ねてhpとcylの散布図と
回帰直線を描画するには？mpgとcylではポイントと直線の色を変えてください。

ggplot2を使ってhpとmpg、hpとcylの散布図と回帰直線を重ねて描画するには、aes()のcolorオプションを利用してポイントと回帰線の色を変更すると便利です。以下のRスクリプトで実現できます。

# 必要なパッケージを読み込む
library(ggplot2)

# グラフの作成
ggplot(mtcars, aes(x = hp)) +
  # mpgの散布図と回帰直線（青色）
  geom_point(aes(y = mpg, color = "mpg"), alpha = 0.7) +
  geom_smooth(aes(y = mpg, color = "mpg"), method = "lm", se = FALSE) +
  # cylの散布図と回帰直線（赤色）
  geom_point(aes(y = cyl, color = "cyl"), alpha = 0.7) +
  geom_smooth(aes(y = cyl, color = "cyl"), method = "lm", se = FALSE) +
  # カスタムの凡例タイトル
  scale_color_manual(values = c("mpg" = "blue", "cyl" = "red")) +
  labs(title = "Horsepower vs. MPG and Cylinders",
       x = "Horsepower (hp)",
       y = "Values",
       color = "Variable") +
  theme_minimal()

このスクリプトでは：

  • aes(y = mpg, color = "mpg") と aes(y = cyl, color = "cyl") を使い、mpgを青、cylを赤で描画しています。
  • geom_smooth() でそれぞれ回帰直線を追加。
  • scale_color_manual(values = c("mpg" = "blue", "cyl" = "red")) により、凡例の色を指定。

このコードを実行すると、hpをX軸に、mpgとcylをY軸にした散布図が表示され、色分けされた回帰直線が重ねられます！

他にカスタマイズしたい点があれば、遠慮なく聞いてくださいね。

Copilotが提示したスクリプトをRのエディタにコピー・貼り付けしてそのまま実行すると以下のグラフが作成されました。

馬力hpが大きくなるほど1ガロン当たりの走行距離mpgが短くなり（青）、馬力hpが大きいほどエンジンのシリンダー数cylが多くなる関係があることが分かります。

それでは、直線回帰の切片と係数を出力するにはどのようなスクリプトを書けばいいのか聞いてみたいと思います。まず、hpとmpgの2つの変数についてです。

hpとmpgの間の直線回帰の切片の値と傾きの係数の値を得るにはどの
ようなスクリプトを書けばいいですか？X軸にhp、Y軸にmpgとしてください。

lm() 関数を使って、hp を X 軸、mpg を Y 軸とする線形回帰の切片（intercept）と傾き（coefficient）を求めるには、以下のスクリプトを使用できます。

# mtcarsデータセットをロード
data(mtcars)

# 線形回帰モデルを作成（mpgを目的変数、hpを説明変数）
model <- lm(mpg ~ hp, data = mtcars)

# 回帰係数（切片と傾き）を取得
coefficients <- coef(model)

# 係数の表示
print(coefficients)

説明：

  • lm(mpg ~ hp, data = mtcars) で mpg を目的変数（Y 軸）、hp を説明変数（X 軸）として線形回帰モデルを作成します。
  • coef(model) を使って切片（intercept）と傾きの係数（slope）を取得します。
  • print(coefficients) で結果を表示。

このスクリプトを実行すると、出力される値は以下の形式になります。

(Intercept)          hp  
  X.XXXXXXX      Y.YYYYYY

ここで (Intercept) は切片の値、hp は傾きの係数を示します。

もし、回帰モデルの詳細な統計情報も確認したい場合は、summary(model) を実行すると、p 値や決定係数 (R-squared) なども確認できます。

summary(model)

カスタマイズや可視化の方法など、さらに詳しく知りたいことがあれば、お気軽に聞いてください！

lm()関数を用いて、直線回帰の分析結果を変数modelに格納し、coef(model)で切片と傾きが得られることが分かります。このスクリプトをそのまま実行すると以下の結果がコンソールに出力されました。summary(model)でp値、決定係数も出力されています。

それでは、燃費mpgを従属変数として、それ以外の変数を説明変数にして重回帰分析をするにはどのようなスクリプトを書けばいいかCopilotに聞いてみましょう。

mpgを従属変数として、それ以外の変数を説明変数にして重回帰分析を
するにはどのようなスクリプトを書けばいいですか？切片、係数、p値、
決定係数の値も出力するようにしてください。

選lm() 関数を使って mpg を従属変数（目的変数）として、それ以外の変数を説明変数（独立変数）に設定し、重回帰分析を行い、切片・回帰係数・p値・決定係数を出力するスクリプトは以下のようになります。

# 必要なパッケージをロード
data(mtcars)

# 線形回帰モデルを作成（mpgを従属変数、その他の変数を説明変数に設定）
model <- lm(mpg ~ ., data = mtcars)

# 回帰モデルの概要を表示
summary_model <- summary(model)

# 切片と回帰係数を取得
coefficients <- coef(model)

# p値を取得
p_values <- summary_model$coefficients[,4]

# 決定係数を取得
r_squared <- summary_model$r.squared

# 結果を出力
print("回帰係数（切片と傾き）:")
print(coefficients)

print("p値:")
print(p_values)

print("決定係数 (R-squared):")
print(r_squared)

説明：

  • lm(mpg ~ ., data = mtcars) で mpg を従属変数とし、他の変数を説明変数として重回帰分析を実行します。
  • summary(model) を取得し、回帰係数・p値・決定係数 (R-squared) を抽出します。
  • summary_model$coefficients[,4] で p値を取得できます（各説明変数に対する統計的有意性）。
  • summary_model$r.squared で決定係数を取得し、モデルの適合度を確認できます。
  • print() により結果を表示。

このスクリプトを実行すると、各変数の回帰係数、p値、決定係数の値が出力されます。モデルの評価や変数の選択についても相談できますので、気軽に聞いてくださいね！択についても相談できますので、気軽に聞いてくださいね！

これもそのままコピー・貼り付けしてRで実行した結果が以下の通りです。決定係数は0.869です。

それでは、最後に燃費mpgを従属変数として、説明変数にhp、wt、qsec、amの4つだけでモデルを作る場合のスクリプトをCopilotに書いてもらいます。これら4つの説明変数は上記のp値が小さいものを選びました。

mpgを従属変数として、説明変数にhp、wt、qsec、amの4つだけで
モデルを作る場合のスクリプトは？

lm()関数を使い、mpgを従属変数として、hp、wt、qsec、amの4つを説明変数に設定して重回帰分析を行う場合、以下のスクリプトを使用できます。

# mtcarsデータセットをロード
data(mtcars)

# 線形回帰モデルを作成（mpgを従属変数、hp, wt, qsec, am を説明変数として設定）
model <- lm(mpg ~ hp + wt + qsec + am, data = mtcars)

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

説明：

  • lm(mpg ~ hp + wt + qsec + am, data = mtcars) で mpg を従属変数、hp、wt、qsec、am を説明変数に設定し、重回帰分析を実行します。
  • summary(model) を実行すると、回帰係数（切片と各変数の傾き）、p値、決定係数などの統計情報が表示されます。

モデルの評価や変数の選択についてさらに詳しく知りたい場合は、お気軽に聞いてください！

ここでは、model = lm(mpg ~ hp + wt + qsec + am, data = mtcars)という記述方法に注目してください。多変量のモデルの指定方法です。これらのスクリプトを、Rで実行した結果は以下の通りでした。

各変数の意味を表にしてもらうのをCopilotに頼んでみましょう。

mtcarsの各変数が何を表しているか表を作成してください。

もちろんです！mtcars データセットに含まれる各変数の意味を表にまとめました。

変数名説明
mpg燃費 (Miles Per Gallon)
cylシリンダー数
disp排気量 (cubic inches)
hp馬力 (Horsepower)
dratリアアクスル比 (Rear Axle Ratio)
wt車両重量 (1000ポンド単位)
qsec1/4マイル走行時間（秒）
vsエンジンの形状 (0 = V型, 1 = 直列型)
am変速機の種類 (0 = オートマ, 1 = マニュアル)
gearギア数
carbキャブレター数

このデータセットは、自動車の性能を分析するために広く使われています。特に、mpg (燃費) と他の変数の関係を調べることで、燃費に影響を与える要因を分析できます。

他に詳しく知りたいことがあれば、気軽に聞いてくださいね！

これらの変数の意味から結果を考察すると、馬力と車両重量に対する係数はマイナスなので、これらの値が大きいほど燃費は悪くなる、オートマ（0）に対してマニュアル（1）は係数がプラスなので、マニュアルの方が燃費がいい、1/4マイル走行時間、すなわち0-4の時間が長く、加速が悪いほど燃費がいい、ということが分かります。調整決定係数が0.8368でこれら4つの因子で燃費の大部分が決まるということが分かります。

医学研究でも相関、直線回帰、重回帰分析が適用される例は多いと思います。また、散布図を描いて視覚的に検討することも多いと思います。今回の例でも、Copilotに聞きながらRで解析を進めていくのは有用と思わされました。今回の例は医学の例ではありませんが、応用は容易と思います。

CopilotはRのスクリプトの部分は、コピー用のボタンがついており、それをクリックして、Rで貼り付けるだけですぐ実行できます。つまり、スクリプトを自分で書き込む必要が無くなり、効率化が図れます。

多項ロジスティック回帰分析 – 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()関数を使用して参照基準を明示的に設定することができます。
—————

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

さらに、名義変数の場合、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に分からないことを聞きながら、そのまま実行可能なスクリプトが提示されるので、素晴らしいと思います。

多項ロジスティック回帰分析 – 1 – 臨床医学研究における活用

今回もAIについての話題です。ブラウザのChromeからGoogleのGemini 2.5 Flashを開き、「臨床医学研究で多項ロジスティック回帰分析はどのように利用されていますか？」と聞いてみました。その回答をそのままGammaにテキストとして入力して、8枚のスライドを作成し、PDFファイルとしてダウンロードしました。GeminiとGamma、2つのAIであっという間にこのような資料が作成できます。PPTファイルも作成できます。

最初のスライドです：

多重ロジスティク回帰分析では従属変数が二値である事象が起きる/起きないですが、多項ロジスティック回帰分析 Multinomial logistic regression analysisは従属変数が3つ以上の値をとる名義変数で、説明変数は複数の連続変数、名義変数を設定でき、それらの説明変数の値により3つ以上の名義変数のそれぞれの起きる確率を計算できます。臨床医学ではさまざまな活用シーンが考えられます。

PDFファイルは、ここから開いてみることができますLink 。

多項ロジスティック回帰分析には、Rのパッケージnnetなどを使うことができます。

次の投稿では、スクリプトやデータ形式などについてGeminiに聞きながらその実用性を試そうと思います