Stanを使うベイジアンアプローチのメタアナリシス:Copilotに聞きながら

Copilotに聞きながら、ベイジアンアプローチのメタアナリシスを実行するスクリプトやコードを書けるか試してみたいと思います。

ベイジアンアプローチによる統計解析は、OpenBUGSJAGSでMarkov Chain Monte Carlo (MCMC) simulationを行って、データから事後分布を推定する解析手法です。頻度論派frequentistの伝統的な統計学による方法よりも、より簡単に複雑なモデルを扱うことができるとされています。MCMCの方法はいくつか考案されていますが、OpenBUGSやJAGSの方法はGibbs samplerです。

メタアナリシスにおいても、今までOpenBUGSやJAGSを使う方法がペア比較のメタアナリシス、ネットワークメタアナリシスで利用されてきました。例えば、ネットワークメタアナリシス用のRのパッケージであるgemtcpcnetmetaはJAGSを用いています。Rとの連携はrjagsが使われており、フォレストプロットなども出力できます。

それではまず、OpenBUGSやJAGS以外にRでMCMCを実行する方法があるか聞いてみましょう。

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

RでOpenBUGSやJAGS以外にベイジアンアプローチの解析でMarkov Chain Monte Carlo (MCMC) simulationを実行する方法はありますか?

はい、OpenBUGSやJAGS以外でも、RでMCMCシミュレーションを実行する方法はいくつかあります。代表的なものを紹介します:

  1. Stan (rstan パッケージ)
    • Hamiltonian Monte Carlo (HMC) を使用し、より効率的なMCMCサンプリングを実現します。
    • rstan を使えば、StanのモデルをRで簡単に扱えます。
    • 収束が速く、複雑なベイジアンモデルにも適しています。
  2. MCMCpack (MCMCpack パッケージ)
    • 一般的なベイジアン統計モデルのMCMC推定を簡単に実行できます。
    • 事前分布を指定し、MCMCを用いた推定を行うことが可能です。
  3. nimble (nimble パッケージ)
    • BUGSのようなモデル記述をサポートしつつ、カスタムMCMCアルゴリズムを作成できます。
    • 高速な計算が可能で、柔軟なモデリングができます。
  4. BayesianTools (BayesianTools パッケージ)
    • MCMCのほか、SMC(Sequential Monte Carlo)などのベイジアン推論手法を提供します。
    • モデル適合性の評価や事後分布の解析が容易です。

Stanは特に強力で、最近のベイジアン解析では広く使われています。どれか試してみたいものはありますか?
詳細なチュートリアルやサンプルコードが必要なら、お知らせください!

Rでrstanを介して、Stanを動かし、Gibbs samplerではなく、Hmiltonian Monte Carlo (HMC)を使ったMCMCで、より効率的なサンプリングができるらしいことが分かりました。

そこで、Stanをメタアナリシスに使用した報告があるか聞いてみましょう。

メタアナリシス、ネットワークメタアナリシスでStanを使った報告は多いですか?

はい、Stanはメタアナリシスやネットワークメタアナリシスで広く使用されています。Stanの公式ユーザーガイドにもメタアナリシスのセクションがあり、ベイジアンアプローチの適用方法が詳しく説明されています。

また、Stanを用いたネットワークメタアナリシスに関する報告も増えており、特に従来のWinBUGSやJAGSと比較して、より高速で収束しやすい点が評価されています。日本語の解説もあり、Stanを使ったネットワークメタアナリシスの実装例が紹介されています。

StanはHMC(Hamiltonian Monte Carlo)を採用しているため、従来のMCMCよりも効率的にサンプリングできるのが特徴です。これにより、複雑なメタアナリシスモデルでも安定した推定が可能になります。Stanを使ったメタアナリシスに興味があれば、公式ドキュメントや関連論文をチェックしてみると良いですね!
何か具体的なモデルについて知りたいことはありますか?

参照できるような情報がかなりあるようです。そして、StanはOpenBUGSやJAGSよりも高速で収束しやすいという利点があるようです。そこで、RでStanを動かすために、必要なパッケージをインストールするためのスクリプトを提案してもらおうと思います。すでにインストール済みのパッケージは除き、OSはWindows11にします。

RでStanを使ったベイジアンアプローチのメタアナリシスをしたいと思います。 インストール済みのパッケージ以外のインストールをJapanのミラーサイトから行うRのスクリプトは? OSはWindows11です。

Windows 11環境で、すでにインストール済みのパッケージを確認しつつ、足りないものを日本のCRANミラーサイトからインストールするRのスクリプトはこちらです:

# 必要なパッケージのリスト
packages <- c("rstan", "Rcpp", "ggplot2", "coda")

# CRANの日本ミラーサイトを指定
chooseCRANmirror(graphics = FALSE, ind = 68)  # 68番はJapanのミラー

# インストールされていないパッケージを探す
missing_packages <- setdiff(packages, rownames(installed.packages()))

# 足りないパッケージだけインストール
if(length(missing_packages) > 0) {
    install.packages(missing_packages)
}

# Stanコンパイラの設定(初回のみ)
Sys.setenv(USE_CXX14 = 1)

このスクリプトでは、installed.packages()を使って既にインストール済みのパッケージを確認し、未インストールのものだけをCRANの日本ミラーサイトからインストールするようにしています。
これでスムーズにStanを使ったベイジアンメタアナリシスの環境が整うはずです。何か追加の設定が必要になれば、気軽に聞いてください!

提案されたスクリプトをそのままコピーして、実行したところ、問題なく必要なパッケージをインストールできました。

Stanコンパイラの設定(初回のみ)というスクリプトがあります。このスクリプトの意味を聞いてみましょう。

Stanコンパイラの設定(初回のみ)のスクリプトの目的は?

Sys.setenv(USE_CXX14 = 1) の設定は、Stanを正しくコンパイルするために必要な環境変数を設定する役割を果たします。

StanはCで書かれた統計モデリング言語であり、rstan パッケージを使用する際には、適切なCコンパイラの設定が重要になります。特に、Windows環境では適切なコンパイルオプションを指定しないと、Stanモデルのコンパイル時にエラーが発生する可能性があります。

このスクリプトの目的は:

  1. C++14を有効化 – Stanの最新バージョンではC++14以降が必要なため、環境変数 USE_CXX14=1 を設定します。
  2. モデルのスムーズなコンパイル – この設定を行うことで、Stanコードを正しくコンパイルできるようになります。

この設定は通常、Rセッションの最初に一度だけ行えば十分です(毎回実行する必要はありません)。ただし、環境設定がリセットされた場合や、新しいRセッションを開始した場合には再設定が必要になることもあります。

もしコンパイル時に問題が発生する場合は、rstan の公式ドキュメントを確認するか、追加の設定が必要かどうか調べるのも良いですね。何か困ったことがあれば、いつでも聞いてください!

Copilotの回答は、StanはC言語で書かれており、StanのモデルをコンパイルするためのC++14を有効化することだという説明です。また、このスクリプトの実行は初回のみ必要だということです。Rでrstanをインストールすると、Stanも一緒にインストールされるようです。私の場合は、rstan以外のパッケージはインストール済みでしたので、rstanだけが新たにインストールされました。

それでは、二値変数アウトカムで、ペア比較メタアナリシスを実行し、リスク比の統合値と95%信頼区間を計算するスクリプトを提案してもらいましょう。変数名は指定します。

rstanをインストール後、二値変数アウトカムの場合のペア比較メタアナリシスを実行したいと思います。簡単なサンプルデータと、リスク比の統合値と95%信頼区間を計算するスクリプトを提案してください。研究IDの変数はstudy、対照群の症例数はnc、イベント数はrc、介入群の症例数はnt、イベント数はrtとします。

了解しました。以下のスクリプトでは、rstan を使用してペア比較メタアナリシスを実行し、リスク比 (Risk Ratio, RR) の統合推定値とその95%信頼区間を計算します。


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

まず、指定された変数 study, nc, rc, nt, rt を含むデータを作成します。

# 必要なパッケージをロード
library(rstan)

# サンプルデータ
meta_data <- list(
  N = 5,  # 研究数
  study = 1:5,  # 研究ID
  nc = c(100, 200, 150, 120, 250),  # 対照群の症例数
  rc = c(10, 30, 20, 15, 25),  # 対照群のイベント数
  nt = c(100, 200, 150, 120, 250),  # 介入群の症例数
  rt = c(15, 40, 22, 18, 35)  # 介入群のイベント数
)

2. Stanモデルの定義

Stanのコードを記述し、リスク比のログ値 (logRR) を推定します。

data {
  int<lower=1> N;  
  int nc[N];  
  int rc[N];  
  int nt[N];  
  int rt[N];  
}

parameters {
  real logRR[N];  # リスク比の対数値
  real mu;        # リスク比の統合値(対数)
  real<lower=0> sigma; # 研究間のばらつき
}

model {
  # 事前分布
  mu ~ normal(0, 10);
  sigma ~ cauchy(0, 2);

  # 各研究のリスク比推定
  for (i in 1:N) {
    logRR[i] ~ normal(mu, sigma);
    rc[i] ~ binomial(nc[i], inv_logit(0));  # 対照群
    rt[i] ~ binomial(nt[i], inv_logit(logRR[i]));  # 介入群
  }
}

3. Stanモデルの実行

RからStanモデルを実行し、リスク比とその95%信頼区間を計算します。

# Stanモデルのコンパイルと実行
fit <- stan(
  model_code = stan_model_code, 
  data = meta_data,
  iter = 2000, chains = 4
)

# 推定結果の抽出
posterior_samples <- extract(fit)
mu_samples <- posterior_samples$mu

# リスク比 (exp(mu)) と 95% 信頼区間の計算
risk_ratio <- exp(mean(mu_samples))
conf_int <- exp(quantile(mu_samples, probs = c(0.025, 0.975)))

# 結果の表示
cat("推定されたリスク比:", risk_ratio, "\n")
cat("95% 信頼区間:", conf_int[1], "-", conf_int[2], "\n")

まとめ

このスクリプトでは、
ペア比較メタアナリシスをStanで実行
リスク比 (RR) の統合値を推定
95% 信頼区間の計算

Stanは高速なMCMCを用いるため、収束が速く、従来のBUGS系よりも精度の高い推定が可能です。
さらにカスタマイズしたい点があれば、教えてください!

まずは、提案されたRのスクリプトをそのまま実行してみようと思いますが、実際やってみると、2. Stanモデルの定義のスクリプトの実行でエラーが起きました。そのエラーの最初の部分だけを下に示します。

この部分のスクリプトはモデルの記述で、おそらくStanで実行されるdataとparametersとmodelが記述されている部分だろうと思われます。Stanが実行するコードですから、それをそのままRで実行しようとしてもエラーが起きるのは当然と言えるでしょう。

そこで、正しい記述法をCopilotに提案してもらいましょう。

Stanのモデルの定義をRでそのまま実行するとエラーが起きます。data, parameters, modelの正しい記述法は?

StanのコードをRで正しく扱うためには、Stanのモデルを.stanファイルに保存するか、文字列としてR内に記述する必要があります。
また、stan()関数でモデルを実行する際には、適切なデータ構造を指定することも重要です。


1. Stanコードの記述

Stanのモデル定義は 文字列として保存 し、stan_model_code に格納します。

stan_model_code <- "
data {
  int<lower=1> N;
  int nc[N];  
  int rc[N];  
  int nt[N];  
  int rt[N];  
}

parameters {
  real logRR[N];  
  real mu;        
  real<lower=0> sigma;
}

model {
  mu ~ normal(0, 10);
  sigma ~ cauchy(0, 2);
  
  for (i in 1:N) {
    logRR[i] ~ normal(mu, sigma);
    rc[i] ~ binomial(nc[i], inv_logit(0));  
    rt[i] ~ binomial(nt[i], inv_logit(logRR[i]));  
  }
}
"

この方法でStanコードをR内に保存できます。


2. データの準備

Stanモデルに渡すデータは、リスト形式 にする必要があります。

meta_data <- list(
  N = 5,
  nc = c(100, 200, 150, 120, 250),  
  rc = c(10, 30, 20, 15, 25),  
  nt = c(100, 200, 150, 120, 250),  
  rt = c(15, 40, 22, 18, 35)  
)

3. Stanモデルの実行

stan()関数を使ってモデルを実行します。

library(rstan)

fit <- stan(
  model_code = stan_model_code, 
  data = meta_data,
  iter = 2000, chains = 4
)

print(fit)

エラー回避ポイント

  1. Stanコードは文字列として格納 (stan_model_code <- "...")
  2. データはリスト形式 (meta_data <- list(...))
  3. Stanのモデルはstan()関数で実行

この方法でエラーなくStanのモデルを実行できるはずです。
もしまだエラーが出る場合は、詳細を教えてください!

回答では、サンプルデータの作成のスクリプトは同じままでした。Stanコードの記述がstan_model_codeという変数に文字列として格納するように書き換わっています。Stanで実行されるコードはdata, parameters, modelの3つの要素からなり、それぞれテキストとして記述されるということが分かります。またmodelの部分をテキストファイルとして保存して、実行の際にそれを読み込んで使用することもできるようですが、それは後で試すことにして、Rのスクリプトの中で、stan_model_codeという変数に直接格納することにします。

1.Stanコードの記述の部分を見ると、dataの部分はintegerすなわち、整数で配列構造で整数Nが順番を表すという形式であることを示しています。parametersの部分はmodelでの計算に使われる変数を実数として宣言しています。そしてmodelでは、muは平均値0、標準偏差10の正規分布、sigmaは位置0、スケール2のCauchy distributionコーシー分布からサンプリングすることが定義されています。これら事前分布はかなり幅の広い、あいまい分布になります。

さらに、log(RR)が平均値mu, 標準偏差sigmaの正規分布からサンプリングすること、rc[i], rt[i]すなわちイベント数は二項分布からサンプリングすることが定義されています。

Stanモデルの実行のスクリプトは最初の提案と同じですが、今回のCopilotの回答には、推定結果の抽出と結果の表示の部分のスクリプトは含まれていませんでした。

  1. それではまず、Stanコードの記述の部分を実行してみます。今度はエラーは起きません。つまり、Stanで使われる、stan_model_codeという変数に、data, parameters, modelの3要素が格納されたようです。

modelの設定に続いて、最初のCopilotの提案の、Stanモデルの実行のスクリプトを実行してみます。やはりエラーが起きてしまいました。WARNINGが出力され、Rtoolsをインストールする必要があるとのメッセージです。

そこでWARNINGの中ののURL https://cran.r-project.org/bin/windows/Rtools/を開いてみたところが以下の図です。Rのバージョンによって異なるプログラムが用意されています。

今このPCで使用しているRはVersion 4.4.3なので、RTools 4.4をインストールしてみます。RTools 4.4をクリックすると以下の解説があります。

インストールを実行する前にCopilotに聞いて確認しておきましょう。

Stanモデルの実行で、RToolsのインストールが必要だというWARNINGが出ました。

StanのモデルをWindows環境で実行するには、RTools のインストールが必要です。RToolsは、StanのC++コンパイルをサポートするツールチェーンであり、rstan を正しく動作させるために必要になります。

RToolsのインストール方法

  1. RToolsの公式サイトから最新版をダウンロードします → こちら
  2. インストーラーを実行し、デフォルト設定でインストールします。
  3. インストール後、Rを再起動し、以下のコマンドでRToolsが正しく認識されているか確認します: system("g++ -v") # C++コンパイラのバージョン確認 system("where make") # makeコマンドの確認 これらのコマンドが正常に動作すれば、RToolsのインストールは成功しています。

RToolsをStanで使用するための設定

Stanを正しく動作させるために、以下の設定を行います:

install.packages("Rcpp", type = "source")  # C++インターフェース
install.packages("inline", type = "source")  # Stanのコンパイル補助

これでStanのモデルをスムーズに実行できるはずです。
もし追加のエラーが発生した場合は、詳細を教えてください!

以上を確認したうえで、RTools installer 4.4をダウンロードし、デフォルトの設定でインストールを実行しました。431MBのかなり大きなプログラムなのでダウンロードには少し時間がかかりました。RToolsのインストーラーのダウンロードが終了したら、いったんRを閉じて、ダウンロードしたrtools44-6459-6401.exeファイルをダブルクリックしてインストールしました。

Copilotの回答には、Rを再起動し、system(“g++ -v”)およびsystem(“where make”)を実行してインストールが成功したか確認するとあったので、それを実行してみました。

特にエラーメッセージは出ませんでした。

そこで、再度最初から以下を実行してみました。

サンプルデータmeta-dataの作成、stan_model_codeの作成、そして続けて、Stanモデルのコンパイルと実行です。その結果、コンソールにはChain 1~4までのWarmup、Samplingの結果が出力されました。Chain 1の部分を以下に示します。一つのチェーンで最初の1000回はウォームアップ、残り1000回はサンプリングのデータを保存したということです。チェーン1で0.101秒かかっていました。


上記の最初のCopilotの提案の内、3.Stanモデルの実行の# 推定結果の抽出、# リスク比 (exp(mu)) と 95% 信頼区間の計算、# 結果の表示の部分を実行しました。

リスク比の統合値と95%信頼区間の値が出力されていました。

しかし、このリスク比の統合値は正しいとは思えません。なぜなら、解析対象のデータを見ると、nc, ntつまり対照群と介入群の症例数は同数で、どの研究もrt > rcつまりイベント数は介入群の方が対照群より大きい値です。ですから、リスク比は1.0より大きくなるはずですが、統合値は0.185となっているので、明らかに間違っています。

そこで、modelをよく見ると、以下の様にrc[i] ~ binomial(n[i], inv_logit(0));の部分がどうもおかしいようです。

model {
mu ~ normal(0, 10);
sigma ~ cauchy(0, 2);

for (i in 1:N) {
logRR[i] ~ normal(mu, sigma);
rc[i] ~ binomial(nc[i], inv_logit(0));
rt[i] ~ binomial(nt[i], inv_logit(logRR[i]));
}
}

Stanのinv_logit()関数についてCopilotに聞いてみました。

Stanのinv_logit()関数はどのような演算をするのでしょうか?

Stanのinv_logit()関数は、ロジスティック関数(シグモイド関数)を計算するための関数です。これは、入力値を0から1の範囲に変換する役割を持ちます。

定義

inv_logit(x) は、以下の式で計算されます:

[ \text{inv_logit}(x) = \frac{1}{1 + e^{-x}} ]

この関数は、ロジスティック回帰やベイジアン統計モデルで確率の推定に広く使われます。

特徴

入力値が負の値 → 出力は 0 に近づく
入力値が 0 → 出力は 0.5
入力値が正の値 → 出力は 1 に近づく

Stanでは、inv_logit() を使って確率を表す変数を定義することが多いです。例えば、ロジスティック回帰モデルでは以下のように使用されます:

parameters {
  real beta;
}

model {
  beta ~ normal(0, 1);
  real p = inv_logit(beta);  // 確率の計算
}

Stanの公式ドキュメントにも詳しい説明がありますので、さらに詳しく知りたい場合はこちらをチェックしてみてください!💡

つまり、inv_logit()関数はロジスティック回帰で従属変数の値から確率の値を求める際の計算式と同じで、1/(1 + exp(-x))を計算する関数です。inv_logit(0)は0.5になり、対照群のイベント率がいつも0.5に固定されてしまいます。

そこで、この間違いを指摘ながら、正しいと思われるスクリプトをCopilotに提案してもらおうと思います。

以下のコードはStanで二値変数アウトカムのメタアナリシスを実行するコードです。ncは対照群の症例数、rcは対照群のイベント数、ntは介入群の症例数、rtは介入群のイベント数です。効果指標はリスク比を用います。このコードの中で、rc[i] ~ binomial(nc[i], inv_logit(0)); は正しいですか?inv_logit(0)は0.5に固定されるので、対照群のイベント率が0.5に固定されてしまいます。
<以下上記のモデルのコードをここに記述>

はい、ご指摘の通り、inv_logit(0)常に 0.5 となるため、このコードでは対照群のイベント率が 0.5 に固定されてしまいます。これは、研究ごとのイベント率のばらつきを考慮していないため、不適切な設定になっています。

修正方法

対照群のイベント率をデータに基づいて推定できるようにするため、rc[i] の分布を binomial() に渡す確率 ( p_c[i] ) をデータから取得する形に修正しましょう。

修正後のStanモデル

data {
  int<lower=1> N;  
  int nc[N];  
  int rc[N];  
  int nt[N];  
  int rt[N];  
}

parameters {
  real logRR[N];  // リスク比の対数値
  real mu;        // リスク比の統合値(対数)
  real<lower=0> sigma; // 研究間のばらつき
  real p_c[N];    // 対照群のイベント率
}

model {
  # 事前分布
  mu ~ normal(0, 10);
  sigma ~ cauchy(0, 2);
  
  # 対照群のイベント率の事前分布(ベータ分布を使用)
  for (i in 1:N) {
    p_c[i] ~ beta(1, 1);  // 一様な事前分布
    logRR[i] ~ normal(mu, sigma);
    rc[i] ~ binomial(nc[i], p_c[i]);  // 修正
    rt[i] ~ binomial(nt[i], p_c[i] * exp(logRR[i]));  // 介入群
  }
}

修正のポイント

  1. イベント率 p_c[i] をパラメータとして推定
    • p_c[i] ~ beta(1, 1);一様な事前分布 を設定(研究ごとの不確実性を考慮)
  2. 対照群のイベント率 p_c[i]binomial() に適用rc[i] ~ binomial(nc[i], p_c[i]);
    • 各研究のイベント率 を適用することで、固定値 0.5 ではなく、データに基づいた推定が可能に
  3. 介入群のイベント率を p_c[i] * exp(logRR[i])rt[i] ~ binomial(nt[i], p_c[i] * exp(logRR[i]));
    • リスク比 exp(logRR[i]) を適用 してイベント率を調整

この修正のメリット

イベント率のばらつきを考慮した推定が可能
対照群のイベント率を固定値 0.5 にする問題を解消
ベイジアンアプローチに適した柔軟な推定ができる

この修正が適切かどうか、試してみてください!もしさらに調整が必要なら、お知らせください💡

この修正されたコードをもう一度実行してみましょう。修正後のStanモデル全体をダブルクォーテーションで囲み、変数stan_model_codeに格納します。

さて、この修正モデルを実行するとstan()関数の実行の所で、初期値が扱える範囲外であるというエラーが多発しました。

そこで、そのエラーメッセージに対する対処をCopilotに聞いて、2回ほど修正を加え、実行可能とだったスクリプトは以下の通りです。この部分のCopilotとのやりとりは省略します。

#Fixed effects model with binoial distributions:

# 必要なパッケージをロード
library(rstan)

# サンプルデータ
meta_data <- list(
  N = 5,  # 研究数
  study = 1:5,  # 研究ID
  nc = c(100, 200, 150, 120, 250),  # 対照群の症例数
  rc = c(10, 30, 20, 15, 25),  # 対照群のイベント数
  nt = c(100, 200, 150, 120, 250),  # 介入群の症例数
  rt = c(15, 40, 22, 18, 35)  # 介入群のイベント数
)

stan_model_code <- "data {
  int<lower=1> N;  
  int nc[N];  
  int rc[N];  
  int nt[N];  
  int rt[N];  
}

parameters {
  real logRR[N];  // リスク比の対数値
  real mu;        // リスク比の統合値(対数)
  real<lower=0> sigma; // 研究間のばらつき
  real<lower=0, upper=1> p_c[N];  // 確率値の制約を追加
}

model {
  # 事前分布
  mu ~ normal(0, 10);
  sigma ~ cauchy(0, 2);
  
  # 対照群のイベント率の事前分布(ベータ分布を使用)
  for (i in 1:N) {
    p_c[i] ~ beta(1, 1);  // 一様な事前分布
    logRR[i] ~ normal(mu, sigma);
    rc[i] ~ binomial(nc[i], p_c[i]);  // 修正
    rt[i] ~ binomial(nt[i], fmin(p_c[i] * exp(logRR[i]), 1));  // 最大値を1に制限 
  }
}"

rt[i] ~ binomial(nt[i], p_c[i] * exp(logRR[i]));  // 介入群

N=meta_data$N
init_values <- function() {
  list(mu = 0, sigma = 1, logRR = rep(0, N), p_c = rep(0.1, N)) # 確率の初期値を設定
}

# Stanモデルのコンパイルと実行
fit <- stan(model_code = stan_model_code, data = meta_data, init = init_values, iter = 2000, chains = 4)

# 推定結果の抽出
posterior_samples <- extract(fit)
mu_samples <- posterior_samples$mu

# リスク比 (exp(mu)) と 95% 信頼区間の計算
risk_ratio <- exp(mean(mu_samples))
conf_int <- exp(quantile(mu_samples, probs = c(0.025, 0.975)))

# 結果の表示
cat("推定されたリスク比:", risk_ratio, "\n")
cat("95% 信頼区間:", conf_int[1], "-", conf_int[2], "\n")

######

Copilotが提案したコードはFixed effects modelだと記されていますが、modelを見るとlogRR[i] ~ normal(mu, sigma)となっていて、各研究のリスク比の自然対数を平均値mu、標準偏差sigmaの正規分布からサンプリングしているので、これはRandom effects modelのはずです。そして、stan_model_codeに格納されるmodelの記述の中にあるrt[i] ~ binomial(nt[i], fmin(p_c[i] * exp(logRR[i]), 1));がその下に重複して記述されています。これは間違いと思われるので、Copilotに聞いてみましょう。

これはStan, rstanでペア比較メタアナリシスを実行するためのRのスクリプトです。 これはFixed effects modelではなくRandom effects modelではないでしょうか。 ダブルクォーテーションで囲ったmodelの記述の後の以下の部分は必要ですか? rt[i] ~ binomial(nt[i], p_c[i] * exp(logRR[i])); // 介入群
<以下スクリプト全体のコピー>

このStanのコードはRandom effects modelに該当します。理由としては、以下の点が挙げられます:

  • logRR[i] ~ normal(mu, sigma); の部分で、各研究の対数リスク比 logRR[i]全体の平均 mu と研究間のばらつき sigma に従う正規分布から生成されているため、個々の研究のリスク比が異なる可能性を考慮していることになります。これはRandom effects modelの特徴です。
  • Fixed effects modelでは、logRR[i] は各研究で固定された値を持つべきですが、このモデルでは研究間のばらつきを考慮して正規分布に基づく変動を許容しているため、Random effects modelと言えます。

次に、rt[i] ~ binomial(nt[i], p_c[i] * exp(logRR[i])); の部分ですが、これは介入群のイベント発生数 rt[i] を確率 p_c[i] * exp(logRR[i]) に基づく二項分布から生成している処理です。このコードは必要であり、介入群のイベント率をモデル化する重要な式になっています。

ただし、p_c[i] * exp(logRR[i]) が1を超える可能性があるため、後の修正では fmin(p_c[i] * exp(logRR[i]), 1) を使用して確率の上限を1に制限しています。この修正は適切であり、数値的な不安定さを防ぐために重要です。

もしFixed effects modelにしたい場合は、logRR[i] を変数として個別に推定するのではなく、全ての研究で共通の logRR を固定する形に書き換えると良いでしょう。たとえば、logRR[i] ~ normal(mu, 0.001); のように sigma をほぼゼロにすることで、変動を抑えられます。

何か他に修正点や疑問があれば、お気軽に聞いてください!

やはりRandom effects modelです。そして、Fixed effects modelにしたい場合のlogRR[i]のサンプリングの方法が提案されています。

rt[i] ~ binomial(nt[i], fmin(p_c[i] * exp(logRR[i]), 1));の重複については、必要だとの回答ですが、もう一度確認してみましょう。

stan_model_codeに格納される””内のコードにrt[i] ~binomial(nt[i], fmin(p_c[i] * exp(logRR[i]), 1)); // 最大値を1に制限があるので、その下にある同様のコードは不要では?

その通りです!stan_model_code に格納されるStanのコード内で、rt[i] ~ binomial(nt[i], fmin(p_c[i] * exp(logRR[i]), 1)); の記述があるため、 その後の rt[i] ~ binomial(nt[i], p_c[i] * exp(logRR[i])); は不要です

重複したコードを削除しても問題なく動作するはずですし、余計な記述がなくなることでコードの可読性も向上します。特に、確率値の上限を1に制限する修正がすでに適用されているため、後ろのコードをそのまま残すと、制約なしの二項分布に戻ってしまう可能性があるので注意が必要です。

他にもコードの整理や最適化について相談があれば、お気軽に聞いてください!

やはり不要の様です。

また、このコードには初期値に関する設定が含まれていました。それが必要なのか聞いてみましょう。

init_values <- function(){ }の部分は必要ですか?

init_values <- function(){ } の部分は必須ではありませんが、モデルの収束を安定させるために役立つ場合があります。

この関数は、StanのMCMC(マルコフ連鎖モンテカルロ)サンプリングを開始するときに初期値を設定するためのものです。特に、Stanのデフォルトの初期化では、パラメータにランダムな値を割り当てるため、時々収束の問題が発生することがあります。これを防ぐために、init_values を定義して適切な初期値を設定することで、よりスムーズに計算が進む可能性があります。

ただし、もし init_values を削除してもモデルが問題なく収束するなら、必須ではないので省略できます。より安定した推定を得たい場合は、適切な初期値を設定することをおすすめします。

他にも気になる点があれば聞いてください!

確かに、初期値の設定は必須ではないのですが、時にMCMCが収束しない問題を避けるため適切な初期値を設定することが望ましいです。この例では、muに0、sigmaに1、logRRにN個の研究すべて0、p_cすなわち対照群のイベント率にN個の研究すべてに0.1を設定するようになっています。

そのほかスペルミスもあったので、それらも含めて修正したコードを書き出してもらいましょう。

それでは、Random effects modelに修正、binoial distributionsをbinomial distributionsに修正、コードの重複部分を削除したコードを書き出してください。

了解しました!以下のコードでは、Random effects model に修正しbinomial distributions のスペルミスを修正し重複部分を削除しました。

# 必要なパッケージをロード
library(rstan)

# サンプルデータ
meta_data <- list(
  N = 5,  # 研究数
  study = 1:5,  # 研究ID
  nc = c(100, 200, 150, 120, 250),  # 対照群の症例数
  rc = c(10, 30, 20, 15, 25),  # 対照群のイベント数
  nt = c(100, 200, 150, 120, 250),  # 介入群の症例数
  rt = c(15, 40, 22, 18, 35)  # 介入群のイベント数
)

stan_model_code <- "data {
  int<lower=1> N;  
  int nc[N];  
  int rc[N];  
  int nt[N];  
  int rt[N];  
}

parameters {
  real logRR[N];  // リスク比の対数値
  real mu;        // リスク比の統合値(対数)
  real<lower=0> sigma; // 研究間のばらつき
  real<lower=0, upper=1> p_c[N];  // 確率値の制約を追加
}

model {
  // 事前分布
  mu ~ normal(0, 10);
  sigma ~ cauchy(0, 2);
  
  // 対照群のイベント率の事前分布(ベータ分布を使用)
  for (i in 1:N) {
    p_c[i] ~ beta(1, 1);  // 一様な事前分布
    logRR[i] ~ normal(mu, sigma);
    rc[i] ~ binomial(nc[i], p_c[i]);  // 修正
    rt[i] ~ binomial(nt[i], fmin(p_c[i] * exp(logRR[i]), 1));  // 最大値を1に制限 
  }
}" 

N <- meta_data$N
init_values <- function() {
  list(mu = 0, sigma = 1, logRR = rep(0, N), p_c = rep(0.1, N))  # 確率の初期値を設定
}

# Stanモデルのコンパイルと実行
fit <- stan(model_code = stan_model_code, data = meta_data, init = init_values, iter = 2000, chains = 4)

# 推定結果の抽出
posterior_samples <- extract(fit)
mu_samples <- posterior_samples$mu

# リスク比 (exp(mu)) と 95% 信頼区間の計算
risk_ratio <- exp(mean(mu_samples))
conf_int <- exp(quantile(mu_samples, probs = c(0.025, 0.975)))

# 結果の表示
cat("推定されたリスク比:", risk_ratio, "\n")
cat("95% 信頼区間:", conf_int[1], "-", conf_int[2], "\n")

主な修正点:

  1. モデルを Random effects model に変更
    • logRR[i] ~ normal(mu, sigma); を維持し、研究間のばらつきを考慮しました。
  2. スペルミスの修正
    • binoial distributionsbinomial distributions に修正しました。
  3. コードの重複部分を削除
    • rt[i] ~ binomial(nt[i], p_c[i] * exp(logRR[i])); の不要な部分を削除しました。

これで、修正点を反映した適切な Random effects model になっています!
もし追加の調整が必要でしたら、お気軽にご相談ください。

それでは提案されたスクリプトをRで実行してみます。

そして、得られた結果は以下の通りです。

今回のリスク比の統合値は間違いないように見えます。

ただし、# Stanモデルのコンパイルと実行、つまりstan()関数でMCMCを実行した後にサンプリングの数が少ないため平均値や中央値は信頼できませんという警告メッセージが出ていました。

4000のサンプリングでは少ないことは明らかなので、もっと多くしたいと思います。上記のコードではiter 2000, chains 4なのですが、これを大きくすることで対処できるはずです。Warmupの回数も指定したいので、その方法をCopilotに聞いてみましょう。

stan()関数でwarmupの設定はどのようにしますか?

stan() 関数では、warmup の設定を iter と組み合わせて指定します。warmup はMCMCサンプリングの最初の部分で、チェーンの初期化が終わった後に適切なサンプリング状態に到達するまでの試行です。通常、warmup 期間のサンプルは解析には使用されず、モデルの収束を助けるためにのみ使われます。

warmup の設定方法:

fit <- stan(
  model_code = stan_model_code, 
  data = meta_data, 
  init = init_values, 
  iter = 2000,  # 総サンプリング回数
  warmup = 500,  # ウォームアップ期間(初期の500サンプルを破棄)
  chains = 4  # チェーン数
)

この例では iter = 2000 で全体のサンプリングを行い、最初の 500 サンプル (warmup = 500) はウォームアップとして扱われます。通常、warmup の値は 全体の iter の 25~50% 程度に設定するのが推奨 されています。

注意点:

  • warmup を適切に設定することで、収束の質が向上します。
  • warmup が短すぎると、収束が不安定になり、バイアスのある推定になる可能性があります。
  • warmup が長すぎると計算コストが増えますが、適切な値を設定することでパフォーマンスを最適化できます。

最適な warmup の値はモデルによって変わるので、トレースプロット (traceplot(fit)) などを確認して、適切な長さを選ぶのがよいでしょう。

試してみて、もしモデルが収束しない場合は warmup の値を調整しながら最適化するのがオススメです。何か他に調整したい点があれば聞いてください!

Warmupは全体のサンプリング回数の25~50%にすることが推奨されているということなので、iter = 40000, warmup =20000, chains =5で100000個のサンプルを得るようにしてみたいと思います。stan()関数の部分を以下の様に書き換え、 もう一度このスクリプトを実行します。

fit <- stan(model_code = stan_model_code, data = meta_data, init = init_values, iter = 40000, warmup = 20000, chains = 5)

今度はBulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.という警告は出ませんでした。

その結果です。

############################################

おそらくこれであっていると思われますが、頻度論派のメタアナリシスをRのmetaパッケージで実行してみました。結果を比較してみたいと思います。

リスク比の統合値はほぼ同じで、小数点以下1桁まで、1.2で同じですがそれ以下の桁は異なります。95%信頼区間の値は、Stanのベイジアンアプローチの方は0.8854134~1.744159、metaの頻度論派アプローチの方は1.0196~1.6540で前者の方が幅が広めです。Rのパッケージmetaforでもやってみましたが、metaと全く同じ値が得られました。いずれも研究間の分散の計算にはREML (restricted maximum likelihood)法が使われています。

さて、MCMCをstan()関数で実行した結果が変数fitに格納されており、そこからexract(fit)で取り出された値は変数posterior_samplesに格納されています。その構造を調べるため、str(posterior_samples)を実行すると、以下の結果でした。

変数logRR, mu, sigma, p_c, lp_が100000個記録されていることが分かります。

ここにlp_という変数が出てきたのでこれが何か聞いてみましょう。

posterior_samplesのlp_は何ですか?

posterior_samples$lp_ は、Stanの出力に含まれる**対数事後確率(log posterior probability)**の値を指します。

StanのMCMCサンプリングでは、各サンプリングステップで**モデルの対数確率(log probability)**を計算しており、それが lp_ という変数として出力されます。

lp_ の意味:

  • lp_ は、MCMCサンプルごとの**対数確率密度(log probability density)**を表します。
  • lp_ の値がより高いほど、そのサンプルはモデルにとって「より適合した状態」を示します。
  • traceplot(fit, pars = "lp_") を使うことで、lp_ の挙動をプロットし、モデルの収束を確認することができます。

lp_ の使い方:

  1. 収束診断
    lp_ の値が一定の範囲で安定しているかどうかを確認することで、MCMCのチェーンが適切に収束しているか判断できます。
  2. モデル比較
    異なるモデルを比較する際に、lp_ の平均を確認することで、どのモデルの適合度が高いか評価できます。
  3. 異常値の検出
    lp_ の値が極端に低いサンプルは、モデルの異常な状態を示している可能性があるため、適切にフィルタリングすることが重要です。

Stanの出力で lp_ の変動が大きすぎる場合は、モデルの収束が不安定な可能性があるので、ウォームアップの長さ (warmup) や事前分布を調整すると良いでしょう。

もし詳細な活用方法について知りたければ、気軽に聞いてください!

StanのMCMCサンプリングでは、各サンプリングステップで**モデルの対数確率(log probability)**を計算しており、それが lp_ という変数として出力されますとのことでした。収束診断、モデル比較、異常値の検出に使われるとのことです。

##############

今回もCopilotの回答は間違いを含んでいました。Stanをrstanで動作させるには、Rtoolsのインストールが必要なのにそれに言及しなかった、modelはテキストとして変数に格納する必要があることが無視されていた、効果指標として、リスク比を指定したにもかかわらず、オッズ比の例を提案してきた、対照群のイベント率を固定するなどの問題が起きました。しかし、このような間違いや問題はほとんどの場合、実際にRでスクリプトを実行するとエラーが起きるので、それに基づいて修正を提案させることができることも分かりました。Stan, rstanは、今回初めて使用したこと、Stan、rstanのマニュアルや関連した情報を読み込まなくてもできたことを考えると、やはり非常に効率的だと思いました。

Copilotと何回かやりとりしているうちに、Stanのコードの記述法が何となくわかってきました。以下にまとめておきます。

  1. dataブロックでは、解析対象のデータの変数の属性を設定する。
  2. transformed dataブロックでは、各変数間の演算を設定する。
  3. parametersブロックでは平均値などの知りたい値の変数の属性を設定する。
  4. modelブロックでは、パラメータの事前分布の設定と、パラメータのサンプリングを設定する。
  5. MCMCの実行はstan()関数で行い、モデル、データ、サンプリング回数とチェーン回数を指定する。(この実行には少し時間がかかる)。
  6. 結果はstan()関数の実行結果を格納した変数からextract()関数で例えば、posterior_samplesのような変数に取り出し、さらにパラメータの変数名を$で結合して個別に取り出す。
  7. 個別のパラメータのサンプリングに対して、Rの関数を用いて必要な統計値を得る。

R, Stan, brmsを使ったメタアナリシスに関して、以下のサイトから情報が得られます。今回はCopilotに聞きながらどこまでできるかというのがテーマなので、これらは参照しないでやりました:
Stan Development TeamのStan: Software for Bayesian Data AnalysisおよびRSstan: the R interface to Stan
Matti Vuorre氏のMatti’s homepage Bayesian Meta-Analysis with R, Stan, and brms: Meta-analysis is a special case of Bayesian multilevel modeling

それでも何とかベイジアンアプローチのメタアナリシスが実行できました。結果の妥当性についてのチェックは他のデータでも試す必要がありますし、結果をフォレストプロットにしたり、出版バイアスのチェックはどうやるか、など課題は残っています。また、オッズ比、リスク差、平均値差、標準化平均値差の場合やネットワークメタアナリシスへの拡張なども課題です。

追記:
Warn DEらの論文で記述されているランダム効果モデルによるメタアナリシスのためのWinBUGS用のコードを用い、同じサンプルデータでOpenBUGSで分析を実行したところ、RRの統合値は以下のようにStanの場合と近似した値でした。

リスク比の統合値は1.246 (95% CI 0.8901~1.731)で上記の結果に近似していました。こちらのモデルはやはり二項分布を用いており、今回のStanのモデルと同じと言えると思います。

BUGSコードの書き方はStanの場合と同じではありません。BUGSでは正規分布の指定に平均値と精度(分散の逆数)を使いますが、Stanは平均値と標準偏差を使います。それ以外にデータ処理については、例えば、RRをlogRRに相当する変数deltをサンプリングする時点で計算し、RRのサンプリングも記録しています。Stanのモデルの方ではRRはサンプリングはしないで、logRRのサンプルから後で計算しています。しかし、rc、rtを二項分布からサンプリングしていることは共通です。

BUGSのコードをStan用に書き換えるのも今後の課題です。もしCopilotやGeminiが正しく変換してくれると楽で良いのですが。


Warn DE, Thompson SG, Spiegelhalter DJ: Bayesian random effects meta-analysis of trials with binary outcomes: methods for the absolute risk difference and relative risk scales. Stat Med 2002;21:1601-23. doi: 10.1002/sim.1189 PMID: 12111922

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 引数は、フォレストプロットの参照線の位置 を指定するためのものです。通常、効果量の基準値(例えば 10)を設定します。

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)   # 相関係数を計算

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

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を使ってhpmpghpcylの散布図と回帰直線を重ねて描画するには、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軸に、mpgcylを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を従属変数として、hpwtqsecamの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 を従属変数、hpwtqsecam を説明変数に設定し、重回帰分析を実行します。
  • 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)で何種類含まれているかを知ることができるようです。

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

名義変数のそれぞれの値は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を推定する方法、スクリプトについても説明があります。以下がその部分です。臨床の場合であれば、個々の対象者のデータから、診断、あるいは病態を推定することに相当します。ここでは、3つの個体のデータを作成して、推定を行っています。

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