R Markdown

データを準備。

pheno <- read.csv("RiceDiversityPheno.csv")
line <- read.csv("RiceDiversityLine.csv")
line.pheno <- merge(line, pheno, by.x = "NSFTV.ID", by.y = "NSFTVID")
data <- data.frame(
    height = line.pheno$Plant.height,
    flower = line.pheno$Flowering.time.at.Arkansas)
data <- na.omit(data)
x <- data$flower
y <- data$height

回帰係数を計算。

n <- length(x)
ssx <- sum(x^2) - n * mean(x)^2
ssxy <- sum(x * y) - n * mean(x) * mean(y)
b <- ssxy / ssx
b
## [1] 0.6728746
m <- mean(y) - b * mean(x)
m
## [1] 58.05464

MSEを計算しておく。

ssr <- b * ssxy
ssr
## [1] 26881.49
ssy <- sum(y^2) - n * mean(y)^2
sse <- ssy - ssr
sse
## [1] 133903.2
mse <- sse / (n - 2)
mse
## [1] 360.9251

\(H_0: \beta_0 = 0.5\)を検定する。

t.value <- (b - 0.5) / sqrt(mse/ssx)
t.value
## [1] 2.217253
2 * (1 - pt(t.value, n - 2))
## [1] 0.02721132

\(\frac {b - 0.5} {\sqrt{MSE/SSX}}\)が従う自由度\(n-2\)\(t\)分布を描く。

xx <- seq(-5, 5, 0.01)
tt <- dt(xx, n - 2)
plot(xx, tt, type = "l")
# calculate the 100 * (1 - alpha/2) percentile 
t.975 <- qt(1 - 0.025, n - 2)
t.975
## [1] 1.966379
# calculate the 100 * alpha / 2 percentile
t.025 <- qt(0.025, n - 2)
t.025
## [1] -1.966379
# the above value can be calculated also as
- qt(1 - 0.025, n - 2)
## [1] -1.966379
# this is because the shape t-distribution is symmetric
# draw lines of the bounds
abline(v = t.025, col = "green", lty = "dotted")
abline(v = t.975, col = "blue", lty = "dotted")
# draw the line of t.value
abline(v = t.value, col = "red", lty = "dotted")

緑色の線から青色の線までが、この分布が95%の確率で含む値の範囲。赤色の点線は、その外側にある。

つぎに、\(t\)分布の2.5%点(\(t_{n-2, 0.025}\))と97.5%点(\(t_{n-2, 0.975}\))が与えられたときの、\(beta_0\)が95%の確率で収まる範囲を求める。\(beta_0\)がこの範囲に収まるとき、以下の式が成り立つ。 \[ t_{n-2, 0.025} \le \frac {b - \beta_0} {\sqrt{MSE/SSX}} \le t_{n-2, 0.975} \] ここで、2.5%点は、97.5%点に\(-1\)を乗じたものになっているので、 \[ - t_{n-2, 0.975} \le \frac {b - \beta_0} {\sqrt{MSE/SSX}} \le t_{n-2, 0.975} \] \[ - t_{n-2, 0.975} {\sqrt{MSE/SSX}} \le b - \beta_0 \le t_{n-2, 0.975} {\sqrt{MSE/SSX}} \] \[ -b - t_{n-2, 0.975} {\sqrt{MSE/SSX}} \le - \beta_0 \le -b + t_{n-2, 0.975} {\sqrt{MSE/SSX}} \] \[ b - t_{n-2, 0.975} {\sqrt{MSE/SSX}} \le \beta_0 \le b + t_{n-2, 0.975} {\sqrt{MSE/SSX}} \] 上式より、\(\beta_0\)の信頼区間が得られる。