非線形回帰Rプログラミング: confint.profile.nlsで信頼区間を計算する方法
Rプログラミングにおけるconfint.profile.nls
についてですね。これは、非線形最小二乗法 (nls
関数) でフィットさせたモデルのパラメーターに対する信頼区間を、プロファイル尤度 (profile likelihood) の方法を用いて計算するための機能です。
以下で詳しく説明します。
nls
とは?
まず、nls
(Nonlinear Least Squares) はRの標準パッケージ stats
に含まれる関数で、非線形回帰モデルのパラメーターを推定するために使われます。線形回帰とは異なり、応答変数と予測変数の関係が直線的ではない場合に適用されます。
信頼区間とは?
推定されたパラメーターの値がどの程度の範囲に収まる可能性が高いかを示すものです。例えば95%信頼区間であれば、「この区間に真のパラメーター値が含まれる確率が95%である」と解釈されます。
なぜ confint.profile.nls
が必要か?
一般的な線形回帰モデルでは、パラメーターの信頼区間は正規分布の仮定に基づいて比較的容易に計算できます(いわゆる「Wald信頼区間」)。しかし、非線形回帰モデルの場合、パラメーターの推定値の分布は必ずしも正規分布に従いません。特にサンプルサイズが小さい場合やモデルが複雑な場合、Wald信頼区間は信頼性が低いことがあります。
そこで、より信頼性の高い信頼区間を計算するために「プロファイル尤度法」が用いられます。confint.profile.nls
は、nls
オブジェクト(nls
関数でフィットされたモデルの結果)に対してこのプロファイル尤度法を適用するためのS3メソッドです。
プロファイル尤度法は、ある特定のパラメーターの信頼区間を求める際に、そのパラメーターを固定し、残りのパラメーターを最適化するというプロセスを繰り返すことで、尤度関数(データがモデルから生成される確率を表す関数)がどのように変化するかを追跡します。
具体的には、あるパラメーター θiの信頼区間を計算する場合、次のような手順をとります。
- θi以外のすべてのパラメーターを最適化し、そのときの尤度を計算します。
- θiの値を少しずつ変えながら、1.のプロセスを繰り返します。
- 尤度関数が特定のしきい値(通常はカイ二乗分布に基づいて決定される)を下回る θiの値の範囲を特定します。この範囲が信頼区間となります。
この方法により、非線形モデルにおけるパラメーターの非対称な分布や相関を考慮した、より正確な信頼区間が得られます。
confint.profile.nls
の使い方(概念)
通常、nls
モデルをフィットした後、profile()
関数を使ってプロファイルオブジェクトを作成し、そのプロファイルオブジェクトに対してconfint()
を呼び出します。
例:
# 例として、Michaelis-Menten モデルをフィットする
# データはRの組み込みデータセット `BOD` を使用
# 非線形回帰モデルのフィット
fm1 <- nls(demand ~ SSasympOrig(Time, A, lrc), data = BOD)
# プロファイルオブジェクトの作成
# これが 'profile.nls' メソッドを呼び出す
pr1 <- profile(fm1)
# 信頼区間の計算
# これが 'confint.profile.nls' メソッドを呼び出す
confint(pr1)
# または、直接 nls オブジェクトに confint を適用することも可能
# この場合、内部的に profile() が呼ばれる
confint(fm1)
confint.profile.nls
の利点
- バイアス補正
小さなサンプルサイズにおける推定量のバイアスを部分的に補正する効果も期待できます。 - 非対称な信頼区間
パラメーターの信頼区間が非対称になる場合でも、その非対称性を適切に捉えます。これはWald信頼区間では不可能です。 - より正確な信頼区間
非線形モデルにおいて、パラメーターの分布が非正規である場合でも、より正確な信頼区間を提供します。
- nlstools パッケージ
Rの標準的なconfint
関数もnls
オブジェクトに対してプロファイル尤度法を使用しますが、nlstools
パッケージのconfint2
関数は、漸近法(asymptotic method)とプロファイル法(profile method)の両方を選択できるなど、より柔軟なオプションを提供することもあります。 - 収束
nls
モデルのフィット自体が収束しない場合や、プロファイル計算が収束しない場合もあります。その際は、nls
の初期値や制御パラメーターの調整が必要になることがあります。 - 計算コスト
プロファイル尤度法は、各パラメーターの信頼区間を計算するために多くのモデルフィットを内部的に実行するため、計算に時間がかかる場合があります。特にモデルが複雑であったり、データ量が多い場合は顕著です。
confint.profile.nls
の一般的なエラーとトラブルシューティング
エラー: Error in nlsModel(formula, mf, rval, class(data), ...) または Error in nls(formula, data, start, control, algorithm)
これはconfint.profile.nls
が呼び出される前のnls
モデルのフィット自体に問題があることを示しています。
- トラブルシューティング
- 初期値の調整
- データのプロットを確認し、関数形とデータが合うように大まかなパラメーターの値を手動で推測します。
nls
のstart
引数にこれらの初期値を指定します。nls2
パッケージのnls_multstart
やnlstools
パッケージのnls_boot
など、複数の初期値を試すことができる関数も検討してください。- 線形化できる場合は、まず線形回帰で大まかな初期値を得ることも有効です。
- モデル式の確認
モデル式が正しくRの構文で記述されているか、論理的に正しいかを確認します。 - データの確認
欠損値の処理、外れ値の有無、データの分布などを確認します。 - パラメーターの識別性
モデルが複雑すぎる場合、パラメーターの数を減らすことを検討したり、理論的にパラメーターが独立に識別できるかを再考します。
- 初期値の調整
- 考えられる原因
- 不適切な初期値 (start)
nls
は反復的なアルゴリズムでパラメーターを推定するため、良い初期値が非常に重要です。初期値が真の値から離れすぎていると、収束しないか、局所的な最適解に陥ることがあります。 - モデル式の誤り
数式がデータに合っていない、またはRの構文として間違っている可能性があります。 - データの問題
NA
値が含まれている、データが特定のモデルに合わない(例:定数データ、範囲が狭すぎる)、外れ値が多すぎるなど。 - 識別性の問題
モデルのパラメーターがデータから一意に識別できない場合。例えば、異なるパラメーターの組み合わせが同じ出力をもたらす場合など。
- 不適切な初期値 (start)
エラー: Error in profile.nls(fm1) : profiling detected new, lower deviance
これはプロファイリングの過程で、最初にnls
が報告した最小二乗和(または逸脱度)よりも小さい値が見つかったことを意味します。つまり、nls
の初期の最適化が「真の」グローバルな最小値に到達していなかった可能性を示唆しています。
- トラブルシューティング
- nlsのcontrol引数の調整
control = nls.control(maxiter = 1000, tol = 1e-06, minFactor = 1/2048, warnOnly = TRUE)
のように、maxiter
(最大反復回数) を増やしたり、tol
(許容誤差) を小さくしたりして、nls
の収束条件を厳しくしてみます。warnOnly = TRUE
にすると、エラーではなく警告として扱われる場合がありますが、根本的な解決にはなりません。
- nlsのアルゴリズム変更
デフォルトのアルゴリズム ("port"
や"plinear"
) が適切でない場合、他のアルゴリズムを試すこともできます。ただし、これは稀なケースです。 - 初期値の再確認
nls
の初期値が最適解から遠すぎると、収束が不十分になることがあります。
- nlsのcontrol引数の調整
- 考えられる原因
nls
の収束が不十分だった。- 尤度プロファイル計算中に、より良い(小さい)最小二乗和の点が見つかった。
警告メッセージ: Warning: step size became too small, not enough points for profile
これは、プロファイル尤度曲線を構築するために必要なデータポイントが十分に得られなかったことを示唆しています。信頼区間の計算が不正確になる可能性があります。
- トラブルシューティング
- profile関数のstepsizeやmaxsteps引数の調整
pr1 <- profile(fm1, stepsize = 0.01)
のようにstepsize
を小さくして、より細かく尤度プロファイル点を計算するように指示します。pr1 <- profile(fm1, maxsteps = 200)
のようにmaxsteps
を増やし、プロファイル計算の最大ステップ数を増やします。ただし、計算時間が増加します。
- モデルの再評価
もしかしたらモデルがデータに対して過剰に複雑すぎるか、パラメーターの一部が実質的に識別不能である可能性があります。モデルの簡素化を検討します。 - データ量の確認
非常にデータ量が少ない場合、安定したプロファイル尤度を得るのが難しいことがあります。
- profile関数のstepsizeやmaxsteps引数の調整
- 考えられる原因
- モデルのパラメーター空間が非常に平坦である(尤度関数が平坦で、わずかなパラメーターの変化で尤度があまり変わらない)。
- プロファイリングアルゴリズムが、信頼区間の端に到達する前に停止してしまった。
計算時間が非常に長い、またはフリーズする
プロファイル尤度法は、各パラメーターについて最適化を繰り返すため、計算負荷が高いです。
- トラブルシューティング
- 忍耐
まずは待ってみる。特に複雑なモデルや大規模なデータでは時間がかかります。 - プロファイルするパラメーターの選択
全てのパラメーターではなく、特定の関心のあるパラメーターのみをプロファイルする (profile(fm1, parm = "parameter_name")
) ことで計算時間を短縮できます。 - 計算リソースの確認
メモリやCPUの使用状況を確認し、リソースが不足していないか確認します。 - nls.controlの調整
nls
の収束を速くするために、control
引数を調整することを検討します。 - ブートストラップ法
プロファイル尤度法がうまくいかない、または時間がかかりすぎる場合、ブートストラップ法による信頼区間の計算も有効な代替手段です。nlstools
パッケージのnlsBoot
などが利用できます。
- 忍耐
- 考えられる原因
- モデルが非常に複雑(パラメーター数が多い)。
- データセットが大きい。
- プロファイルのステップサイズが小さすぎるか、最大ステップ数が大きすぎる。
信頼区間が非現実的な値になる(非常に広い、負の値になるなど)
計算は完了しても、得られた信頼区間が統計的に意味をなさない場合。
- トラブルシューティング
- モデル診断
plot(fm1)
を使って残差プロットなどを確認し、モデルがデータに適切にフィットしているか診断します。モデルがデータに合っていない場合、信頼区間も意味をなしません。 - プロファイルプロットの確認
plot(profile(fm1))
を実行して、各パラメーターのプロファイル尤度曲線がどのように見えるかを確認します。曲線が非常に平坦であったり、非対称性が強すぎたりする場合、パラメーターの識別性や信頼区間の計算に問題がある可能性があります。 - パラメーター制約の導入
もしパラメーターに物理的な制約(例:パラメーターが常に正であるべき)がある場合、nls
のlower
やupper
引数(一部のnls
実装や、nls.lm
などの代替関数で使用可能)や、モデル式の中で制約を組み込む(例:パラメーターをexp()
で囲むなど)ことを検討します。 - データ量の増加
信頼区間の精度はデータ量に依存します。データが少なすぎる場合、信頼区間は広くなる傾向があります。
- モデル診断
- 考えられる原因
- モデルの適合が悪い。
- モデルのパラメーターがデータからうまく識別できない。
- プロファイル尤度曲線が非常に平坦、または歪んでいる。
- パラメーターに物理的な制約があるが、モデルにその制約が組み込まれていない。
- 統計的アドバイス
解決が困難な場合は、統計専門家や経験豊富なRユーザーに相談することも良いでしょう。 - 代替パッケージの検討
nls
がうまくいかない場合、minpack.lm
パッケージのnls.lm
関数(Levenberg-Marquardtアルゴリズム)や、nlme
パッケージのgnls
関数(より柔軟なエラー構造に対応)など、他の非線形回帰関数を試すことも有効です。これらの関数は、confint
メソッドも提供している場合があります。 - summary(nls_object)
nls
オブジェクトのsummary
は、推定されたパラメーター、標準誤差、残差の自由度などを提供します。これで基本的なフィットの状況を確認できます。 - 初期値の試行錯誤
初期値はnls
成功の鍵です。グラフを見て、データがどのような形状をしているか把握し、それに基づいて初期値を設定します。 - 単純なモデルから始める
複雑なモデルを試す前に、より単純なモデルでデータにフィットすることを試み、基本的な理解を深めます。
基本的な使用例:Michaelis-Menten モデル
酵素反応速度論でよく使われるMichaelis-Mentenモデルを例に見てみましょう。
rate=Km​+concVmax​⋅conc​
ここで、Vmax​は最大反応速度、Km​はMichaelis定数です。
# データの準備(例として架空のデータを作成)
set.seed(123)
conc <- c(0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0)
Vm_true <- 150 # 真のVm値
Km_true <- 2.0 # 真のKm値
# ノイズを加えた反応速度
rate <- (Vm_true * conc) / (Km_true + conc) + rnorm(length(conc), mean = 0, sd = 5)
# データをデータフレームにまとめる
df_mm <- data.frame(conc = conc, rate = rate)
# データのプロット (初期値の目安を立てるために重要)
plot(df_mm$conc, df_mm$rate,
xlab = "Concentration",
ylab = "Reaction Rate",
main = "Michaelis-Menten Data",
pch = 16, col = "blue")
# nls を用いてモデルをフィット
# 初期値 (start) はデータのプロットから大まかに推測する
# Vm は最大のrateのあたり、Km は最大のrateの半分になるconcのあたり
fit_mm <- nls(rate ~ Vm * conc / (Km + conc),
data = df_mm,
start = list(Vm = 140, Km = 1.5))
# フィット結果の概要
summary(fit_mm)
# フィットされたモデルをプロットに追加
curve((coef(fit_mm)["Vm"] * x) / (coef(fit_mm)["Km"] + x),
add = TRUE, col = "red", lwd = 2)
legend("bottomright", legend = c("Data", "Fitted Curve"),
col = c("blue", "red"), pch = c(16, NA), lty = c(NA, 1), lwd = c(NA, 2))
# プロファイルオブジェクトの作成
# これにより、各パラメーターの尤度プロファイルが計算される
# このステップは時間がかかる場合があります
cat("\n--- Profiling Michaelis-Menten model ---\n")
profile_mm <- profile(fit_mm)
# プロファイルプロットの表示 (オプション)
# これにより、各パラメーターの尤度プロファイルの形状を確認できる
# 曲線が対称的であれば、正規近似が妥当な可能性が高い
plot(profile_mm)
# confint.profile.nls を用いて信頼区間を計算
# profile オブジェクトを直接 confint に渡すのが一般的
cat("\n--- Confidence Intervals for Michaelis-Menten model ---\n")
ci_mm <- confint(profile_mm, level = 0.95) # 95%信頼区間
print(ci_mm)
# または、profile()を明示的に呼び出さずに、nlsオブジェクトに直接confintを適用することも可能
# この場合、内部的にprofile()が実行される
cat("\n--- Confidence Intervals (direct call) ---\n")
ci_mm_direct <- confint(fit_mm, level = 0.95)
print(ci_mm_direct)
別の例:ロジスティック成長モデル
生物の個体群成長や、反応の飽和現象などで用いられるロジスティック成長モデルもよく非線形回帰で扱われます。
y=1+exp((xmid−x)/scal)Asym​
ここで、Asymは漸近線(最大値)、xmidは変曲点でのx値、scalはスケールパラメーターです。
# データの準備(例として架空のデータを作成)
set.seed(456)
x_val <- seq(0, 30, by = 1)
Asym_true <- 100 # 真の漸近線
xmid_true <- 15 # 真の変曲点
scal_true <- 2 # 真のスケール
y_true <- Asym_true / (1 + exp((xmid_true - x_val) / scal_true))
y_obs <- y_true + rnorm(length(x_val), mean = 0, sd = 3) # ノイズを加える
df_logistic <- data.frame(x = x_val, y = y_obs)
# データのプロット
plot(df_logistic$x, df_logistic$y,
xlab = "Time",
ylab = "Population/Value",
main = "Logistic Growth Data",
pch = 16, col = "darkgreen")
# nls を用いてモデルをフィット
# 初期値はデータのプロットから推測
# Asym: 最大のy値のあたり
# xmid: yがAsym/2になるx値のあたり
# scal: 傾きに関係、xmidの前後でyが大きく変わる範囲の幅の目安
fit_logistic <- nls(y ~ Asym / (1 + exp((xmid - x) / scal)),
data = df_logistic,
start = list(Asym = 105, xmid = 14, scal = 2.5))
# フィット結果の概要
summary(fit_logistic)
# フィットされたモデルをプロットに追加
curve(coef(fit_logistic)["Asym"] / (1 + exp((coef(fit_logistic)["xmid"] - x) / coef(fit_logistic)["scal"])),
add = TRUE, col = "purple", lwd = 2)
legend("bottomright", legend = c("Data", "Fitted Curve"),
col = c("darkgreen", "purple"), pch = c(16, NA), lty = c(NA, 1), lwd = c(NA, 2))
# プロファイルオブジェクトの作成
cat("\n--- Profiling Logistic Growth model ---\n")
profile_logistic <- profile(fit_logistic)
# プロファイルプロットの表示 (オプション)
plot(profile_logistic)
# confint.profile.nls を用いて信頼区間を計算
cat("\n--- Confidence Intervals for Logistic Growth model ---\n")
ci_logistic <- confint(profile_logistic, level = 0.95)
print(ci_logistic)
特定のパラメーターの信頼区間のみを計算
confint
関数にはparm
引数があり、これにより特定のパラメーターの信頼区間のみを計算することができます。計算時間を短縮したい場合に便利です。
# Michaelis-Menten モデルのKmのみの信頼区間
cat("\n--- Confidence Interval for Km only ---\n")
ci_km_only <- confint(profile_mm, parm = "Km", level = 0.95)
print(ci_km_only)
# ロジスティック成長モデルのAsymとxmidの信頼区間
cat("\n--- Confidence Intervals for Asym and xmid only ---\n")
ci_asym_xmid_only <- confint(profile_logistic, parm = c("Asym", "xmid"), level = 0.95)
print(ci_asym_xmid_only)
profile()
関数には、プロファイリングの挙動を制御するための引数があります。エラーや警告が出る場合に調整すると良いでしょう。
maxsteps
: 各パラメーターについて計算する最大ステップ数。警告が出る場合、これを増やすと解決することがあります。stepsize
: プロファイリングのステップサイズ。小さいほど詳細に、大きいほど大まかに計算されます。デフォルトでは自動的に調整されます。
# より細かいステップでプロファイリングを試す例
# (通常はデフォルトで十分ですが、収束に問題がある場合に試す)
cat("\n--- Profiling with smaller stepsize (example) ---\n")
profile_mm_fine <- profile(fit_mm, stepsize = 0.005) # デフォルトよりも小さく
ci_mm_fine <- confint(profile_mm_fine)
print(ci_mm_fine)
# 最大ステップ数を増やす例
cat("\n--- Profiling with more steps (example) ---\n")
profile_mm_more_steps <- profile(fit_mm, maxsteps = 100) # デフォルトはもっと少ない
ci_mm_more_steps <- confint(profile_mm_more_steps)
print(ci_mm_more_steps)
confint.profile.nls
の代替方法
主な代替方法としては、以下のものが挙げられます。
- Wald信頼区間 (Asymptotic/Normal Approximation)
- ブートストラップ法 (Bootstrap Method)
- デルタ法 (Delta Method)
- ベイズ法 (Bayesian Method)
それぞれについて詳しく見ていきましょう。
Wald信頼区間 (Asymptotic/Normal Approximation)
これは最もシンプルで、計算も速い方法です。nls
オブジェクトのsummary()
から得られる標準誤差と、漸近正規性の仮定に基づいて信頼区間を計算します。
-
Rでの実装例
confint.profile.nls
のように直接的な関数はありませんが、summary(nls_object)
の結果から手動で計算できます。# Michaelis-Menten モデルの例(再掲) set.seed(123) conc <- c(0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0) Vm_true <- 150 Km_true <- 2.0 rate <- (Vm_true * conc) / (Km_true + conc) + rnorm(length(conc), mean = 0, sd = 5) df_mm <- data.frame(conc = conc, rate = rate) fit_mm <- nls(rate ~ Vm * conc / (Km + conc), data = df_mm, start = list(Vm = 140, Km = 1.5)) # summaryから標準誤差と自由度を取得 s_fit <- summary(fit_mm) coefs <- s_fit$coefficients # 推定値、標準誤差、t値、p値を含む行列 df_res <- s_fit$df[2] # 残差の自由度 # Wald信頼区間の計算(95%信頼水準) alpha <- 0.05 t_crit <- qt(1 - alpha/2, df = df_res) wald_ci <- matrix(NA, nrow = nrow(coefs), ncol = 2, dimnames = list(rownames(coefs), c("2.5 %", "97.5 %"))) for (i in 1:nrow(coefs)) { est <- coefs[i, "Estimate"] se <- coefs[i, "Std. Error"] wald_ci[i, 1] <- est - t_crit * se wald_ci[i, 2] <- est + t_crit * se } cat("--- Wald Confidence Intervals ---\n") print(wald_ci) # `nlstools`パッケージの`confint2`関数は、Wald信頼区間も計算できます # install.packages("nlstools") # 必要であればインストール library(nlstools) cat("\n--- Wald CI using nlstools::confint2(method = 'asymptotic') ---\n") print(confint2(fit_mm, method = "asymptotic"))
-
欠点
- 非線形モデルでは、パラメーターの分布が正規分布から大きく外れる場合があるため、信頼性が低いことがあります。特にサンプルサイズが小さい場合や、モデルの非線形性が強い場合に問題となります。
- 信頼区間が非対称になることを考慮できない(常に推定値を中心とした対称な区間になる)。
- 制約のあるパラメーター(例:負の値をとらないはずのパラメーターが負の信頼区間を持つ)の場合に不適切。
-
利点
- 計算が非常に速い。
- 実装が簡単。
-
- 非線形最小二乗法のパラメーター推定量は、十分なサンプルサイズがあれば漸近的に正規分布に従うという性質を利用します。
- 信頼区間は「推定値 ± t値 × 標準誤差」で計算されます。
ブートストラップ法 (Bootstrap Method)
ブートストラップは、観測されたデータから多数の再サンプリングを行い、それぞれの再サンプリングデータセットでモデルをフィットすることで、パラメーターの推定値の経験的分布を構築する方法です。
-
欠点
- 計算負荷が非常に高い(数千回の
nls
フィットが必要)。 nls
のフィットが不安定な場合(収束しないなど)は、ブートストラップの繰り返しがうまくいかないことがある。- 再サンプリングの方法(残差ブートストラップ、ケースブートストラップなど)の選択が結果に影響を与える場合がある。
- 計算負荷が非常に高い(数千回の
-
利点
- 特定の分布(正規分布など)を仮定しないため、非線形モデルに特に適している。
- 信頼区間が非対称になることを自然に考慮できる。
- モデルの仮定からの逸脱に対してロバスト。
-
考え方
- 元のデータから重複を許してN個のサンプルをランダムに抽出する(リサンプリング)。
- このリサンプリングされたデータセットで
nls
モデルをフィットし、パラメーターを推定する。 - このプロセスを数千回繰り返し、パラメーターの推定値の分布を作成する。
- この分布のパーセンタイル(例:2.5パーセンタイルと97.5パーセンタイル)を信頼区間とする。
デルタ法 (Delta Method)
デルタ法は、非線形関数の分散を、一次のテイラー展開近似を用いて線形近似することで推定する方法です。これは、パラメーターの信頼区間だけでなく、パラメーターの非線形関数の信頼区間(例:Vm/Kmなどの比率)を計算する際にもよく用いられます。
-
Rでの実装例
デルタ法を直接実装する関数はRのbaseにはありませんが、car
パッケージのdeltaMethod
関数などが利用できます。# install.packages("car") # 必要であればインストール library(car) # Michaelis-Menten モデルのフィット(再掲) # ... (上記のfit_mmを使用) cat("\n--- Delta Method Confidence Intervals ---\n") # 各パラメーターの信頼区間 # 単純なパラメーターの場合、Wald法とほぼ同じ結果になることが多い print(deltaMethod(fit_mm, "Vm", level = 0.95)) print(deltaMethod(fit_mm, "Km", level = 0.95)) # パラメーターの非線形関数の信頼区間(例:VmとKmの比率) # Michaelis-Mentenモデルでは意味のある比率ではないが、例として cat("\n--- Delta Method for a non-linear function (e.g., Vm/Km) ---\n") print(deltaMethod(fit_mm, "Vm/Km", level = 0.95))
propagate
パッケージのpredictNLS
関数もデルタ法に基づいて予測区間や信頼区間を計算できます。 -
欠点
- 一次のテイラー展開に基づくため、モデルの非線形性が強い場合やサンプルサイズが小さい場合に精度が低下する可能性がある。
- Wald法と同様に、正規性の仮定に依存する。
- 信頼区間が非対称になることを考慮できない。
-
利点
- 計算が比較的速い。
- パラメーターの非線形関数の信頼区間も計算できる。
- 漸近的に正確(サンプルサイズが十分大きい場合)。
-
考え方
- 推定されたパラメーターの共分散行列(
vcov(nls_object)
で取得)を使用します。 - 非線形関数をパラメーターの点で一次のテイラー展開で近似し、線形近似された関数の分散を導出します。
- 推定されたパラメーターの共分散行列(
ベイズ統計学では、パラメーターを確率変数とみなし、データと事前情報(事前分布)に基づいてパラメーターの事後分布を導出します。この事後分布から、パラメーターの信頼区間(信用区間または確信区間と呼ばれる)を直接推定することができます。
-
欠点
- モデルの構築とMCMCサンプリングに専門的な知識と時間が必要。
- MCMCの収束を確認する必要がある。
- 計算負荷が非常に高い。
-
利点
- 信頼区間が非対称になることを自然に考慮できる。
- 事前情報(例:パラメーターは正であるべき、など)をモデルに組み込むことができる。
- 小サンプルサイズでも比較的安定した結果が得られやすい。
- パラメーターの不確実性をより直感的に解釈できる。
-
考え方
- パラメーターの事前分布を設定する。
- 尤度関数(データが与えられたときのパラメーターの確率)と事前分布を掛け合わせ、事後分布を計算する(通常、MCMC法などのサンプリング手法を用いる)。
- 事後分布のサンプリング結果から、信頼区間(例:事後分布の2.5パーセンタイルと97.5パーセンタイル)を直接求める。
方法 | 利点 | 欠点 | 用途 |
---|---|---|---|
Wald信頼区間 | 計算が速く、実装が簡単 | 非線形モデルでは信頼性が低い可能性、非対称な区間を扱えない | 簡易的な確認、線形性が高いと期待される場合 |
ブートストラップ法 | 分布の仮定なし、非対称区間を扱える、ロバスト | 計算負荷が高い、収束問題の可能性 | 信頼性の高い区間が必要な場合、複雑なモデル |
デルタ法 | 計算が比較的速い、非線形関数の区間も扱える | 非線形性が強いと精度低下、正規性の仮定、非対称区間を扱えない | パラメーター関数の区間計算、十分なサンプルサイズの場合 |
ベイズ法 | 事前情報を組み込める、非対称区間を扱える、小サンプルにも対応 | 専門知識が必要、計算負荷が高い、MCMCの収束確認が必要 | 厳密な不確実性評価、複雑なモデル、事前情報がある場合 |