カテゴリー: Systematic Review

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

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

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

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

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を使ったベイジアンメタアナリシスの環境が整うはずです。何か追加の設定が必要になれば、気軽に聞いてください！

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 の公式ドキュメントを確認するか、追加の設定が必要かどうか調べるのも良いですね。何か困ったことがあれば、いつでも聞いてください！

了解しました。以下のスクリプトでは、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系よりも精度の高い推定が可能です。
さらにカスタマイズしたい点があれば、教えてください！

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のモデルを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のモデルをスムーズに実行できるはずです。
もし追加のエラーが発生した場合は、詳細を教えてください！

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(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のコードは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 をほぼゼロにすることで、変動を抑えられます。

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

その通りです！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(){ } の部分は必須ではありませんが、モデルの収束を安定させるために役立つ場合があります。

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

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

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

了解しました！以下のコードでは、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 distributions → binomial distributions に修正しました。
3. コードの重複部分を削除
   • rt[i] ~ binomial(nt[i], p_c[i] * exp(logRR[i])); の不要な部分を削除しました。

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

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 の値を調整しながら最適化するのがオススメです。何か他に調整したい点があれば聞いてください！

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) や事前分布を調整すると良いでしょう。

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

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

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

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

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() を使って 新規ウインドウ にプロットを出力。

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

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

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

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

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

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

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

その通りです！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) を指定すると、フォレストプロットに 予測区間 が追加されます。

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

その通りです！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) で フォレストプロットに予測区間を追加。

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

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 を使用。

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