R - フラグを使った感度分析

## はじめにデータ分析では、特定の条件に該当する観測値を除外して結果が変わるかどうかを確認する「感度分析（sensitivity analysis）」を行うことがあります。感度分析は、分析結果が特定のサブグループに左右されていないかを検証するための手法です。この記事では、BMI（Body Mass Index）に基づくフラグ変数の作成と、そのフラグを利用した感度分析の手順を紹介します。 ## フラグ変数とはフラグ変数（flag variable）は、特定の条件を満たすかどうかを 0/1 の二値で表す変数です。分析の中で次のような用途に使います。 - 外れ値やハイリスク群の識別 - サブグループ解析の対象選定 - 感度分析での除外基準の明示フラグ変数を作成しておくと、除外条件がコード上で明確になり、再現性も高まります。 ## データの準備とフラグ付与身長と体重の架空データ（30名分）を使います。BMI を計算し、BMI 25 以上の観測値に `flag = 1` を付与します。 ```r df <- data.frame( id = 1:30, height_cm = c(172.3, 158.1, 175.6, 162.4, 168.9, 155.0, 180.2, 160.8, 170.1, 164.5, 176.8, 171.3, 182.5, 157.2, 169.4, 163.0, 178.0, 156.5, 173.1, 166.2, 185.3, 159.7, 167.5, 170.8, 174.2, 161.3, 179.5, 153.8, 171.9, 168.0), weight_kg = c(68.5, 52.3, 80.2, 55.1, 72.0, 48.7, 75.3, 58.9, 85.6, 53.2, 92.4, 49.5, 95.1, 51.0, 70.8, 78.4, 88.2, 54.6, 76.9, 57.3, 98.0, 72.1, 69.4, 83.9, 82.7, 56.8, 78.1, 47.2, 73.5, 63.0) ) df$bmi <- df$weight_kg / (df$height_cm / 100)^2 df$flag <- ifelse(df$bmi >= 25, 1, 0) cat("BMI >= 25 該当者:", sum(df$flag == 1), "名 /", nrow(df), "名中\n") ``` ``` BMI >= 25 該当者: 12 名 / 30 名中 ``` 30名中12名が BMI 25 以上に該当しました。フラグが立った観測値を確認します。 ```r df[df$flag == 1, c("id", "height_cm", "weight_kg", "bmi")] ``` ``` id height_cm weight_kg bmi 1 3 175.6 80.2 26.00910 2 5 168.9 72.0 25.23906 3 9 170.1 85.6 29.58456 4 11 176.8 92.4 29.56021 5 13 182.5 95.1 28.55320 6 16 163.0 78.4 29.50807 7 17 178.0 88.2 27.83739 8 19 173.1 76.9 25.66446 9 21 185.3 98.0 28.54140 10 22 159.7 72.1 28.26998 11 24 170.8 83.9 28.75982 12 25 174.2 82.7 27.25268 ``` ## 散布図でフラグを可視化するフラグの有無を色分けして散布図を描くと、どの観測値が該当しているかを直感的に把握できます。 ```r col <- ifelse(df$flag == 1, "red", "blue") plot(df$height_cm, df$weight_kg, col = col) legend("topleft", c("BMI < 25", "BMI >= 25"), col = c("blue", "red")) ``` ![[fig_scatter_bmi_flag.png]] 赤い点（BMI 25 以上）は、同じ身長の他の観測値と比べて体重が大きい傾向にあります。これらの点が回帰直線の傾きに影響を与えている可能性があるため、除外した場合に結果がどう変わるかを調べます。 ## 感度分析の実施身長から体重を予測する単回帰モデルを、全データとフラグ除外後のデータそれぞれで推定し、結果を比較します。 ### 全データでの回帰 ```r model_all <- lm(weight_kg ~ height_cm, data = df) summary(model_all) ``` ``` Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -173.4851 32.0903 -5.406 9.15e-06 *** height_cm 1.4431 0.1903 7.585 2.91e-08 *** Residual standard error: 8.716 on 28 degrees of freedom Multiple R-squared: 0.6726 ``` 身長が 1 cm 増えると、体重が平均 1.44 kg 増加するという推定結果です。 ### BMI 25 以上を除外した回帰 ```r df_excl <- df[df$flag == 0, ] model_excl <- lm(weight_kg ~ height_cm, data = df_excl) summary(model_excl) ``` ``` Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -114.1150 28.7288 -3.972 0.00109 ** height_cm 1.0542 0.1736 6.073 1.61e-05 *** Residual standard error: 5.653 on 16 degrees of freedom Multiple R-squared: 0.6975 ``` フラグ該当者を除外すると、傾きは 1.05 に減少しました。残差の標準誤差も 8.72 から 5.65 に改善されています。 ### 結果の比較 | 指標 | 全データ | BMI 25 以上を除外 | |---|---|---| | 回帰係数（身長） | 1.443 | 1.054 | | 標準誤差 | 0.190 | 0.174 | | $R^2$ | 0.673 | 0.698 | | 残差標準誤差 | 8.716 | 5.653 | | p値 | 2.91e-08 | 1.61e-05 | 回帰係数は全データで 1.44、除外後は 1.05 と、約 0.4 の差が生じました。一方で、いずれも統計的に有意（$p < 0.001$）であり、身長と体重の正の関連という結論自体は変わりません。 ## 回帰直線の比較 2つの回帰直線を同じ散布図に重ねると、傾きの違いが視覚的に確認できます。 ```r plot(df$height_cm, df$weight_kg, col = col) abline(model_all) abline(model_excl, lty = 2) legend("topleft", c("All data", "BMI<25 only"), lty = c(1, 2)) ``` ![[fig_regression_bmi_flag.png]] 実線（全データ）の方が破線（BMI 25 未満のみ）よりも傾きが急になっています。BMI 25 以上の観測値が、身長と体重の関連を実際より強く見せる方向に働いていることがわかります。 ## 感度分析の報告感度分析の結果は、メインの解析結果と並べて報告します。論文では次のような記述が一般的です。 > 感度分析として、BMI 25 以上の 12 名を除外して同様の回帰分析を行った。身長の回帰係数は 1.05（SE = 0.17, p < 0.001）であり、全データでの推定値 1.44 と比較して小さくなったが、身長と体重の正の関連は維持された。感度分析で結論が変わらなければ、結果の頑健性（robustness）を主張する根拠になります。逆に結論が変わる場合は、その要因を考察する必要があります。