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

第4回から続く回帰分析の解説の続きです。今回は夜偵を複数詰んだ場合について扱います。夜偵を複数詰んだ場合の集計結果は以下のようになっています。

夜偵を複数詰んだ場合の元データを以下におきました。csvファイルとしてダウンロードしてください。

• 艦これで学ぶ統計学 - Google スプレッドシート

夜偵の搭載スロット数も考慮して回帰分析を行います。以下のような推定式を立てます。

夜偵発動率 = a + b × √(レベル) + c × 搭載スロット数

レベルに加えて搭載スロット数も推定式に加えました。レベルと搭載スロット数といったように複数の要素から行う回帰分析を重回帰分析といいます。

Rでは重回帰分析も通常の回帰分析と同様に行うことができます。Rに以下のように入力します。

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

result <- lm(発動可否 ~ sqrt(レベル) + 搭載スロット数, data)

summary(result)

推定式は以下のようになります。

夜偵発動率 = (-0.190640±0.019333) + (0.078152±0.001287) × √(レベル) + (0.097281 ± 0.009084) × 搭載スロット数

搭載スロット数が1の場合の推定式のグラフを書いてみます。Rに以下のように入力します。

result.x <- 1:150

result.newdata.1 <- data.frame(レベル = result.x, 搭載スロット数 = 1)

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

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

result.newdata.1 <- data.frame(レベル = result.x, 搭載スロット数 = 1)とすることで、搭載スロット数が1の場合のグラフを書くことができます。

同様に、搭載スロット数が2, 3の場合の推定式のグラフを重ねて書いてみます。Rに以下のように入力します。

par(new = TRUE)

result.newdata.2 <- data.frame(レベル = result.x, 搭載スロット数 = 2)

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

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

par(new = TRUE)

result.newdata.3 <- data.frame(レベル = result.x, 搭載スロット数 = 3)

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

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

par(new = TRUE)

legend("bottomright", legend = c("搭載スロット数1", "搭載スロット数2", "搭載スロット数3"), lty = c(1, 2, 3), col = c(1, 2, 3))

plotする際にlty = 2, lty = 3とすることで線の種類を変え、col = 2, col = 3とすることで線の色を変えることができます。

legend("bottomright", legend = c("搭載スロット数1", "搭載スロット数2", "搭載スロット数3"), lty = c(1, 2, 3), col = c(1, 2, 3))とすることで判例を書くことができます。"bottomright"は判例を書く位置を表します。"lefttop"で左上に、"righttop"で右上に、"leftbottom"で左下に、"rightbottom"で右下に書くことができます。lty = c(1, 2, 3)とcol = c(1, 2, 3)にはplotした際に指定したltyとcolを順番に指定します。

元データの値も図に書き込みます。Rに以下のように入力します。

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.1 <- subset(data, 搭載スロット数 == 1 & レベル == 1)

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

data.1.50.55 <- subset(data, 搭載スロット数 == 1 & 50 <= レベル & レベル <= 55)

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

data.1.99 <- subset(data, 搭載スロット数 == 1 & レベル == 99)

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

data.1.116.125 <- subset(data, 搭載スロット数 == 1 & 116 <= レベル & レベル <= 125)

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

data.1.150 <- subset(data, 搭載スロット数 == 1 & レベル == 150)

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

data.2.1 <- subset(data, 搭載スロット数 == 2 & レベル == 1)

error.bar(data.2.1$レベル, data.2.1$発動可否, col = 2)

data.2.99 <- subset(data, 搭載スロット数 == 2 & レベル == 99)

error.bar(data.2.99$レベル, data.2.99$発動可否, col = 2)

data.2.150 <- subset(data, 搭載スロット数 == 2 & レベル == 150)

error.bar(data.2.150$レベル, data.2.150$発動可否, col = 2)

data.3.1 <- subset(data, 搭載スロット数 == 3 & レベル == 1)

error.bar(data.3.1$レベル, data.3.1$発動可否, col = 3)

data.1.1 <- subset(data, 搭載スロット数 == 1 & レベル == 1)において& 搭載スロット数 == 1とすることで、搭載スロット数が1のデータを抽出できます。error.bar(data.1.1$レベル, data.1.1$発動可否, col = 1)においてcol = 1と指定することで、エラーバーの色を指定できます。

AICの計算もしてみます。Rに以下のように入力します。

AIC(result)

ここで計算したAICは、前回計算したAICと比較することはできません。AICは元データが同じ場合にのみ比較できます。前回は搭載スロット数が1のみのデータで、今回は搭載スロット数が2, 3のデータも含まれているため、比較することができません。

最後に今回使用したRのコードを載せておきます。

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

result <- lm(発動可否 ~ sqrt(レベル) + 搭載スロット数, data)

summary(result)

result.x <- 1:150

result.newdata.1 <- data.frame(レベル = result.x, 搭載スロット数 = 1)

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

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

par(new = TRUE)

result.newdata.2 <- data.frame(レベル = result.x, 搭載スロット数 = 2)

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

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

par(new = TRUE)

result.newdata.3 <- data.frame(レベル = result.x, 搭載スロット数 = 3)

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

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

par(new = TRUE)

legend("bottomright", legend = c("搭載スロット数1", "搭載スロット数2", "搭載スロット数3"), lty = c(1, 2, 3), col = c(1, 2, 3))

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.1 <- subset(data, 搭載スロット数 == 1 & レベル == 1)

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

data.1.50.55 <- subset(data, 搭載スロット数 == 1 & 50 <= レベル & レベル <= 55)

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

data.1.99 <- subset(data, 搭載スロット数 == 1 & レベル == 99)

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

data.1.116.125 <- subset(data, 搭載スロット数 == 1 & 116 <= レベル & レベル <= 125)

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

data.1.150 <- subset(data, 搭載スロット数 == 1 & レベル == 150)

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

data.2.1 <- subset(data, 搭載スロット数 == 2 & レベル == 1)

error.bar(data.2.1$レベル, data.2.1$発動可否, col = 2)

data.2.99 <- subset(data, 搭載スロット数 == 2 & レベル == 99)

error.bar(data.2.99$レベル, data.2.99$発動可否, col = 2)

data.2.150 <- subset(data, 搭載スロット数 == 2 & レベル == 150)

error.bar(data.2.150$レベル, data.2.150$発動可否, col = 2)

data.3.1 <- subset(data, 搭載スロット数 == 3 & レベル == 1)

error.bar(data.3.1$レベル, data.3.1$発動可否, col = 3)

AIC(result)