「データ解析」（下平英寿） 講義資料7 検定と信頼区間 1 回帰係数のt-検定

Copyright (c) 2004,2005 Hidetoshi Shimodaira
2005-04-13 17:37:08 shimo
「データ解析」（下平英寿）
• 不偏推定量 S²2 を用いると，
V̂ (β̂i ) = S²2 aii
講義資料７
である．したがって標準誤差の推定は
検定と信頼区間
se(
ˆ β̂i ) =
• 目標： 検定と信頼区間について理解する．
q
√
V̂ (β̂i ) = S² aii
Ｒの組み込み関数では，この式を用いている．
1. 回帰係数の t-検定
• t-統計量，p-値は
2. 回帰係数，予測値の信頼区間
ti =
3. 回帰モデルの F 検定
4. 予測式の信頼区間（同時信頼区間）
β̂i
se(
ˆ β̂i )
pi = 2 Pr{T > |ti |}
,
で与えられる．T は自由度 n − p − 1 の t-分布に従う確率変数．確率値は，仮に βi = 0 で
あると仮定したとき，実際に観測した |ti | より絶対値の大きな t-統計量を観測する確率を
表す．これが小さいほど，仮定した βi = 0 と矛盾するので，βi 6= 0 と考えられる．
（背理
回帰係数の t-検定
1
法と同様のロジック．
） 重回帰分析の復習
1.1
• 統計的仮説検定では，あらかじめ閾値 α をさだめておく（通常 α = 0.05 または 0.01 を使
うことが多い）．そして，pi < α ならば βi 6= 0，pi ≥ α ならば βi = 0 と「判定」する．
• 重回帰モデルは
y = β0 + β1 x1 + · · · + βp xp + ²
• もし本当に βi = 0 だとすると，pi は区間 [0, 1] の一様分布に従う（確率変数をその分布関
数で変換してるから）．従って，pi < α となり誤って βi 6= 0 と判定してしまう確率は α
ベクトル表現すると
である．
y = Xβ + ²
• １００％正確な判断はありえない．α は判断の正確さを調整するパラメタと考えてよい．
• 回帰係数ベクトル β を最小二乗法で推定すると
0
α を小さくすればするほど，βi = 0 を誤って βi 6= 0 と判定する確率は小さくなる．しか
0
−1
β̂ = (X X) X y
し逆に，もし βi 6= 0 が真実のときに誤って βi = 0 と誤判定する確率は増える．この両者
はトレードオフの関係にある．
で与えられる．
• 正規モデルを仮定すると，t-統計量は，自由度 n − p − 1 の t-分布に従う確率変数の実現
• 正規モデル
値である．このことを次で示す．
² ∼ Nn (0, σ 2 I n )
を仮定すると，回帰係数ベクトルは多変量正規分布に従う
# run0073.R
β̂ ∼ Np+1 (β, σ 2 (X 0 X)−1 )
# 回帰係数の t-検定
# あらかじめ x に説明変数の行列，y に目的変数のベクトルを設定しておく
ここでサイズ (p + 1) × (p + 1) の行列 A を
fit <- lm(y~.,data.frame(x,y)) # 線形モデルの当てはめ
A = (aij ) = (X 0 X)−1 ,
be <- coef(summary(fit)) # 結果の取り出し
i, j = 0, . . . , p
cat("# 回帰係数，標準誤差，t-統計量，p-値\n"); print(be)
とおけば，
tval <- be[,"t value"] # t-統計量
β̂i ∼ N (βi , σ 2 aii )
mx <- max(abs(tval))+0.5; x0 <- seq(-mx,mx,len=300) # プロットの範囲を決める
• データから V (β̂i ) = σ 2 aii を推定するには，誤差分散 σ̂ 2 をその推定量で置き換える．誤差
² の分散 σ 2 の不偏推定は
S²2 =
1
n−p−1
n
X
plot(x0,dt(x0,df.residual(fit)),type="l",
xlab="t-statistic",ylab="density") # t-分布の密度関数
rug(tval); text(tval,0.01,names(tval),srt=90,adj=0)
e2i =
i=1
kek2
n−p−1
dev.copy2eps(file="run0073-d.eps")
pv0 <- 2*pt(abs(x0),df.residual(fit),lower.tail=F)
回帰係数 β0 , β1 , . . . , βp の個数は p+1 個であることに注意．したがって，自由度は n−(p+1)．
1
2
• 一般に n 個の独立な正規分布 Z1 , . . . , Zn ∼ N (0, 1) の二乗和は自由度 n のカイ二乗分布に
abline(h=c(0.01,0.05),lty=3)
従う
rug(tval); text(tval,pv0[1]*1.1,names(tval),srt=90,adj=0)
n
X
dev.copy2eps(file="run0073-t.eps")
i=1
従って，誤差の二乗和
> dat <- read.table("dat0002.txt") # データの読み込み (47 x 10 行列)
> source("run0073.R")
である．残差二乗和のほうでは，ei ∼ N (0, σ 2 (1 − hii )) だったので，大雑把にいって ei /σ
Pn e2
は近似的に N (0, 1) であるから， i=1 σi2 ∼ χ2n でもよさそうな気もちょっとするが，実際
# 回帰係数，標準誤差，t-統計量，p-値
Std. Error
t value
Pr(>|t|)
(Intercept) 15.89979396 6.030963886
2.6363603 0.012175936
Zouka
2.5211448 0.016136162
1.35299378 0.536658497
Ninzu
-3.26224134 1.066976341 -3.0574636 0.004131840
Kaku
-0.02794050 0.036376197 -0.7680984 0.447303442
Tomo
-0.01586652 0.009506475 -1.6690220 0.103554493
Tandoku
-0.11120432 0.052370545 -2.1234134 0.040470167
X65Sai
0.03597994 0.035336821
Kfufu
には自由度が n − p − 1 になる．
• この事実をキチンと証明するには，まず
e = (I n − H)²
kek2 = ²0 (I n − H)2 ² = ²0 (I n − H)² = k²k2 − ²0 H²
に注意する．H は Span(x0 , . . . , xp ) への射影行列であり，I n − H はその直交補空間への
1.0181996 0.315195078
射影行列である．適当に座標系をとりなおすことにより，n × n の直交行列
-0.20127472 0.058708226 -3.4283904 0.001504385
Ktan
0.10625525 0.053686459
Konin
Zi2 ∼ χ2n
n
X
²2i
∼ χ2n
σ2
i=1
> x <- dat[,-10]; y <- dat[,10]
Estimate
# p-値
plot(x0,pv0,type="l",log="y",xlab="t-statistic",ylab="p-value")
1.9791816 0.055273358
U = (u1 , u2 , . . . , un )
-0.08330936 0.113092312 -0.7366492 0.465981104
を使って，
0.4
H = (u1 , . . . , up+1 )(u1 , . . . , up+1 )0
5e−01
I n − H = (up+2 , . . . , un )(up+2 , . . . , un )0
0.3
とかける．
p−value
−2
0
2
4
−4
t−statistic
−2
0
2
2
kek2 /σ 2 = zp+2
+ · · · + zn2 ∼ χ2n−p−1
Zouka
(Intercept)
Ktan
X65Sai
Kaku
Konin
Tomo
Tandoku
Ninzu
Kfufu
5e−04
Zouka
(Intercept)
Ktan
X65Sai
Kaku
Konin
Tomo
Tandoku
Ninzu
Kfufu
0.1
0.0
−4
i = 1, . . . , n
2
²0 H²/σ 2 = z12 + · · · + zp+1
∼ χ2p+1
5e−03
0.2
density
5e−02
zi = u0i ²/σ,
と置けば，これは互いに独立に N (0, 1) に従う確率変数である．従って，
である．さらに，β̂ は z1 , . . . , zp+1 だけに関係していて zp+2 , . . . , zn とは無関係であるから，
4
t−statistic
run0073-d
run0073-t
β̂ と S²2 = kek2 /(n − p − 1) は互いに独立である．
• ここまでの結果をまとめると，
1. 回帰係数の最小二乗推定 β̂ の従う分布は
1.2
β̂ ∼ Np+1 (β, σ 2 (X 0 X)−1 )
カイ二乗分布と t-分布
• 正規モデルを仮定すると以下の結果が得られる．
ここでサイズ (p + 1) × (p + 1) の行列 A を
A = (aij ) = (X 0 X)−1 ,
2
• 残差二乗和を誤差分散 σ で割ったものは，自由度 n − p − 1 のカイ二乗分布に従う．
n
X
e2i
∼ χ2n−p−1
2
σ
i=1
3
i, j = 0, . . . , p
とおけば，
β̂i ∼ N (βi , σ 2 aii )
4
2. 誤差分散の不偏推定 S²2 の従う分布は
(n − p −
σ2
#n <- 11 # データサイズ
1)S²2
#be <- c(0,1) # 真の回帰係数 (beta0,beta1)
∼ χ2n−p−1
#x <- seq(0,1,len=n) # x を決める
#filename <- "run0074-"
3. β̂ と
S²2
は互いに独立
p <- length(be)-1 # p+1 が回帰係数の個数
• ところで一般に Z ∼ N (0, 1) と S 2 ∼ χ2n が互いに独立なら，これらから作られる確率変
p
数 T = Z/ S 2 /n の従う確率分布は自由度 n の t-分布であることが知られている．
Z
T =p
∼ t-分布n
S 2 /n
確率密度関数は
Γ((n + 1)/2)
f (t) = √
nπΓ(n/2)
µ
¶−(n+1)/2
t2
1+
n
X <- cbind(1,x) # データ行列
colnames(X) <- names(be) <- c("beta0","beta1")
A <- solve(t(X) %*% X) # A=(X’X)^-1
B <- A %*% t(X) # B = (X’X)^-1 X’
sqrA <- sqrt(diag(A)) # A の対角項の平方根
func0074 <- function(y) {
be <- B %*% y # 回帰係数の推定
pred <- X %*% be # 予測値
である．
resid <- y-pred # 残差
se2 <- sum(resid^2)/(n-p-1) # sigma^2 の不偏推定
• これを回帰分析に利用する
se <- sqrt(se2) # sigma の推定
β̂i − βi
∼ N (0, 1)
√
σ aii
tval <- be/(se*sqrA) # t-統計量
pval <- 2*pt(abs(tval),n-p-1,lower.tail=F) # p-値
(n − p − 1)S²2
∼ χ2n−p−1
σ2
なので，
β̂i − βi
√
σ aii
式を整理すると，
.r
list(be=be,se2=se2,tval=tval,pval=pval)
}
S²2
β̂i − βi
= √
∼ t-分布n−p−1
σ2
aii S²
y0 <- X %*% be # 真の y（誤差＝０)
sigma0 <- 0.3 # 誤差の標準偏差の真値
cat("# start simulation: ",date(),"\n")
β̂i − βi
se(
ˆ β̂i )
∼ t-分布n−p−1
b <- 10000 # シミュレーション繰り返し数
simy <- matrix(0,n,b) # y を格納するアレイ
simbe <- matrix(0,length(be),b) # 回帰係数を格納するアレイ
• ところで R で計算している t-統計量は
ti =
β̂i
se(
ˆ β̂i )
であるから，もし βi = 0 ならば
simse2 <- rep(0,b) # se2 を格納するアレイ
simtval <- matrix(0,length(be),b) # t 統計量を格納するアレイ
simpval <- matrix(0,length(be),b) # p-値を格納するアレイ
for(i in 1:b) {
simy[,i] <- y0 + rnorm(n,mean=0,sd=sigma0)
ti ∼ t-分布n−p−1
fit <- func0074(simy[,i])
simbe[,i] <- fit$be; simse2[i] <- fit$se2
である．もし βi > 0 ならば，t-分布から予想されるよりも大きな ti が得られる可能性が
simtval[,i] <- fit$tval; simpval[,i] <- fit$pval
高く，逆にもし βi < 0 ならば，t-分布から予想されるよりも小さな ti が得られる可能性
}
が高い．
cat("# end simulation: ",date(),"\n")
cat("# １回目のシミュレーション結果\n"); print(func0074(simy[,1]))
シミュレーションで確認
1.3
cat("# pval0 < 0.05 の回数 = ",sum(simpval[1,]<0.05),"\n")
cat("# pval1 < 0.05 の回数 = ",sum(simpval[2,]<0.05),"\n")
# run0074.R
if(!is.null(filename)) {
# 回帰係数の分布：シミュレーション
plot(x,y0); abline(be)
source("run0044.R") # drawhist のロード
dev.copy2eps(file=paste(filename,"s0.eps",sep=""))
5
6
plot(x,simy[,1]); abline(simbe[,1])
1.0
dev.copy2eps(file=paste(filename,"s1.eps",sep=""))
1.0
plot(simbe[1,],simbe[2,],pch=".",xlab="beta0",ylab="beta1")
0.8
dev.copy2eps(file=paste(filename,"th1.eps",sep=""))
drawhist(simtval[i,],30,paste("tval",i-1,sep=""))
simy[, 1]
0.4
0.4
y0
for(i in 1:2) {
0.6
dev.copy2eps(file=paste(filename,"th2.eps",sep=""))
0.6
0.8
plot(simse2,simbe[2,],pch=".",xlab="se2",ylab="beta1")
0.2
0.2
t0 <- seq(min(simtval[i,]),max(simtval[i,]),len=300)
lines(t0,dt(t0,n-p-1),col=4,lty=2)
0.0
0.0
dev.copy2eps(file=paste(filename,"tval",i-1,".eps",sep=""))
drawhist(simpval[i,],20,paste("pval",i-1,sep=""),filename)
0.0
0.2
0.4
}
0.6
0.8
1.0
0.0
0.2
0.4
0.6
x
}
0.8
1.0
x
run0074-s0
run0074-s1
2.0
> be <- c(0,1) # 真の回帰係数 (beta0,beta1)
2.0
> n <- 11 # データサイズ
> filename <- "run0074-"
1.5
1.5
> x <- seq(0,1,len=n) # x を決める
1.0
beta1
6 11:45:50 2004
1.0
Wed Oct
Wed Oct
6 11:45:52 2004
# １回目のシミュレーション結果
0.5
# end simulation:
0.5
# start simulation:
beta1
> source("run0074.R")
[,1]
beta0 -0.01592933
beta1
0.0
0.0
$be
−0.6
−0.4
−0.2
0.0
0.2
0.4
0.6
0.00
0.10
0.15
0.35
0.30
0.3
0.25
[,1]
beta0 -0.1181108
0.30
tval1
0.4
$tval
0.25
run0074-th2
tval0
[1] 0.05716638
0.20
se2
run0074-th1
$se2
0.20
4.3051007
0.2
$pval
0.15
beta1
0.05
beta0
0.98142543
[,1]
0.1
0.10
beta0 0.908573939
504
# pval1 < 0.05 の回数 =
8748
0.0
# pval0 < 0.05 の回数 =
0.00
0.05
beta1 0.001975822
−5
0
5
0
mean= −0.0082902 , sd= 1.1296
run0074-tval0
7
5
10
mean= 3.83 , sd= 1.5606
run0074-tval1
8
15
pval0
2.2
pval1
確率値
v̂ ∼ N (v, σv2 )
15
1.0
• 結局
0.8
ただし
0.6
10
v = w0 β,
2
0.4
• σ を
S²2
σv2 = σ 2 w0 (X 0 X)−1 w
で推定すると，
0.2
5
σ̂v2 = S²2 w0 (X 0 X)−1 w = S²2
0
0.0
つまり
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
mean= 0.50308 , sd= 0.28884
run0074-pval0
0.4
0.6
0.8
(n − p − 1)
1.0
mean= 0.024496 , sd= 0.058195
run0074-pval1
• 以上より，
• b = 10000 回のシミュレーション．デスクトップＰＣで計算して 2 秒程度．Pentium 1G
式を整理すると，
のノートパソコン（vmware 上の linux 環境) で実行しても，３秒程度．
2.1
s
σ̂v2
∼ χ2n−p−1
σv2
σ̂v2
∼ t-分布n−p−1
σv2
v̂ − v
∼ t-分布n−p−1
σ̂v
• R では自由度 n − p − 1 の t-分布に従う確率変数 T が t 以下の値をとる確率は
回帰係数，予測値の信頼区間
2
v̂ − v .
σv
σv2
kek2 /σ 2 2
=
σ
σ2
n−p−1 v
Pr{T ≤ t} = pt(t,n-p-1)
回帰係数の線形結合
とかける．これを pt(t) とあらわしておく．
• 回帰係数の線形結合を考え，その性質を調べる
Pr{T > t} = 1-pt(t,n-p-1) = pt(t,n-p-1,lower.tail=F)
v = w0 β0 + w1 β1 + · · · + wp βp
というオプションも用意されていて，確率値の計算では良く用いる．
v̂ = w0 β̂0 + w1 β̂1 + · · · + wp β̂p
• 回帰係数の t-検定では，βi = 0 という仮説を考えた．より一般的に，帰無仮説 v = v0 を
• w = (w0 , w1 , . . . , wp ) を設定することにより，
wj = 1, wk = 0, k 6= j
⇒
v̂ = β̂j
w0 = xi0 , . . . wp = xip
⇒
v̂ = ŷi
対立仮説 v 6= v0 に対して検定することを考える．p-値は
¯
¯¾
¯¶¶
½
µ
µ¯
¯ v̂ − v0 ¯
¯
¯
¯ = 2 × 1 − pt ¯ v̂ − v0 ¯
p-値 (v0 ) = Pr |T | > ¯¯
¯ σ̂v ¯
σ̂v ¯
ここで帰無仮説 v = v0 の確率値を p-値 (v0 ) とあらわした．これがもし p-値 (v0 ) < α なら
などが表現できるので便利．
ば，この v = v0 という仮説が棄却され v 6= v0 と判断される．
• ベクトル表現
v = w0 β,
• 特に真値 v に関しては，
v̂ = w0 β̂
• 回帰係数は多変量正規分布に従う
v̂ − v .
σv
より
β̂ ∼ Np+1 (β, σ 2 (X 0 X)−1 )
s
σ̂v2
∼ t-分布n−p−1
σv2
p-値 (v) ∼ 一様分布 [0, 1]
となり，
従って，v̂ は正規分布に従う
Pr{p-値 (v) < α} = α
v̂ ∼ N (w0 β, σ 2 w0 (X 0 X)−1 w)
である．
9
2.3
10
信頼区間
2.4
• 先ほど，帰無仮説 v = v0 の確率値を p-値 (v0 ) とあらわした．これがもし < α ならば，こ
回帰係数の信頼区間
• 回帰係数 βj は w を次のようにすればよい．
の v = v0 という仮説が棄却され v 6= v0 と判断される．
wj = 1, wk = 0, k 6= j
• この方法で棄却されない v0 の集合を考える
• V (v̂) の推定は
信頼区間 (1 − α) = {v0 | p-値 (v0 ) ≥ α}
これを信頼係数（または信頼度）1 − α の信頼区間と呼ぶ．たとえば α = 0.05 ならば，信
⇒
σ̂v2 = S²2 w0 (X 0 X)−1 w = S²2 ajj
• βj の信頼区間は
頼係数 0.95 の信頼区間である．
√
信頼区間 (1 − α) = [β̂j − S² ajj qt(1 − α/2),
• 真値 v が信頼区間に入る確率は 1 − α（以上）である．
Pr{v ∈ 信頼区間 (1 − α)} = 1 − α
v̂ = β̂j
2.5
（証明）「v ∈ 信頼区間 (1 − α)」と「p-値 (v) ≥ α」は同値である．ところが p-値 (v) は区
√
β̂j + S² ajj qt(1 − α/2)]
予測値の信頼区間
• 説明変数が x = (1, x1 , . . . , xp )0 の時の予測値は ŷ = x0 β̂ である．この「真値」は y = x0 β
間 [0,1] の一様分布に従うから，この事象が起こる確率は 1 − α である．
である．v = y とするには w = x とおけばよい．
• 信頼区間を計算するために，まず a = pt(b) の逆関数として，b = qt(a) を用意する．つま
• V (v̂) の推定は
り，自由度 n − p − 1 の t-分布に従う確率変数 T が t 以下の値をとる確率がちょうど a に
σ̂v2 = S²2 w0 (X 0 X)−1 w = S²2 x0 (X 0 X)−1 x
なるような t の値を qt(a) と書く．
• y の信頼区間は
Pr{T ≤ qt(a)} = a
信頼区間 (1−α) = [x0 β̂−S²
R では
q
x0 (X 0 X)−1 x qt(1−α/2), x0 β̂+S²
q
x0 (X 0 X)−1 x qt(1−α/2)]
qt(a) = qt(a,n-p-1)
# run0075.R
という関数が用意されている．0 < a < 1 である．
• 信頼区間は
# 回帰係数と予測値の信頼区間：多項式回帰
¯¶¶
µ¯
¯ v̂ − v0 ¯
¯
2 × 1 − pt ¯¯
≥α
σ̂v ¯
µ
を整理して
¯¶
µ¯
¯ v̂ − v0 ¯
¯ ≤ 1 − α/2
pt ¯¯
σ̂v ¯
従って，
これを整理すると，
¯
¯
¯ v̂ − v0 ¯
¯
¯
¯ σ̂v ¯ ≤ qt(1 − α/2)
# x=説明変数のベクトル，y=目的変数のベクトル，p=多項式の次数
# a <- func0075a(x,p); b <- func0075b(y,0.05,a)
xpow <- function(a,p) a^(0:p) # c(a^0,a^1,...,a^p)
calcX <- function(x,p) { # デザイン行列 X を作る
X <- matrix(0,length(x),p+1)
for(i in 1:length(x)) X[i,] <- xpow(x[i],p)
X
}
calcq0 <- function(n,p,alpha) qt(1-alpha/2,n-p-1) # 個別の信頼区間
v̂ − σ̂v qt(1 − α/2) ≤ v0 ≤ v̂ + σ̂v qt(1 − α/2)
calcq1 <- function(n,p,alpha) sqrt((p+1)*qf(1-alpha,p+1,n-p-1)) # 同時信頼区間
func0075a <- function(x,p) { # 回帰分析の準備
X <- calcX(x,p) # データ行列
つまり
信頼区間 (1 − α) = [v̂ − σ̂v qt(1 − α/2),
v̂ + σ̂v qt(1 − α/2)]
colnames(X) <- paste("beta",0:p,sep="")
A <- solve(t(X) %*% X) # A=(X’X)^-1
B <- A %*% t(X) # B = (X’X)^-1 X’
sqrA <- sqrt(diag(A)) # A の対角項の平方根
x0 <- seq(min(x),max(x),len=300) # x の範囲を 300 等分しておく
11
12
X0 <- calcX(x0,p) # x0 に相当するデータ行列
a <- func0075a(x,i)
sqrXAX <- apply(X0,1,function(x) sqrt(t(x) %*% A %*% x))
c1 <- func0075c(x,y,a,func0075b(y,0.05,a))
# sqrt(x’Ax) のベクトル
c2 <- func0075c(x,y,a,func0075b(y,0.01,a),col=3,add=T)
list(X=X,A=A,B=B,sqrA=sqrA,x0=x0,X0=X0,sqrXAX=sqrXAX)
colnames(c1$coef)[5:6] <- c("Lo05","Up05")
}
colnames(c2$coef)[5:6] <- c("Lo01","Up01")
func0075b <- function(y,alpha,a,calcq=calcq0) { # 信頼区間の計算
coef <- cbind(c1$coef,c2$coef[,5:6,drop=F])
n <- nrow(a$X); p <- ncol(a$X)-1
cat("RSS=",c1$rss,", RSQ=",c1$rsq,"\n")
q0 <- calcq(n,p,alpha)
print(round(coef,3))
be <- a$B %*% y # 回帰係数の推定
dev.copy2eps(file=paste("run0076-s",i,".eps",sep=""))
pred <- a$X %*% be # 予測値
}
resid <- y-pred # 残差
rss <- sum(resid^2) # 残差平方和
> dat <- read.table("dat0001.txt") # データの読み込み (47 x 2 行列)
se2 <- rss/(n-p-1) # sigma^2 の不偏推定
> x <- dat[,1]/100; y <- dat[,2]
se <- sqrt(se2) # sigma の推定
> p <- 3
bese <- se*a$sqrA # 回帰係数の標準誤差
> source("run0076.R")
rsq <- 1-rss/sum((y-mean(y))^2) # 決定係数
# 次数= 0
tval <- be/(se*a$sqrA) # t-統計量
RSS= 0.815383 , RSQ= 0
pval <- 2*pt(abs(tval),n-p-1,lower.tail=F) # p-値
Estimate StdErr t-value p-value
beconf <- cbind(be-q0*se*a$sqrA,be+q0*se*a$sqrA) # 信頼区間
beta0
pred0 <- a$X0 %*% be # 予測値 (x0)
# 次数= 1
pred0conf <- cbind(pred0-q0*se*a$sqrXAX,pred0+q0*se*a$sqrXAX)
RSS= 0.3812667 , RSQ= 0.5324078
list(be=be,bese=bese,se=se,rss=rss,rsq=rsq,tval=tval,pval=pval,
1.473
0.019
75.848
Estimate StdErr t-value p-value
beconf=beconf,pred0=pred0,pred0conf=pred0conf)
Lo05
Up05
Lo01
Up01
0 1.434 1.512 1.421 1.525
Lo05
Up05
Lo01
Up01
1.662
1.823
1.635
1.850
beta0
1.742
0.040
43.592
0
}
beta1
-2.825
0.395
-7.158
0 -3.620 -2.030 -3.886 -1.763
func0075c <- function(x,y,a,b,col=2,lty=2,add=F) {
# 次数= 2
if(!add) plot(x,y)
RSS= 0.3786836 , RSQ= 0.5355758
lines(a$x0,b$pred0)
Lo05
Up05
Lo01
Up01
lines(a$x0,b$pred0conf[,1],col=col,lty=lty)
beta0
Estimate StdErr t-value p-value
1.681
0.120
14.043
0.000
1.440
1.922
1.359
2.003
lines(a$x0,b$pred0conf[,2],col=col,lty=lty)
beta1
-1.672
2.141
-0.781
0.439
-5.987
2.642
-7.436
4.091
coef <- cbind(b$be,b$bese,b$tval,b$pval,b$beconf)
beta2
-4.698
8.576
-0.548
0.587 -21.983 12.586 -27.788 18.391
colnames(coef) <- c("Estimate","StdErr",
# 次数= 3
"t-value","p-value","Lower","Upper")
RSS= 0.3747561 , RSQ= 0.5403926
invisible(list(coef=coef,se=b$se,rss=b$rss,rsq=b$rsq))
Estimate
}
beta0
1.462
StdErr t-value p-value
0.347
Lo05
Up05
Lo01
Up01
4.212
0.000
0.762
2.162
0.527
2.398
0.638
-14.381
23.211
-20.704
29.534
beta1
4.415
9.320
0.474
# run0076.R
beta2
-56.680
77.913
-0.727
0.471 -213.807 100.447 -266.664 153.304
# 回帰係数と予測値の信頼区間：多項式回帰 p 次まで
beta3
135.499 201.844
0.671
0.506 -271.559 542.557 -408.492 679.491
# x=説明変数ベクトル，y=目的変数ベクトル，p=多項式次数をセットしておく
source("run0075.R")
for(i in 0:p) {
cat("# 次数=",i,"\n")
13
14
3. 次数 p = 2 のとき，β0 = 0 を検定する確率値は，ほぼゼロ．つまり β0 6= 0 と結論．
1.8
結論．これは p = 0 を示唆？ 矛盾？
4. 次数 p = 3 のときは，β0 6= 0，β1 = β2 = β3 = 0 と結論．これは p = 0 を示唆？ 1.6
矛盾？
y
• この例からも分かるように，回帰係数の t-検定は解釈に注意する．
1.4
1.4
y
1.6
1.8
ところが β1 = 0 を検定する確率値は 0.439 なので β1 = 0 と結論．同様に β2 = 0 と
• この場合は，p = 1 と判断するのが適切．次数 p = 2 のチェックのときは，β2 = 0 or β2 6= 0
1.2
1.2
だけをチェックすべき．
0.05
0.10
0.15
0.20
0.05
0.10
x
0.15
0.20
回帰モデルの F 検定
3
x
run0076-s0
3.1
run0076-s1
部分モデル
1.8
y = β0 + β1 x1 + · · · + βp xp + ²
ベクトル表現すると
1.6
1.6
1.8
• これまで p 個の説明変数の重回帰モデルを考えた
y
1.4
1.4
y
y = Xβ + ²
• 最初の k 個の説明変数だけを用いる重回帰モデルを考える（例：多項式回帰で次数を p か
y = β0 + β1 x1 + · · · + βk xk + ²
1.2
1.2
ら k に変更する）
ベクトル表現すると
y = X (k) β (k) + ²
0.05
0.10
0.15
0.20
0.05
0.10
x
0.15
0.20
x
run0076-s2
ただし
X (k) = (x0 , x1 , . . . , xk ),
run0076-s3
• 各次数 p = 0, 1, 2, 3 のグラフを描いた．予測値とその信頼区間も示した．赤は α = 0.05，
• なお，
緑は α = 0.01 である．
• 次数を上げると RSS = kek2（残差平方和）は小さくなる．決定係数 R2 は大きくなる．い
X (−k) = (xk+1 , xk+2 , . . . , xp ),
るのが適切なのか？
• 回帰係数の t 検定の確率値を順番にみていく
1. 次数 p = 0 のとき，β0 = 0 を検定する確率値は，ほぼゼロ．つまり β0 6= 0 と結論．
つまり，p ≥ 0 が示唆される．
X = (X (k) , X (−k) ),
β=
"
β (k)
β
(−k)
#
と分割してかける．
• 回帰係数 β (k) の最小二乗推定は，
β̂
i = 0, . . . , k について，一般に β̂
(k)
(k)
0
0
= (X (k) X (k) )−1 X (k) y
の第 i 成分と β̂ の第 i 成分は一致しない．
2. 次数 p = 1 のとき，β0 = 0 を検定する確率値は，ほぼゼロ．つまり β0 6= 0 と結論．
同様に β1 6= 0 と結論．つまり，p ≥ 1 が示唆される．
15
β (−k) = (βk+1 , βk+2 , . . . , βp )0
とおけば
ずれも，次数の増加とともに，回帰の当てはまりがよくなることを示唆する．
• ところがグラフを見ると，次数 p = 1 くらいで十分な感じ．いったい，次数はいくつにす
β (k) = (β0 , β1 , . . . , βk )0
16
• p 個の説明変数をもつモデルにおいて，
• F 分布の分布関数は R の関数 pf で計算できる．一般に X が自由度 m, n の F 分布に従う
とき，この分布関数は
βk+1 = βk+2 = · · · = βp = 0
Pr{X ≤ x} = pf(x, m, n)
と設定すれば，k 個の説明変数をもつモデルになる．したがって，後者は前者の「部分モ
Pr{X > x} = 1 − pf(x, m, n) = pf(x, m, n, lower.tail = F)
デル」または「部分回帰」と呼ばれる．
pf の逆関数は qf である．
• ここでは最初の k 個の説明変数としたが，添え字を付け替えることにより，実際にはどの
Pr{X ≤ qf(a, m, n)} = a
k 個を選んでも，同様な議論ができる．
である．
• さらに以降の議論で本質的なのは，
Span(x0 , . . . , xk ) ⊂ Span(x0 , . . . , xp )
ということなので，説明変数の線形結合をとっても良い．
# run0077.R
# Ｆ検定
fkentei <- function(fitk,fitp) {
rssk <- sum(resid(fitk)^2) # RSS^(k)
3.2
残差平方和と F 検定
rssp <- sum(resid(fitp)^2) # RSS^(p)
dfk <- df.residual(fitk) # n-k-1
• 回帰モデル (k) のハット行列
H
(k)
=X
(k)
(X
(k)0
X
(k) −1
) X
dfp <- df.residual(fitp) # n-p-1
(k)0
bunshi <- (rssk - rssp)/(dfk-dfp) # (RSS^(k)-RSS^(p))/(p-k)
予測値のベクトル
bunbo <- rssp/dfp # RSS^(p) / (n-p-1)
ŷ (k) = H (k) y
fvalue <- bunshi/bunbo # F 統計量
pvalue <- pf(fvalue,dfk-dfp,dfp,lower.tail=F)
残差ベクトル
e(k) = y − ŷ (k) = (I n − H (k) )y
• 回帰モデル (k) の残差平方和は
RSS(k) = ke(k) k2
である．もしモデル (k) が正しければ，
である．さらに残差平方和の差
RSS
> dat <- read.table("dat0001.txt") # データの読み込み (47 x 2 行列)
> x <- dat[,1]/100; y <- dat[,2]
> fit0 <- lm(y ~ 1, data.frame(x,y)) # k=0
(p)
− RSS
σ2
=
ke(k) k2 ke(p) k2
−
∼ χ2p−k
σ2
σ2
である．これと
RSS(p)
∼ χ2n−p−1
σ2
> fit1 <- lm(y ~ x, data.frame(x,y)) # k=1
> fit2 <- lm(y ~ x + I(x^2), data.frame(x,y)) # k=2
> summary(fit2) # k=2 のサマリ
Call:
は互いに独立である．
• もしモデル (k)
> source("run0077.R")
> # 多項式回帰
RSS(k)
∼ χ2n−k−1
σ2
(k)
list(fvalue=fvalue,df1=dfk-dfp,df2=dfp,pvalue=pvalue)
}
(k)
が正しくない場合， RSS
σ2
• モデル (k) が正しければ，
lm(formula = y ~ x + I(x^2), data = data.frame(x, y))
は「非心カイ二乗分布」という分布に従う．
Residuals:
(RSS(k) − RSS(p) )/(p − k)
F =
RSS(p) /(n − p − 1)
は互いに独立な自由度 p − k のカイ二乗と自由度 n − p − 1 のカイ二乗の比であり，した
がって，自由度 p − k, n − p − 1 の F 分布に従う．
F ∼ Fp−k,n−p−1
Median
3Q
Max
-0.29411 -0.04875 -0.00797
Min
1Q
0.04600
0.32282
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)
1.6807
0.1197
14.043
17
<2e-16 ***
18
x
-1.6725
2.1408
-0.781
0.439
I(x^2)
-4.6985
8.5762
-0.548
0.587
• パラメタの最尤推定を θ̂1 , θ̂2 ，最大対数尤度を `1 (θ̂1 )，`2 (θ̂2 ) と書く．
（重回帰モデルの例）
`1 (θ̂1 ) = `(β̂0 , β̂1 , . . . , β̂k , σ̂) = −
n
2
(
RSS(k)
)
n
`2 (θ̂1 ) = `(β̂0 , β̂1 , . . . , β̂p , σ̂) = −
n
2
(
)
RSS(p)
)
n
--Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.09277 on 44 degrees of freedom
Multiple R-Squared: 0.5356,Adjusted R-squared: 0.5145
F-statistic: 25.37 on 2 and 44 DF,
1 + log(2π
1 + log(2π
)
p-value: 4.7e-08
• モデル 1 が正しいとき，対数尤度の差の２倍は，近似的に カイ二乗分布に従う．自由度
は p2 − p1 である．
（n が十分大きいとき，近似がよくなる）．つまり
> unlist(fkentei(fit0,fit2)) # 上の F 検定はコレ!
fvalue
df1
df2
pvalue
2 × (`2 (θ̂1 ) − `1 (θ̂2 )) ∼ χ2p2 −p1
2.537048e+01 2.000000e+00 4.400000e+01 4.700319e-08
（重回帰モデルの例）
> unlist(fkentei(fit1,fit2)) # 上の t 検定 (x^2 の項) はコレ!
fvalue
0.3001376
df1
df2
pvalue
1.0000000 44.0000000
0.5865649
2 × (`2 (θ̂1 ) − `1 (θ̂2 )) = n log RSS(k) − n log RSS(p) ∼ χ2p−k
> sqrt(0.3001376) # これが t 統計量 (x^2 の項) の絶対値
[1] 0.5478482
# run0081.R
• F 検定は，モデル (k) からモデル (p) に変更したときに，当てはまり（ＲＳＳ値）がどれ
だけ改善したかをチェックしている．
# 尤度比検定
yudohikentei <- function(fitk,fitp) {
rssk <- sum(resid(fitk)^2) # RSS^(k)
• 特によく利用されるのは次の２パターン
rssp <- sum(resid(fitp)^2) # RSS^(p)
dfk <- df.residual(fitk) # n-k-1
1. k = 0 とする．これは y = 定数 というモデルと比較して，モデル (p) の回帰が意味あ
dfp <- df.residual(fitp) # n-p-1
るのかどうかを調べる．これで F 検定の p-value が < α となり有意だとしても，p 個
n <- length(resid(fitk)) # n
の説明変数すべてが必要であるという意味ではない．
（その一部で十分かもしれない．
）
chisqvalue <- n*(log(rssk)-log(rssp)) # 尤度比統計量
2. k = p − 1 とする．これは説明変数 xp が必要かどうかをチェックし，回帰係数 βp = 0
pvalue <- pchisq(chisqvalue,dfk-dfp,lower.tail=F)
or βp 6= 0 を判定する．数学的に βp の t- 検定と等価．(おなじ確率値になる．) また，
list(chisqvalue=chisqvalue,df=dfk-dfp,n=n,pvalue=pvalue)
2
統計量も |t| = F の関係がある．したがって，通常 F 検定ではなくて，t-検定とし
}
て実装されている．
> dat <- read.table("dat0001.txt") # データの読み込み (47 x 2 行列)
3.3
尤度比検定＊
> x <- dat[,1]/100; y <- dat[,2]
• 一般の確率モデルの比較では，F 検定の変わりに尤度比検定が適用できる．
> fit0 <- lm(y ~ 1, data.frame(x,y)) # k=0
> fit1 <- lm(y ~ x, data.frame(x,y)) # k=1
• 尤度比検定を重回帰モデルに適用すると，F 検定とほぼ同じ結果になる．
> fit2 <- lm(y ~ x + I(x^2), data.frame(x,y)) # k=2
• モデル１とモデル２は互いに包含関係（ネスト）であり，モデル２はモデル１を一般化し
> unlist(fkentei(fit0,fit2)) # F 検定
たものとする．モデル１はモデル２の特殊な場合になる．
（重回帰モデルの例）モデル１＝モデル (k)，モデル２＝モデル (p)．
• まず二つの確率モデルの対数尤度を `1 (θ1 )，`2 (θ2 ) とする．ここで θ1 と θ2 はパラメタベ
クトルで，次元はそれぞれ p1 = dim θ1 ,p2 = dim θ2 とおく．
> source("run0081.R")
fvalue
df1
df2
pvalue
2.537048e+01 2.000000e+00 4.400000e+01 4.700319e-08
> unlist(yudohikentei(fit0,fit2)) # 尤度比検定
chisqvalue
df
n
pvalue
（重回帰モデルの例）θ1 = (β0 , β1 , . . . , βk , σ)，θ2 = (β0 , β1 , . . . , βp , σ)，p1 = k + 2，p2 =
3.604697e+01 2.000000e+00 4.700000e+01 1.487646e-08
p + 2．
> unlist(fkentei(fit1,fit2)) # F 検定
19
20
fvalue
0.3001376
df2
pvalue
[1,] -0.4472136 -0.4262868
0.4925791
1.0000000 44.0000000
df1
0.5865649
[2,] -0.4472136 -0.3443086
0.1516455
> unlist(yudohikentei(fit1,fit2)) # 尤度比検定
[3,] -0.4472136 -0.1803521 -0.3301785
chisqvalue
0.3195130
n
pvalue
[4,] -0.4472136
0.1475608 -0.6936976
1.0000000 47.0000000
df
0.5719005
[5,] -0.4472136
0.8033866
0.3796515
> Q %*% R
3.4
[,1] [,2] [,3]
QR 分解
• 最小二乗法の計算では，通常，QR 分解と呼ばれる行列の分解をしている．
X = QR
ここで，X はサイズ n × (p + 1)，Q はサイズ n × (p + 1)，R はサイズ (p + 1) × (p + 1)．
Q の列は互いに直交していて，Q0 Q = I p+1 である．R は上三角行列である．
[1,]
1
2
[2,]
1
4
3
9
[3,]
1
8
27
[4,]
1
16
81
[5,]
1
32
243
> t(Q) %*% Q
[,1]
X 0 X = R0 Q0 QR = R0 R
(X 0 X)−1 = (R0 R)−1 = R−1 R−1
[1,]
0
[,2]
[,3]
1.000000e-00 -8.180305e-17
1.924459e-17
[2,] -8.180305e-17
1.000000e-00 -1.428436e-17
[3,]
1.924459e-17 -1.428436e-17
3.5
F 検定のつづき＊
1.000000e-00
• QR 分解は，β̂ の計算時に利用される．
0
β̂ = (X 0 X)−1 X 0 y = R−1 R−1 R0 Q0 y = R−1 Q0 y
まず Q0 y を計算した後，β̂ = R−1 (Q0 y) の計算時には R−1 を計算する必要はない．R が
• まえに議論した直交基底
U = (u1 , u2 , . . . , un )
上三角であることを利用して，Q0 y と R から直接 β̂ を計算するアルゴリズムがあり，こ
のほうが R−1 を経由するより演算数が少ないし，数値的にも安定する．
をうまく定義すると，
H (k) = (u1 , . . . , uk+1 )(u1 , . . . , uk+1 )0 ,
> X <- matrix(c(1^(1:5),2^(1:5),3^(1:5)),5,3)
> X
k = 0, . . . , p
とできる．x0 , x1 , . . . , xp を順番にグラム・シュミットの直交化すればよい．じつは QR
[,1] [,2] [,3]
分解 X = QR は，これと同等の計算をしている（ここでは述べないが，グラム・シュ
[1,]
1
2
3
ミットの直交化よりも効率の良いアルゴリズム）．ui = q i , i = 1, . . . , p + 1，のこりの
[2,]
1
4
9
up+2 , . . . , un はこれらと直交する適当な基底を選べばよい．
[3,]
1
8
27
[4,]
1
16
81
[5,]
1
32
243
• 回帰モデル (p) で正規モデルを仮定する．
y = Xβ + ²
> QR <- qr(X)
² ∼ Nn (0, σ 2 I n )
> R <- qr.R(QR)
これは
> Q <- qr.Q(QR)
y = X (k) β (k) + X (−k) β (−k) + ²
> R
[,1]
[,2]
[,3]
ともかけるので，回帰モデル (k) の残差ベクトルは
[1,] -2.236068 -27.72724 -162.33854
[2,]
0.000000
24.39672
197.92824
[3,]
0.000000
0.00000
29.99355
e(k) = (I n − H (k) )y = (I n − H (k) )(X (−k) β (−k) + ²)
もし β
> Q
[,1]
[,2]
(−k)
= 0 ならば
(k)
e
[,3]
= (I n − H (k) )y = (I n − H (k) )² = (uk+2 , . . . , un )(uk+2 , . . . , un )0 ²
22
21
したがって，zi = u0 ²/σ, i = 1, . . . , n とおけば，
• 以上より，
ke(k) k2
2
= zk+2
+ · · · + zn2
σ2
F =
kzk2 /(p + 1)
∼ Fp+1,n−p−1
S²2 /σ 2
は自由度 p + 1, n − p − 1 の F 分布に従う．kzk2 = kX(β̂ − β)k2 /σ 2 を代入すると，
• 結局，モデル (k) が正しいとき，
RSS(k)
ke(k) k2
2
=
= zk+2
+ · · · + zn2 ∼ χ2n−k−1
σ2
σ2
RSS(p)
ke(p) k2
2
=
= zp+2
+ · · · + zn2 ∼ χ2n−p−1
σ2
σ2
(k)
(p)
RSS − RSS
ke(k) k2 ke(p) k2
2
2
=
−
= zk+2
+ · · · + zp+1
∼ χ2p−k
σ2
σ2
σ2
(p)
(k)
(p)
そして，RSS と RSS − RSS は，互いに共通の zi を含まないので，独立．
F =
kX(β̂ − β)k2
(p + 1)S²2
である．
• β̂ の信頼域（信頼係数 or 信頼度＝ 1 − α）を作るには次のように考える．
qf(a) = qf(a, p + 1, n − p − 1)
と置けば，
Pr{F ≤ qf(1 − α)} = 1 − α
予測式の信頼区間（同時信頼区間）
4
4.1
である．そこで
n
信頼域 (1 − α) = β
回帰係数ベクトルの同時信頼域
• 復習
(β − β̂)0 (X 0 X)(β − β̂) ≤ (p + 1)S²2 qf(1 − α)
:
o
と定めれば，これは β̂ を中心とする楕円体である．そして
1. 正規回帰モデルを仮定する．回帰係数の最小二乗推定 β̂ は多変量正規分布に従う．
Pr{β ∈ 信頼域 (1 − α)} = 1 − α
β̂ ∼ Np+1 (β, σ 2 (X 0 X)−1 )
である．
2. 誤差分散 σ 2 の不偏推定 S²2 は
(n − p − 1)S²2
∼ χ2n−p−1
σ2
• 信頼域を実際にグラフ描いてみるには，
γ = Rβ,
であり，これは β̂ とは独立である．
• ここで，適当な (p + 1) × (p + 1) 行列 R を選ぶと，
戻せばよい．γ の信頼域は
X 0 X = R0 R
©
とできることに注意する．具体的には，前に説明した QR 分解 X = QR から得られる R
でよい．
(X 0 X)−1 = (R0 R)−1 = R−1 R−1
0
• p + 1 ベクトル z = (z1 , z2 , . . . , zp+1 )0 を
γ
:
ª
kγ − γ̂k2 ≤ (p + 1)S²2 qf(1 − α)
であるから γ̂ を中心とする半径 S²
p
(p + 1)qf(1 − α) の球体．
# run0078.R
1
R(β̂ − β)
σ
で定義する．z は多変量正規分布に従い，平均は E(z) = 0，分散は
z=
V (z) =
γ̂ = Rβ̂
と回帰係数を変数変換して γ のほうで考えたほうが分かりやすい．それを β = R−1 γ で
1
1
0
Rσ 2 (X 0 X)−1 R0 = RR−1 R−1 R0 = I p+1
σ
σ
従って，
# 回帰係数の信頼域：単回帰
# x=説明変数のベクトル，y=目的変数のベクトル
source("run0075.R")
p <- 1 # 単回帰なので次数=1
n <- length(y) # サンプルサイズ
alpha <- 0.05 # 信頼係数=1-alpha
z ∼ Np+1 (0, I p+1 )
成分で書けば，
a <- func0075a(x,p) # 回帰分析の準備
a$QR <- qr(a$X); a$R <- qr.R(a$QR) # QR 分解
z1 , z2 , . . . , zp+1 ∼ N (0, 1)
でこれらは互いに独立．
23
a$IR <- solve(a$R) # R^-1
b0 <- func0075b(y,alpha,a) # 回帰係数など
24
Pr{β1 ∈ 信頼区間1 } = 1 − α
th <- seq(0,2*pi,len=300) # 描画の範囲（0..2*pi を 300 等分）
ga <-
b0$se*calcq1(n,p,alpha)*rbind(cos(th),sin(th)) # gamma でみた円
であるが，一般に
be <- as.vector(b0$be) + a$IR %*% ga # beta の信頼域（楕円）
Pr{β0 ∈ 信頼区間0 ,
plot(t(be),type="l") # 楕円
points(t(b0$be)) # 中心
β1 ∈ 信頼区間1 } < 1 − α
となる可能性がある．
abline(v=b0$beconf[1,],lty=2,col=2) # beta0 の信頼区間（赤）
abline(h=b0$beconf[2,],lty=2,col=2) # beta1 の信頼区間（赤）
• 各係数について，同時信頼域の最小値，最大値を直線（青）で示した．これを「同時信頼
b1 <- func0075b(y,alpha,a,calcq1) # calcq1 を使ってもう一度
区間」と呼ぶと
abline(v=b1$beconf[1,],lty=3,col=4) # beta0 の同時信頼区間（青）
Pr{β0 ∈ 同時信頼区間0 ,
abline(h=b1$beconf[2,],lty=3,col=4) # beta1 の同時信頼区間（青）
dev.copy2eps(file="run0078-b.eps")
β1 ∈ 同時信頼区間1 } ≥ 1 − α
である．β の集合としてみると，同時信頼域より大きくなってしまっている．
coef <- cbind(b0$be,b0$beconf,b1$beconf)
colnames(coef) <- c("Estimate","Lo","Up","JointLo","JointUp")
print(coef)
• 「回帰式の信頼区間」，言い換えると，
「予測値の同時信頼区間」を求める． 一般には，回
> dat <- read.table("dat0001.txt") # データの読み込み (47 x 2 行列)
帰係数の線形結合 v = w0 β の同時信頼区間を考える．
> x <- dat[,1]/100; y <- dat[,2]
• 同時ではない，個別の信頼区間は以前求めた．
> source("run0078.R")
beta0
回帰式の信頼区間
4.2
Estimate
Lo
Up
JointLo
JointUp
1.742483
1.661974
1.822993
1.641291
1.843676
信頼区間 (1 − α) = [v̂ − σ̂v qt(1 − α/2),
beta1 -2.824869 -3.619719 -2.030019 -3.823917 -1.825821
ただし
σ̂v = S²
q
v̂ + σ̂v qt(1 − α/2)]
w0 (X 0 X)−1 w
−2.0
qt(a) = qt(a,n-p-1)
この信頼区間は
−2.5
Pr{v ∈ 信頼区間 (1 − α)} = 1 − α
−3.0
beta1
を満たす．
• ここで求めたい「同時信頼区間」は，w に関する任意の集合 W に対して，
−3.5
Pr{w0 β ∈ 同時信頼区間w(1 − α),
w ∈ W} ≥ 1 − α
となるようなものである．
1.65
1.70
1.75
1.80
1.85
beta0
• (例) 前節で求めた，回帰係数 β0 と β1 は集合 W の要素が二つの場合に相当する．そこで
run0078-b
求めた同時信頼区間は
• 単回帰の例：β = (β0 , β1 )0 の同時信頼域（楕円）を描いた．以下すべて α = 0.05 とした．
Pr{β0 ∈ 同時信頼区間0 ,
Pr{β ∈ 同時信頼域 } = 1 − α
β1 ∈ 同時信頼区間1 } ≥ 1 − α
を満たしていた．
• 各回帰係数 βi , i = 0, 1 の信頼区間を直線で示した（赤）．これは各係数を別々に見て信頼
区間を t-検定から得ている．つまり，
Pr{β0 ∈ 信頼区間0 } = 1 − α
26
25
• 導出は後ほど述べるが，結論としては，次のように計算すればよい．
同時信頼区間w(1 − α) = [v̂ − σ̂v
ただし
p
(p + 1)qf(1 − α),
v̂ + σ̂v
p
(p + 1)qf(1 − α)]
qf(a) = qf(a, p + 1, n − p − 1)
である．
Estimate StdErr t-value p-value
beta0
1.473
0.019
75.848
Lo
Up JointLo JointUp
0 1.434 1.512
1.434
1.512
# 次数= 1
RSS= 0.3812667 , RSQ= 0.5324078
Estimate StdErr t-value p-value
Lo
Up JointLo JointUp
beta0
1.742
0.040
43.592
0
1.662
1.823
1.641
1.844
beta1
-2.825
0.395
-7.158
0 -3.620 -2.030
-3.824
-1.826
# 次数= 2
# run0079.R
RSS= 0.3786836 , RSQ= 0.5355758
# 回帰係数と予測値の同時信頼区間：多項式回帰
Estimate StdErr t-value p-value
Lo
Up JointLo JointUp
# x=説明変数のベクトル，y=目的変数のベクトル，p=多項式の次数, alpha も．
beta0
1.681
0.120
14.043
0.000
1.440
1.922
1.333
source("run0075.R")
beta1
-1.672
2.141
-0.781
0.439
-5.987
2.642
-7.895
4.550
func0079 <- function(x,y,a,alpha,output=T) {
beta2
-4.698
8.576
-0.548
0.587 -21.983 12.586 -29.628
20.231
b0 <- func0075b(y,alpha,a)
# 次数= 3
b1 <- func0075b(y,alpha,a,calcq1)
RSS= 0.3747561 , RSQ= 0.5403926
coef <- cbind(b0$be,b0$bese,b0$tval,b0$pval,b0$beconf,b1$beconf)
colnames(coef) <- c("Estimate","StdErr","t-value","p-value",
"Lo","Up","JointLo","JointUp")
if(output) {
func0075c(x,y,a,b0); func0075c(x,y,a,b1,col=4,lty=3,add=T)
Estimate
StdErr t-value p-value
2.029
Lo
Up
beta0
1.462
0.347
4.212
0.000
0.762
2.162
JointLo JointUp
0.345
2.579
beta1
4.415
9.320
0.474
0.638
-14.381
23.211
-25.578
34.407
beta2
-56.680
77.913
-0.727
0.471 -213.807 100.447 -307.403 194.042
beta3
135.499 201.844
0.671
0.506 -271.559 542.557 -514.031 785.029
cat("RSS=",b0$rss,", RSQ=",b0$rsq,"\n")
1.8
1.8
print(round(coef,3))
}
list(b0=b0,b1=b1,coef=coef)
func0079b <- function(x,y,p,alpha,filename="run0079-") {
1.6
1.6
}
y
1.4
cat("# 次数=",i,"\n")
1.4
y
for(i in 0:p) {
1.2
func0079(x,y,a,alpha)
1.2
a <- func0075a(x,i)
if(!is.null(filename))
dev.copy2eps(file=paste(filename,"s",i,".eps",sep=""))
}
0.05
0.10
0.15
0.20
0.05
x
}
0.15
x
run0079-s0
run0079-s1
> source("run0079.R")
> dat <- read.table("dat0001.txt") # データの読み込み (47 x 2 行列)
> x <- dat[,1]/100; y <- dat[,2]
> p <- 3; alpha <- 0.05
> func0079b(x,y,p,alpha)
# 次数= 0
RSS= 0.815383 , RSQ= 0
27
0.10
28
0.20
1.8
1.8
4.4
シミュレーションで確認
# run0080.R
# run0074.R のシミュレーション結果を利用する．
1.6
1.6
# 信頼区間，同時信頼区間：シミュレーション
y
source("run0079.R")
1.4
1.4
y
# run0075.R,run0079.R の関数をつかって信頼区間を計算する．
func0080a <- function(v) { # check if v1 in [v2,v3]
1.2
1.2
(v[1] >= v[2]) && (v[1] <= v[3])
}
func0080b <- function(be,y00,d) {
0.05
0.10
0.15
0.20
0.05
0.10
0.15
x
coef <- cbind(be,d$coef)
0.20
be0 <- apply(coef[,c("be","Lo","Up")],1,func0080a)
x
run0079-s2
be1 <- apply(coef[,c("be","JointLo","JointUp")],1,func0080a)
run0079-s3
pred0 <- all(apply(cbind(y00,d$b0$pred0conf),1,func0080a))
• すべて α = 0.05，つまり信頼係数 95%である．
pred1 <- all(apply(cbind(y00,d$b1$pred0conf),1,func0080a))
list(be0=be0,be1=be1,pred0=pred0,pred1=pred1)
• 回帰係数の信頼区間 (Up,Lo) と同時信頼区間 (JointLo,JointUp)
}
• 予測値の信頼区間（赤）と同時信頼区間（青）
a <- func0075a(x,p) # 準備
y00 <- a$X0 %*% be # 300 分割した点に対応する真の y
cat("# 一回目のシミュレーション結果のチェック\n")
同時信頼区間の導出＊
4.3
d <- func0079(x,simy[,1],a,alpha) # 信頼区間等の計算
• v = w0 β, w ∈ W の同時信頼区間を求める
abline(be,col=3,lty=4) # 真の回帰直線
0
• γ = Rβ, a = R −1 w と変換すると，v = w0 β = a0 γ とかける．このほうが，扱いやすい．
cat("# 信頼区間に入っていたか？\n")
• シュバルツ (Schwartz) の不等式より
dev.copy2eps(file="run0080-s1.eps")
|v̂ − v| = |a0 (γ̂ − γ) ≤ kak · kγ̂ − γk
cat("# シミュレーション結果のチェック\n")
yesno <- matrix(0,length(yesno1),b) # 結果をしまうアレイ
で等号は a と γ̂ − γ が平行のときのみ．
• γ の同時信頼域（球体）は
n
γ :
rownames(yesno) <- names(yesno1)
kγ − γ̂k ≤ S²
p
(p + 1)qf(1 − α)
であったから，γ がこの同時信頼域に入っていれば，
p
|v̂ − v| ≤ kakS² (p + 1)qf(1 − α)
for(i in 1:b) {
o
d <- func0079(x,simy[,i],a,alpha,output=F)
yesno[,i] <- unlist(func0080b(be,y00,d))
}
print(yesno[,1:5])
である．γ が同時信頼域に入る確率は 1 − α なので，
n
o
p
Pr |v̂ − v| ≤ kakS² (p + 1)qf(1 − α) ≥ 1 − α
である．
√
p
0
w0 R−1 R −1 w = w0 (X 0 X)−1 w に注意すれば，
½
¾
q
p
同時信頼区間w = v : |v̂ − v| ≤ S² w0 (X 0 X)−1 w (p + 1)qf(1 − α)
• ここで kak =
yesno1 <- unlist(func0080b(be,y00,d)); print(yesno1)
cat("# シミュレーションで信頼区間に入っていた回数\n")
print(apply(yesno,1,sum))
> n <- 11 # データサイズ
> be <- c(0,1) # 真の回帰係数 (beta0,beta1)
> x <- seq(0,1,len=n) # x を決める
> filename <- NULL
> source("run0074.R")
29
# start simulation:
# end simulation:
Wed Oct
Wed Oct
30
6 11:46:43 2004
> sum(yesno[1,] & yesno[2,]) # joint0
6 11:46:45 2004
[1] 9235
# １回目のシミュレーション結果
> sum(yesno[3,] & yesno[4,]) # joint1
$be
[1] 9750
[,1]
beta0 -0.1491011
1.3627077
0.8
beta1
$se2
0.6
[1] 0.04218202
0.4
$tval
y
[,1]
6.958817
0.0
beta1
0.2
beta0 -1.287003
−0.2
$pval
[,1]
beta0 2.302066e-01
0.0
beta1 6.619553e-05
0.2
0.4
0.6
0.8
1.0
x
run0080-s1
# pval0 < 0.05 の回数 =
508
# pval1 < 0.05 の回数 =
8751
• シミュレーションなので回帰係数の真値 β0 = 0, β1 = 1 が分かっている．これが信頼区
間に何回入ったかを数える．シミュレーションは 10000 回の繰り返しなので，信頼係数
=95%ならば，9500 回程度が理想的．
> alpha <- 0.05
> source("run0080.R")
• 個別の（同時でない）信頼区間
# 一回目のシミュレーション結果のチェック
#{β0 ∈ 信頼区間0 } = 9498
RSS= 0.5859556 , RSQ= 0.6594943
Estimate StdErr t-value p-value
beta0
-0.208
0.144
-1.447
beta1
1.016
0.243
4.175
Lo
0.182 -0.534 0.117
0.002
#{β1 ∈ 信頼区間1 } = 9453
Up JointLo JointUp
0.465 1.566
-0.628
0.212
0.306
1.726
であるから，個別に β0 と β1 を見るかぎりほぼ 95%．しかし，
#{β0 ∈ 信頼区間0 ,
# 信頼区間に入っていたか？
be0.beta0 be0.beta1 be1.beta0 be1.beta1
TRUE
TRUE
TRUE
TRUE
pred0
pred1
FALSE
TRUE
# シミュレーション結果のチェック
• 同時信頼区間
#{β0 ∈ 同時信頼区間0 ,
[,1] [,2] [,3] [,4] [,5]
be0.beta0
1
1
1
1
1
be0.beta1
1
1
1
0
1
be1.beta0
1
1
1
1
1
be1.beta1
1
1
1
1
1
pred0
0
1
1
0
1
pred1
1
1
1
1
1
• 予測値の信頼区間各 x で個別に見れば，
#{β0 + β1 x ∈ 信頼区間 (x), } = 約 9500,
9835
すべての x
のはずである．ところが
be0.beta0 be0.beta1 be1.beta0 be1.beta1
9453
β1 ∈ 同時信頼区間1 } = 9750
となり，同時に入るのは約 98%になる．
#{β0 + β1 x ∈ 信頼区間 (x),
# シミュレーションで信頼区間に入っていた回数
9498
β1 ∈ 信頼区間1 } = 9235
となり，同時に入るのは約 92%に低下する．
9825
31
pred0
pred1
8736
9519
すべての x} = 8736
なので，予測式としてみると信頼区間に入る確率は約 87%となり，95%よりかなり低下
する．
32
• 予測値の同時信頼区間
#{β0 + β1 x ∈ 同時信頼区間 (x),
すべての x} = 9519
なので，予測式としてみると信頼区間に入る確率は約 95%となり理想的．
課題
5
5.1
課題 7-1
X2000 データセットから自由に変量を選び多項式回帰を次数 p = 1, 2, 3 について行え．Ｒに
組み込み関数や講義で説明した関数など自由に用いてよい．自分で書いたプログラムはレポー
トに添付すること．
• それぞれの次数について，決定係数，および各回帰係数の推定値，標準誤差，t 統計量，
p-値を示せ．
• 結局，どの次数 p を用いるのが最も適切かを判断せよ．
• 上で選んだ次数について x,y の散布図上に推定した回帰直線（曲線）とその 95%信頼区間
（赤線），および 95%同時信頼区間（青線）を重ねて描け．
5.2
課題 7-2
run0075.R では次の二つの関数を定義している．
calcq0 <- function(n,p,alpha) qt(1-alpha/2,n-p-1) # 個別の信頼区間
calcq1 <- function(n,p,alpha) sqrt((p+1)*qf(1-alpha,p+1,n-p-1)) # 同時信頼区間
n = 30, α = 0.05 とおき，calcq0(n,p,alpha) と calcq1(n,p,alpha) を p = 0, 1, . . . , 10 の範
囲で計算して重ねてプロットせよ（横軸=p，縦軸=関数値とする)．これらの関数値を比較せ
よ．この違いは回帰直線 (曲面) の信頼区間，同時信頼区間ついてどのような結果をもたらすか
を説明せよ．
33