【Rで医療統計】ロジスティック回帰とは：共変量が1つか２つの例

ロジスティック回帰とは？

　ロジスティック回帰とは、ある出来事が起こる確率を予測するための統計モデル です。たとえば、次のような予測をしたいとします：

ある人が肺がんになるかどうか（「なる」か「ならない」）
送られてきたメールがスパムかどうか（「そう」か「そうじゃない」）
ある商品が売れるかどうか（「売れる」か「売れない」）

　このように 結果が「Yesか、Noか」のような分類問題を考えるときに、ロジスティック回帰が使われます。

　「Yes」と「No」のような言葉のままでは、数式が使いづらいですし、はっきりと白黒つけることが難しい場合もありますので、ロジスティック回帰では、Yesである確率を予測するモデルを構築します。

　普通の線形回帰（直線の関係を求める回帰）は、予測値が0以下になったり、1を超えたりすることがあります。でも、確率は 必ず0から1の間 なので、ロジスティック回帰はシグモイド関数という曲線を使って確率を表現 します。

ロジスティック回帰とは？
ロジスティック回帰の実際 (架空の喫煙者の分析)
- 共変量が1つの例（喫煙期間が長いと、肺がんになりやすい？）
  - Rでロジスティック回帰
  - Rスクリプトの説明
- 共変量が2つの例（喫煙期間だけでなく、年齢も肺がんになりやすさと関係する？）
  - Rでロジスティック回帰
  - Rスクリプトの説明
まとめ

ロジスティック回帰の実際 (架空の喫煙者の分析)

　以下では、架空の話として「喫煙期間が肺がんの発症リスクに与える影響を評価したい」と考えたとします。なお、本分析で使用するデータは、そもそもの仮定として、喫煙期間が長いほど肺がんの発症確率が高くなるように作成されており、必ずしも実際のデータを反映しているわけではありません。

共変量が1つの例（喫煙期間が長いと、肺がんになりやすい？）

　ということで、「喫煙期間が長いと、肺がんになる確率はどれくらい上がるのか？」をロジスティック回帰で調べてみます。

　説明変数（共変量、独立変数とも言います）は、喫煙期間（年数）を表す $\displaystyle{ x }$ です。

　目的変数（従属変数）は $\displaystyle{ y }$ で、肺がんになったかどうか (1: なった, 0: ならなかった)を表します。 $\displaystyle{ y }$ は確率ではなく、0か1の値をとります。

　ロジスティック回帰の式は次のようになります：

$\displaystyle{ P(y=1 | x) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x)}} }$

この式では、

$\displaystyle{ P(y=1 | x) }$ は「喫煙年数 $\displaystyle{ x }$ の人が肺がんになる確率」
$\displaystyle{ \beta _ 0 }$ は「切片（intercept）」で、基本的な確率の基準
$\displaystyle{ \beta _ 1 }$ は「回帰係数（coefficient）」で、喫煙年数が1年増えると確率がどのくらい変わるかを表す
$\displaystyle{ e }$ は自然対数の底

Rでロジスティック回帰

　ここでは、Rを使って架空のデータを生成し、喫煙年数と肺がんの関係を調べてみます。

# データ数
n <- 100

# 喫煙期間をランダムに生成（0〜40年）
smoking_years <- runif(n, 0, 40)

# 肺がんになる確率をロジスティック関数で計算
prob <- 1 / (1 + exp(-(-4 + 0.2 * smoking_years)))

# 確率に基づいて肺がんになったかどうかを決定（0 or 1）
lung_cancer <- rbinom(n, 1, prob)

# ロジスティック回帰モデルの構築
model1 <- glm(lung_cancer ~ smoking_years, family = binomial)

# 結果の表示
print(summary(model1))

# 予測用データを作成
new_smoking_years <- seq(0, 40, length.out = 100)
new_prob <- predict(model1, newdata = data.frame(smoking_years = new_smoking_years), type = "response")

# プロット
plot(smoking_years, lung_cancer, main = "喫煙期間と肺がんの関係", las=1,
     xlab = "喫煙期間（年）", ylab = "肺がんの発生確率", pch = 19, col = 4-lung_cancer*2)
lines(new_smoking_years, new_prob, col = "blue", lwd = 2)

Rスクリプトの説明

　上のRスクリプトは、喫煙期間（年数）と肺がんの発生確率の関係をロジスティック回帰で分析し、その結果をプロットするものです。以下で、各コマンドと変数につい説明します。

1. データの準備

# データ数
n <- 100

これはデータ数（サンプル数）を指定しています。
n = 100 としているので、100人分のデータ をシミュレーションします。

# 喫煙期間をランダムに生成（0〜40年）
smoking_years <- runif(n, 0, 40)

喫煙期間（年数）を0〜40の範囲でランダムに生成しています。
runif(n, 0, 40) は、一様分布（Uniform Distribution） に従う乱数を生成する関数です。

2. 肺がん発生確率の計算

# 肺がんになる確率をロジスティック関数で計算
prob <- 1 / (1 + exp(-(-4 + 0.2 * smoking_years)))

ロジスティック関数（シグモイド関数）を使用して、肺がんになる確率を計算 しています。

$\displaystyle{ P(y=1 | x) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x)}} }$

ここでは、

切片 $\displaystyle{ \beta _ 0 }$ : -4
喫煙年数の影響 $\displaystyle{ \beta _ 1 }$ : 0.2
説明変数 $x$ : smoking_years（喫煙年数）
exp() で指数関数を計算

3. 確率に基づく肺がん発生データの作成

# 確率に基づいて肺がんになったかどうかを決定（0 or 1）
lung_cancer <- rbinom(n, 1, prob)

rbinom(n, 1, prob) は、二項分布に基づいて0または1をランダムに生成 する関数です。

4. ロジスティック回帰モデルの構築

# ロジスティック回帰モデルの構築
model1 <- glm(lung_cancer ~ smoking_years, family = binomial)

glm() は一般化線形モデル（Generalized Linear Model） を作成する関数です。

lung_cancer ~ smoking_yearsで、目的変数（従属変数）をlung_cancer（肺がんの有無）、説明変数（独立変数）をsmoking_years（喫煙期間）に指定
family = binomial: ロジスティック回帰を指定（目的変数が0または1のカテゴリデータであるため）

5. 結果の表示

# 結果の表示
print(summary(model1))

ロジスティック回帰の結果を表示しています。

出力例： - 切片（Intercept）: $\displaystyle{ \beta _ 0 }$ の値 - 喫煙期間の係数（Estimate）: $\displaystyle{ \beta _ 1 }$ の値（喫煙期間が1年増えるごとに肺がんの確率がどれくらい変わるか）

6. 予測データの作成

# 予測用データを作成
new_smoking_years <- seq(0, 40, length.out = 100)

0 から 40 までの値を 100個 均等に作成（滑らかな曲線を描くため）

new_prob <- predict(model1, newdata = data.frame(smoking_years = new_smoking_years), type = "response")

model1 を使って new_smoking_years の肺がん確率を予測。
type = "response" を指定することで、確率として出力。

7. グラフの描画

# プロット
plot(smoking_years, lung_cancer, main = "喫煙期間と肺がんの関係", las=1,
     xlab = "喫煙期間（年）", ylab = "肺がんの発生確率", pch = 19, col = 4-lung_cancer*2)
lines(new_smoking_years, new_prob, col = "blue", lwd = 2)

plotで元のデータのばらつきを描いています。lineでフィットした曲線を重ね書きしています。

　この架空のデータでは、「喫煙期間が長いと肺がんのリスクが上がる」 ことになっていますが、現実の傾向については、論文などを参照してください。

共変量が2つの例（喫煙期間だけでなく、年齢も肺がんになりやすさと関係する？）

　次に、「喫煙期間だけでなく、年齢も考慮したら肺がんの確率はどうなるか？」をロジスティック回帰で分析します。ここでも、データはあくまで架空のものですので、信じないでください。

　ここでは、以下の２つ説明変数 (共変量)を考えます。

$x _ 1$ ：喫煙期間（年数）
$x _ 2$ ：年齢（歳）

　目的変数は、肺がんになったかどうか (1: なった, 0: ならなかった)を表す， $y$ です。

このとき、ロジスティック回帰の式は次のようになります：

$\displaystyle{ P(y=1 | x_1, x_2) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1 + \beta_2 x_2)}} }$

$\displaystyle{ \beta _ 1 }$ は喫煙年数の影響
$\displaystyle{ \beta _ 2 }$ は年齢の影響

Rでロジスティック回帰

　Rを使って、架空のデータを生成し、喫煙年数、年齢と肺がんの関係を調べてみます。

# データ数
n <- 100

# 年齢をランダムに生成（30〜80歳）
age <- runif(n, 30, 80)

# 喫煙期間（年数）をランダムに生成（0〜40年）
smoking_years <- runif(n, 0, 40)
# 喫煙期間を年齢を考慮して修正
smoking_years[smoking_years > age - 20] <- smoking_years[smoking_years > age - 20] - 20

# 肺がんになる確率をロジスティック関数で計算（喫煙年数と年齢を考慮）
prob <- 1 / (1 + exp(-(-6 + 0.1 * smoking_years + 0.08 * age)))

# 確率に基づいて肺がんになったかどうかを決定（0 or 1）
lung_cancer <- rbinom(n, 1, prob)

# ロジスティック回帰モデルの構築（共変量2つ）
model2 <- glm(lung_cancer ~ smoking_years + age, family = binomial)

# グリッド状のデータを作成（喫煙期間と年齢の範囲）
smoking_grid <- seq(0, 40, length.out = 50)  # 喫煙期間のグリッド（0〜40年を50ステップ）
age_grid <- seq(30, 80, length.out = 50)  # 年齢のグリッド（30〜80歳を50ステップ）

# ヒートマップ用の確率行列を作成
prob_matrix <- matrix(0, nrow = length(smoking_grid), ncol = length(age_grid))

# 各グリッド点での肺がん発生確率を計算
for (i in 1:length(smoking_grid)) {
  for (j in 1:length(age_grid)) {
    prob_matrix[i, j] <- predict(model2, 
                                 newdata = data.frame(smoking_years = smoking_grid[i], age = age_grid[j]), 
                                 type = "response")
  }
}

# ヒートマップを描画
image(smoking_grid, age_grid, prob_matrix, col = rev(heat.colors(100)), 
      xlab = "喫煙期間（年）", ylab = "年齢（歳）", main = "喫煙期間・年齢と肺がん発生確率")

# 肺がんが発生した人（1）と発生しなかった人（0）をプロット
points(smoking_years[lung_cancer == 1], age[lung_cancer == 1], pch = 17, col = 1)  # 肺がん発生者
points(smoking_years[lung_cancer == 0], age[lung_cancer == 0], pch = 16, col = 4)  # 肺がん非発生者

# 凡例
legend("bottomright", legend = c("肺がん発生", "肺がんなし"), pch = c(17, 16), col = c(1,4))

Rスクリプトの説明

　上のRスクリプトでは、喫煙期間（年数）と年齢が肺がん発生確率にどのように影響するかをロジスティック回帰で分析し、ヒートマップで可視化しています。以下で、各コマンドと変数について解説します。

1. データの準備

# データ数
n <- 100

100人分の仮想データを作成するため、サンプル数を 100 に設定。

# 年齢をランダムに生成（30〜80歳）
age <- runif(n, 30, 80)

runif(n, 30, 80) は、30〜80歳の範囲で100個のランダムな値（年齢）を一様分布に従い生成しています。

# 喫煙期間（年数）をランダムに生成（0〜40年）
smoking_years <- runif(n, 0, 40)
smoking_years[smoking_years > age - 20] <- smoking_years[smoking_years > age - 20] - 20

runif(n, 0, 40) は、0〜40年の範囲で100個のランダムな値（喫煙年数）を一様分布に従うように生成します。ただし、喫煙年数が、実際の年齢よりも長いのは不自然ですので、喫煙年数が年齢よりも長いものを修正しています。

2. 肺がん発生確率の計算

# 肺がんになる確率をロジスティック関数で計算（喫煙年数と年齢を考慮）
prob <- 1 / (1 + exp(-(-6 + 0.1 * smoking_years + 0.08 * age)))

ロジスティック関数を使い、喫煙年数と年齢が肺がんの発生確率に与える影響を計算しています。

$\displaystyle{ P(y=1 | x_1, x_2) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1 + \beta_2 x_2)}} }$

3. 確率に基づく肺がん発生データの作成

# 確率に基づいて肺がんになったかどうかを決定（0 or 1）
lung_cancer <- rbinom(n, 1, prob)

rbinom(n, 1, prob) で、確率 prob に基づいて肺がんの有無（0: なし, 1: あり）を２項分布に従って決定しています。

4. ロジスティック回帰モデルの構築

# ロジスティック回帰モデルの構築（共変量2つ）
model2 <- glm(lung_cancer ~ smoking_years + age, family = binomial)

glm()（Generalized Linear Model：一般化線形モデル）を使用してロジスティック回帰を実行しています。

lung_cancer ~ smoking_years + age:
- 目的変数（従属変数）: lung_cancer（肺がんの有無）
- 説明変数（独立変数）: smoking_years（喫煙期間）と age（年齢）
family = binomial:
- ロジスティック回帰（0-1の分類問題）を指定

5. ヒートマップの作成

# グリッド状のデータを作成（喫煙期間と年齢の範囲）
smoking_grid <- seq(0, 40, length.out = 50)
age_grid <- seq(30, 80, length.out = 50)

グリッドデータの作成しています。

seq(0, 40, length.out = 50): 0〜40年の喫煙期間を50ステップで分割
seq(30, 80, length.out = 50): 30〜80歳の年齢を50ステップで分割

# ヒートマップ用の確率行列を作成
prob_matrix <- matrix(0, nrow = length(smoking_grid), ncol = length(age_grid))

確率を格納するために、50×50 の空の行列を作成しています。

# 各グリッド点での肺がん発生確率を計算
for (i in 1:length(smoking_grid)) {
  for (j in 1:length(age_grid)) {
    prob_matrix[i, j] <- predict(model2, 
                                 newdata = data.frame(smoking_years = smoking_grid[i], age = age_grid[j]), 
                                 type = "response")
  }
}

prob_matrix の各要素に予測確率を代入しています。predict() を使い、喫煙期間と年齢の各グリッド点で肺がん発生確率を計算しています。

6. ヒートマップの描画

# ヒートマップを描画（パッケージ不要のbase R）
image(smoking_grid, age_grid, prob_matrix, col = rev(heat.colors(100)), 
      xlab = "喫煙期間（年）", ylab = "年齢（歳）", main = "喫煙期間・年齢と肺がん発生確率")

rev(heat.colors(100)) で色をつけ、発生確率をヒートマップとして描いています。

まとめ

ロジスティック回帰は確率を予測するための回帰モデル
シグモイド関数を使って、確率を0～1の範囲に変換する
Rの glm() 関数を使えば簡単に実装できる！

ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【Rで医療統計】ロジスティック回帰とは：共変量が1つか２つの例

ロジスティック回帰とは？

ロジスティック回帰の実際 (架空の喫煙者の分析)

共変量が1つの例（喫煙期間が長いと、肺がんになりやすい？）

Rでロジスティック回帰

Rスクリプトの説明

1. データの準備

2. 肺がん発生確率の計算

3. 確率に基づく肺がん発生データの作成

4. ロジスティック回帰モデルの構築

5. 結果の表示

6. 予測データの作成

7. グラフの描画

共変量が2つの例（喫煙期間だけでなく、年齢も肺がんになりやすさと関係する？）

Rでロジスティック回帰

Rスクリプトの説明

1. データの準備

2. 肺がん発生確率の計算

3. 確率に基づく肺がん発生データの作成

4. ロジスティック回帰モデルの構築

5. ヒートマップの作成

6. ヒートマップの描画

まとめ