[艦これで学ぶ統計学] 第5回 AIC

前回に引き続き回帰分析の解説を行います。あらかじめ前回のRのコードをRに入力したうえで以下の記事をお読みください。

前回立てた推定式とそのグラフを再掲します。

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

これを見るとレベルが低いときに推定式のグラフが元データからずれているようです。そこでより適した推定式がほかにないか考えてみます。ここで夜偵発動率がレベルの平方根に応じて変化するという仮説を立て、以下の推定式をもとに考えてみます。

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

この推定式をもとに回帰分析を行います。Rに以下のように入力します。

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

summary(sqrt.result)

sqrt(レベル)の部分が平方根を表しています。この結果から推定式は以下のようになります。

夜偵発動率 = (-0.041632±0.016555) + (0.071799±0.001727) × √(レベル)

この推定式のグラフを前回書いた図に重ねてみます。Rに以下のように入力します。

par(new = TRUE)

sqrt.result.x <- 1:150

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

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

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

1行目のようにpar(new = TRUE)と入力することで、グラフを重ねることができます。これを入力しないと今まで書いたグラフが消えてしまうため注意してください。

平方根による推定式は低レベルでもよくあてはまっているように見えます。これを目測ではなく、統計を用いてどちらの推定式があてはまっているかを調べてみます。Rに以下のように入力します。

AIC(result)

AIC(sqrt.result)

ここで表示した数値をAICといい、このAICの数値が小さいほどあてはまりがよいことを表しています。平方根の推定式のほうがAICが小さいため、よりあてはまりがよいことがわかります。このように、推定式を比較する場合はそのAICを計算し、AICがより小さい推定式がないか探します。

最後に、今回使用したRのコードを以下に再掲します。

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$発動可否)

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

summary(sqrt.result)

par(new = TRUE)

sqrt.result.x <- 1:150

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

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

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

AIC(result)

AIC(sqrt.result)