［R］回帰モデルの信頼区間と予測区間

lm関数により求めた回帰モデルの信頼区間と予測区間はpredict関数で求めることができる。以下の計算により、赤実線が求めた回帰直線。桃破線がその95%信頼区間、橙点線がその95%予測区間。

> # 説明変数xと目的変数y
> x <- c(1, 2, 3, 5, 6, 7, 8, 10, 11, 12)
> y <- c(0, 1, 1, 2, 4, 4, 6,  5,  8,  8)
> # 回帰直線を求める
> r <- lm(y ~ x)
> print(r)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept)            x  
    -0.8567       0.7318  
> # yの推定値
> yest <- fitted(r)
> # 信頼区間と予測区間の計算(後続の手順のためdata.frame化)
> dtf <- data.frame(x)
> rcon <- predict(r, dtf, interval = "confidence", level = 0.95)
> rpre <- predict(r, dtf, interval = "prediction", level = 0.95)
> rcon <- data.frame(rcon)
> rpre <- data.frame(rpre)
> # 図の作成(黒丸:観測値,赤実線:回帰直線,桃破線:信頼区間,橙点線:予測区間)
> plot(x, y, type = "n", asp = 1.0)
> lines(x, rpre$lwr, lty = "dotted", col = "orange", lwd = 2.0)
> lines(x, rpre$upr, lty = "dotted", col = "orange", lwd = 2.0)
> lines(x, rcon$lwr, lty = "dashed", col = "pink", lwd = 2.0)
> lines(x, rcon$upr, lty = "dashed", col = "pink", lwd = 2.0)
> lines(x, yest, col = "red", lwd = 2.0)
> points(x, y, pch = 20, col = "black")

> # 以下はpredict関数と手計算の計算結果の比較
> ah <- 0.95
> n <- length(x)
> k <- 2
> mxxx <- matrix(c(rep(1, n), x), ncol = 2)
> mxy <- matrix(y, ncol = 1)
> # 最小二乗推定量
> mxb <- solve(t(mxxx) %*% mxxx) %*% (t(mxxx) %*% mxy)
> print(mxb)
           [,1]
[1,] -0.8567050
[2,]  0.7318008
> # yの推定値
> yest <- mxxx %*% mxb
> # 残差分散
> s <- sqrt(sum((y - yest) ^ 2) / (n - k))
> #
> s0 <- sh0 <- double(n)
> for (i in 1:n) {
+ mxx0 <- matrix(c(1, x[i]), ncol = 1)
+ s0[i] <- s * sqrt(t(mxx0) %*% solve(t(mxxx) %*% mxxx) %*% mxx0)
+ sh0[i] <- s * sqrt(1 + t(mxx0) %*% solve(t(mxxx) %*% mxxx) %*% mxx0)
+ }
> # 95%信頼区間
> ycon1 <- yest - s0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> ycon2 <- yest + s0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> # 95%予測区間
> ypre1 <- yest - sh0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> ypre2 <- yest + sh0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> # predict関数による信頼区間と予測区間
> print(data.frame(rcon$lwr, rcon$upr, rpre$lwr, rpre$upr))
     rcon.lwr  rcon.upr   rpre.lwr  rpre.upr
1  -1.1763909 0.9265825 -2.2315171  1.981709
2  -0.3152118 1.5290049 -1.4382140  2.652007
3   0.5349489 2.1424457 -0.6558465  3.333241
4   2.1772621 3.4273356  0.8728263  4.731771
5   2.9513451 4.1168542  1.6179064  5.450293
6   3.6831458 4.8486549  2.3497072  6.182094
7   4.3726644 5.6227379  3.0682286  6.927174
8   5.6575543 7.2650511  4.4667589  8.455846
9   6.2709951 8.1152118  5.1479929  9.238214
10  6.8734175 8.9763909  5.8182913 10.031517
> # 手計算による信頼区間と予測区間
> print(data.frame(ycon1, ycon2, ypre1, ypre2))
        ycon1     ycon2      ypre1     ypre2
1  -1.1763909 0.9265825 -2.2315171  1.981709
2  -0.3152118 1.5290049 -1.4382140  2.652007
3   0.5349489 2.1424457 -0.6558465  3.333241
4   2.1772621 3.4273356  0.8728263  4.731771
5   2.9513451 4.1168542  1.6179064  5.450293
6   3.6831458 4.8486549  2.3497072  6.182094
7   4.3726644 5.6227379  3.0682286  6.927174
8   5.6575543 7.2650511  4.4667589  8.455846
9   6.2709951 8.1152118  5.1479929  9.238214
10  6.8734175 8.9763909  5.8182913 10.031517

［R］母比率p（バストが85cm以上の女子大生の割合）を95%信頼区間で区間推定（「統計解析のはなし」（東京図書）pp.xxx-xxx）

> x <- c(80, 83, 83, 83, 81, 78, 83, 83, 80, 85)
> x <- c(x, 82, 80, 90, 85, 80, 82, 87, 81, 83, 85)
> nn <- length(x)
> ah <- 0.05
> m <- 5
> z <- qnorm(ah / 2, lower.tail = FALSE)
> a <- m / nn - z * sqrt(m / nn * (1 - m / nn) / nn)
> b <- m / nn + z * sqrt(m / nn * (1 - m / nn) / nn)
> # N=20におけるバスト85cm以上の女子大生の母比率pの95%信頼区間
> # (標本数Nが大きい場合の公式による)
> cat(sprintf("%.2f≦p≦%.2f\n", a, b))
0.06≦p≦0.44
> d1 <- 2 * (nn - m + 1)
> d2 <- 2 * m
> e1 <- 2 * (m + 1)
> e2 <- 2 * (nn - m)
> fd <- qf(ah / 2, d1, d2, lower.tail = FALSE)
> fe <- qf(ah / 2, e1, e2, lower.tail = FALSE)
> a <- d2 / (d1 * fd + d2)
> b <- e1 * fe / (e1 * fe + e2)
> # N=20におけるバスト85cm以上の女子大生の母比率pの95%信頼区間
> # (標本数Nが小さい場合の公式による)
> cat(sprintf("%.2f≦p≦%.2f\n", a, b))
0.09≦p≦0.49
> nn <- 200
> m <- 45
> a <- m / nn - z * sqrt(m / nn * (1 - m / nn) / nn)
> b <- m / nn + z * sqrt(m / nn * (1 - m / nn) / nn)
> # N=200におけるバスト85cm以上の女子大生の母比率pの95%信頼区間
> cat(sprintf("%.3f≦p≦%.3f\n", a, b))
0.167≦p≦0.283

［R］スコップの柄の直径の2つの母分散の比の90%信頼区間（「統計解析のはなし」（東京図書）pp.167-168）

> x1 <- c(5.47, 5.62, 5.52, 5.77, 5.53, 5.39)
> x2 <- c(5.63, 5.46, 5.53, 5.54, 5.59)
> x1m <- mean(x1)
> x2m <- mean(x2)
> s12 <- var(x1)
> s22 <- var(x2)
> nn1 <- length(x1)
> nn2 <- length(x2)
> degf1 <- nn1 - 1
> degf2 <- nn2 - 1
> ah <- 0.1
> ff1 <- qf(1 - ah / 2, degf1, degf2, lower.tail = FALSE)
> ff2 <- qf(ah / 2, degf1, degf2, lower.tail = FALSE)
> a <- s12 / s22 / ff2
> b <- s12 / s22 / ff1
> # スコップの柄の直径の2つの母分散の比の90%信頼区間
> cat(sprintf("%.2f≦σ1^2/σ2^2≦%.2f\n", a, b))
0.67≦σ1^2/σ2^2≦21.67

［R］自動車の燃費の母分散の95%信頼区間（「統計解析のはなし」（東京図書）pp.165-166）

> x <- c(15.4, 16.1, 15.7, 16.6, 14.9, 15.5, 16.2)
> nn <- length(x)
> xm <- mean(x)
> s2 <- var(x)
> degf <- nn - 1
> ah <- 0.05
> ca <- qchisq(1 - ah / 2, degf, lower.tail = FALSE)
> cb <- qchisq(ah / 2, degf, lower.tail = FALSE)
> a <- (nn - 1) * s2 / cb
> b <- (nn - 1) * s2 / ca
> # 自動車の燃費の母分散の95%信頼区間
> cat(sprintf("%.2f≦σ^2≦%.2f\n", a, b))
0.14≦σ^2≦1.58

［R］女子大生200人の身長の母分散の90%信頼区間（「統計解析のはなし」（東京図書）pp.163-164）

> x <- c(159, 158, 151, 167, 151)
> x <- c(x, 160, 160, 158, 160, 158)
> xm <- mean(x)
> nn <- length(x)
> s2 <- var(x)
> degf <- nn - 1
> ah <- 0.1
> ca <- qchisq(1 - ah / 2, degf, lower.tail = FALSE)
> cb <- qchisq(ah / 2, degf, lower.tail = FALSE)
> a <- (nn - 1) * s2 / cb
> b <- (nn - 1) * s2 / ca
> # 女子大生200人の身長の母分散の90%信頼区間
> cat(sprintf("%.2f≦σ^2≦%.2f\n", a, b))
11.32≦σ^2≦57.62

［R］新型ナットの平均内径の99%信頼区間（「統計解析のはなし」（東京図書）pp.154-155）

> x <- c(1.29, 1.34, 1.31, 1.30, 1.30)
> nn <- length(x)
> degf <- nn - 1
> xm <- mean(x)
> s <- 0.015
> ah <- 0.01
> z <- qnorm(ah / 2, mean = 0, sd = 1, lower.tail = FALSE)
> a <- xm - z * s / sqrt(nn)
> b <- xm + z * s / sqrt(nn)
> # 新型ナットの平均内径の99%信頼区間
> cat(sprintf("%.3f≦μ≦%.3f\n", a, b))
1.291≦μ≦1.325

［R］T大学の女子大生の平均I.Q.の95%信頼区間（「統計解析のはなし」（東京図書）p.154）

> nn <- 500
> ah <- 0.05
> xm <- 112
> s <- 9.5
> z <- qnorm(ah / 2, mean = 0, sd = 1, lower.tail = FALSE)
> a <- xm - z * s / sqrt(nn)
> b <- xm + z * s / sqrt(nn)
> # T大学の女子大生の平均I.Q.の95%信頼区間
> cat(sprintf("%.1f≦μ≦%.1f\n", a, b))
111.2≦μ≦112.8

［R］サラダ油800g缶の平均内容量の99%信頼区間（「統計解析のはなし」（東京図書）p.153）

> x <- c(807, 811, 801, 798, 798, 795, 803, 805, 804)
> nn <- length(x)
> degf <- nn - 1
> ah <- 0.01
> xm <- mean(x)
> s <- sd(x)
> tval <- qt(ah / 2, degf, lower.tail = FALSE)
> a <- xm - tval * s / sqrt(nn)
> b <- xm + tval * s / sqrt(nn)
> # サラダ油800g缶の平均内容量の99%信頼区間
> cat(sprintf("%.1f≦μ≦%.1f\n", a, b))
796.8≦μ≦808.0

［Anaconda］仮想環境の一覧を得る

condaコマンドを使う。

>conda env list
# conda environments:
#
base                  *  C:\Users\○○○\anaconda3

［R］乗用車のレギュラーガソリンの平均売り上げ量の95%信頼区間（「統計解析のはなし」（東京図書）pp.152-153）

> x <- c(45, 39, 42, 57, 28, 33, 40, 51)
> nn <- length(x)
> degf <- nn - 1
> ah <- 0.05
> xm <- mean(x)
> s <- sd(x)
> tval <- qt(ah / 2, degf, lower.tail = FALSE)
> a <- xm - tval * s / sqrt(nn)
> b <- xm + tval * s / sqrt(nn)
> # 乗用車のレギュラーガソリンの平均売り上げ量の95%信頼区間
> cat(sprintf("%.1f≦μ≦%.1f\n", a, b))
34.1≦μ≦49.6

［R］アオスジアゲハ夏型の体長の90%信頼区間（「統計解析のはなし」（東京図書）pp.151-152）

> x <- c(76, 85, 82, 83, 76, 78)
> nn <- length(x)
> xm <- mean(x)
> s <- sd(x)
> degf <- nn - 1
> ah <- 0.1
> tval <- qt(ah / 2, degf, lower.tail = FALSE)
> a <- xm - tval * s / sqrt(nn)
> b <- xm + tval * s / sqrt(nn)
> # アオスジアゲハ夏型の体長の90%信頼区間
> cat(sprintf("%.3f≦μ≦%.3f\n", a, b))
76.835≦μ≦83.165

［R］溶液のpH値の99%信頼区間（「統計解析のはなし」（東京図書）pp.150-151）

> x <- c(7.86, 7.89, 7.84, 7.90, 7.82)
> nn <- length(x)
> xm <- mean(x)
> s <- sd(x)
> degf <- nn - 1
> ah <- 0.01
> tval <- qt(ah / 2, degf, lower.tail = FALSE)
> a <- xm - tval * s / sqrt(nn)
> b <- xm + tval * s / sqrt(nn)
> # 溶液のpHの99%信頼区間
> cat(sprintf("%.3f≦μ≦%.3f\n", a, b))
7.793≦μ≦7.931

［R］文字列の文字コードを変換する

iconv関数を使う。以下の例では、シフトSJISの環境下で、以下の3行からなるUTF-8（改行コードはCR+LF）で保存されたテキストファイル（ファイル名はutf8.txt）を直接読み込み、それをUTF-8からシフトJISに保存をしている。

ABC
あいう
阿位宇

ファイルを読み込み、文字列を変換する。

> scan("utf8.txt", what = character(), fileEncoding = "UTF-8")
Read 3 items
[1] "ABC"    "あいう" "阿位宇"
> ch <- readChar("utf8.txt", 27, useBytes = TRUE)
> print(ch)
[1] "ABC\r\n縺ゅ＞縺<86>\r\n髦ｿ菴榊ｮ<87>\r\n"
> iconv(ch, from = "UTF-8", to = "SJIS")
[1] "ABC\r\nあいう\r\n阿位宇\r\n"

［R］環境変数を取得する

Sys.getenv関数を使う。

> Sys.getenv("HOMEDRIVE")
[1] "C:"
> Sys.getenv("COMSPEC")
[1] "C:\\Windows\\system32\\cmd.exe"