[艦これで学ぶ統計学] 第4回 回帰分析

今回から何回かに分けて回帰分析について解説します。以下のブロマガを題材にして解説します。

• ［艦これ］夜偵検証結果と発動率の推定式:うみゃのブロマガ - ブロマガ

このブロマガにあるデータを用いて、夜偵を1隊だけ積んだ場合の夜偵発動率を夜偵を積んだ艦娘のレベルごとに集計すると以下のようになっていました。

このデータを用いて、各レベルごとに夜偵発動率がどの程度になるか推定してみたいと思います。回帰分析という処理を行うことで、レベルごとの夜偵発動率の推定式を立てることができます。

Rでの回帰分析は、データをまとめたcsvファイルを用いて行うことができます。csvファイルは表計算ソフトを用いて作ることができます。本記事での解説に使用するデータを以下から開いてください。

• [艦これで学ぶ統計学] 第4回 単回帰分析 - Google スプレッドシート

「レベル」列が夜偵を積んだ艦娘のレベルを表し、「発動可否」列が1なら夜偵発動を、0なら夜偵不発を表しています。

ページ内から「ファイル」→「形式を指定してダウンロード」→「カンマ区切りの値(.csv、現在のシート)」を選び、csvファイルをダウンロードしてください。

ダウンロードしたcsvファイルをRで読み込みます。まずRで使用するフォントを日本語に適したものに変更します。Rから「編集」→「GUI プリファレンス...」を選択して「Rgui 設定エディター」を開き、「Font」を「MS Gothic」に変更し、「Save...」したのち「OK」を選択してください。

Rでcsvファイルを開きます。csvファイルはフルパスで指定してください。また、「\」は二つ重ねて入力する必要があるので注意してください。

data <- read.csv(fileEncoding = "UTF-8", "C:\\Users\\Master\\Downloads\\艦これで学ぶ統計学 - 第4回 回帰分析.csv")

head(data)

1行目でcsvファイルを読み込んでいます。fileEncodingではcsvファイルの文字コードを指定しています。Google スプレッドシートからダウンロードしたcsvファイルは文字コードがUTF-8なので、fileEncoding="UTF-8"と指定する必要があります。Excelから作成したcsvファイルの場合は、単にdata <- read.csv("<csvファイルのフルパス>")とすれば読み込むことができます。

2行目でcsvファイルの先頭部分を表示しています。csvファイルのデータすべてを表示する場合は単にdataと入力します。

このデータを用いて回帰分析を行います。推定式として以下の式を立てます。

夜偵発動率 = a + b × レベル

このようにして、レベルによって夜偵発動率が変わる式を立てることができます。aが定数項、bがレベルの係数になります。回帰分析をすることで、元データに最も適したa, bを計算することができます。

Rで回帰分析を行います。Rに以下のように入力します。

result <- lm(発動可否 ~ レベル, data)

summary(result)

1行目で回帰分析を行い、2行目でその結果を表示しています。

この結果から、レベルごとの夜偵発動率が以下のように推定されます。

夜偵発動率 = 0.1766427 + 0.0046722 × レベル

画像の赤枠にある数値を用いてこのように推定式を立てることができます。(Intercept)にある数値が推定式のaにあたる定数項を表し、レベルにある数値が推定式のbにあたるレベルの係数を表しています。

緑枠にある数値は標準誤差とよばれるもので、各係数がどの程度ぶれているかを表しています。この標準誤差を推定式に組み込むと以下のようになります。

夜偵発動率 = (0.1766427±0.0120462) + (0.0046722±0.0001158) × レベル

標準誤差は推定式がどの程度ぶれているかの目安になるため、なるべく書き添えるようにしましょう。

この標準誤差を約1.959964倍すると95%信頼区間を計算できます。レベルの係数の場合を例に説明すると、95%信頼区間の下限は0.0046722 - 0.0001158 × 1.959964 = 0.004445236、上限は0.0046722 + 0.0001158 × 1.959964 = 0.004899164となります。

得られた推定式のグラフを書いてみます。Rに以下のように入力します。

result.x <- 1:150

result.newdata <- data.frame(レベル = result.x)

result.y <- predict(result, result.newdata)

plot(result.x, result.y, xlim = c(1, 150), ylim = c(0, 1), xlab = "レベル", ylab = "夜偵発動率", type = "l")

これでレベル1以上150以下の場合の推定式のグラフを書くことができます。

この推定式のグラフがどの程度当てはまっているか確認するために、上の表にある各レベルごとのデータを重ねて表示してみます。グラフを表示したのちに、以下のように入力します。

error.bar <- function (x, y, ...) {

    x.mean <- mean(x)

    y.mean <- mean(y)

    y.se <- sqrt(var(y) / length(y))

    points(x.mean, y.mean, ...);

    arrows(x.mean, y.mean - y.se, x.mean, y.mean + y.se, angle = 90, code = 3, length = 0.05, ...)

}

data.1 <- subset(data, レベル == 1)

error.bar(data.1$レベル, data.1$発動可否)

data.50.55 <- subset(data, 50 <= レベル & レベル <= 55)

error.bar(data.50.55$レベル, data.50.55$発動可否)

data.99 <- subset(data, レベル == 99)

error.bar(data.99$レベル, data.99$発動可否)

data.116.125 <- subset(data, 116 <= レベル & レベル <= 125)

error.bar(data.116.125$レベル, data.116.125$発動可否)

data.150 <- subset(data, レベル == 150)

error.bar(data.150$レベル, data.150$発動可否)

data.1 <- subset(data, レベル == 1)でcsvデータからレベル1のデータのみを抜き出しています。途中、レベル == 1の部分は「=」を重ねていることに注意してください。同様に、data.50.55 <- subset(data, 50 <= レベル & レベル <= 55)でレベル50以上55以下のデータを抜き出しています。

このようにして抜き出したのち、error.bar(data.1$レベル, data.1$発動可否)のようにすることで抜き出したデータをグラフに書き込んでいます。

ここで表示したものをエラーバーといいます。この縦の幅が各点の標準誤差を表しています。

今回は以上になります。次回以降で、よりよい推定式が他にないか模索していきます。最後に、本記事で使用したコードの全文を以下に書いておきます。

data <- read.csv(fileEncoding = "UTF-8", "C:\\Users\\Master\\Downloads\\艦これで学ぶ統計学 - 第4回 回帰分析.csv")

head(data)

result <- lm(発動可否 ~ レベル, data)

summary(result)

result.x <- 1:150

result.newdata <- data.frame(レベル = result.x)

result.y <- predict(result, result.newdata)

plot(result.x, result.y, xlim = c(1, 150), ylim = c(0, 1), xlab = "レベル", ylab="夜偵発動率", type = "l")

error.bar <- function (x, y, ...) {

    x.mean <- mean(x)

    y.mean <- mean(y)

    y.se <- sqrt(var(y) / length(y))

    points(x.mean, y.mean, ...);

    arrows(x.mean, y.mean - y.se, x.mean, y.mean + y.se, angle = 90, code = 3, length = 0.05, ...)

}

data.1 <- subset(data, レベル == 1)

error.bar(data.1$レベル, data.1$発動可否)

data.50.55 <- subset(data, 50 <= レベル & レベル <= 55)

error.bar(data.50.55$レベル, data.50.55$発動可否)

data.99 <- subset(data, レベル == 99)

error.bar(data.99$レベル, data.99$発動可否)

data.116.125 <- subset(data, 116 <= レベル & レベル <= 125)

error.bar(data.116.125$レベル, data.116.125$発動可否)

data.150 <- subset(data, レベル == 150)

error.bar(data.150$レベル, data.150$発動可否)