線形回帰分析 - 久留米大学

線形回帰分析
久留米大学バイオ統計センター
第 1 回公開セミナー
川口 淳
kawaguchi [email protected]
2005 年 7 月 2 日
線形回帰分析 – p.1/81
アウトライン
2.1 線形回帰モデル
2.2 傾きと切片に関する推測
2.2.3 回帰直線全体に関する推測
2.3 出力値からの入力値の推定
応用例 薬物用量と AUC
4.4 量的実験要因
線形回帰分析 – p.2/81
線形回帰モデル
反応変数 X を入力変数 Z に関連付けるために仮定
された統計モデル
X = α∗ + βZ + e
ただし，
α∗ : Z = 0 での X の期待値 (直線の切片)
β : Z の単位変化あたりの X の平均変化 (直線の傾き)
e : 平均 0，分散 σ 2 の正規分布に従う確率変数
線形回帰分析 – p.3/81
表 2.1
Z
X
1
0.045
2
0.114
4
0.215
6
0.346
7
0.41
8
0.52
10
0.67
15
0.942
Z : 溶液に溶解したアンチモン
のマイクログラム数
X : Beckmann の分光光度計
の示度
線形回帰分析 – p.4/81
R コード (データプロット)
データの読み込み．
z<-c(1, 2, 4, 6, 7, 8, 10, 15)
x<-c(0.045, 0.114, 0.215, 0.346,
0.410, 0.520, 0.670, 0.942)
横軸 z, 縦軸 x の二次元プロット．
plot(z,x)
線形回帰分析 – p.5/81
0.2
0.4
x
0.6
0.8
データプロット
2
4
6
8
10
12
14
z
線形回帰分析 – p.6/81
最小二乗推定量
n 組の観測値 (Z1 , X1 ), . . . , (Zn , Xn )
n
X
S2 =
(Xi − α∗ − βZi )2
i=1
を最小にするような a∗ ，b は次のように与えら
れる．
P
(Zi − Z̄)(Xi − X̄)
P
b =
(Zi − Z̄)2
a∗ = X̄ − bZ̄
線形回帰分析 – p.7/81
R コード (最小二乗推定量)
組み込み関数を使用
lm(x˜z)
アウトプット
Call:
lm(formula = x ˜ z)
Coefficients:
(Intercept)
z
-0.02777
0.06574
線形回帰分析 – p.8/81
SAS コード (最小二乗推定量)
data beckmann;
input z x @@;
datalines;
1 0.045 2 0.114 4 0.215 6 0.346
7 0.410 8 0.520 10 0.670 15 0.942
run;
proc reg data=beckmann;
model x=z ;
run;
線形回帰分析 – p.9/81
SAS 結果 (最小二乗推定量)
SAS アウトプットの一番下の表
線形回帰分析 – p.10/81
同値な表現
統計モデルを
X = α + β(Z − Z̄) + e
と表したとき，α と β の最小二乗推定量は
P
(Zi − Z̄)(Xi − X̄)
P
a = X̄,
b=
(Zi − Z̄)2
回帰直線は Z と X の平均を通る．
a と b は無相関．
線形回帰分析 – p.11/81
共分散と相関係数
Z と X の共分散の推定量 sZX
P
(Zi − Z̄)(Xi − X̄)
sZX =
n−1
⇒
sZX
b=
sX
Z と X の相関係数の推定量 r (−1 ≤ r ≤ 1)
sZX
r=
sZ sX
⇒
sX
b=r
sZ
ただし，
P
1
sZ = n−1 (Zi − Z̄)2 ，sX =
1
n−1
P
(Xi − X̄)2
線形回帰分析 – p.12/81
決定係数
2
R =
S12 −S22
S12
のことを決定係数と呼ぶ．ただし，
X
(Xi − X̄)2 ,
S12 =
S22 =
X
[Xi − X̄ −b(Zi − Z̄)]2
決定係数は，データの変動のうち求めた回帰直線
が説明できる割合．
他方，
R2 = r 2
が成立し，相関係数の二乗と決定係数は一致する．
線形回帰分析 – p.13/81
測定の信頼性
b は βU,T RZ の推定となる．
ただし，
X = U +f,
Z = T +e,
σT2
RZ = 2
,
2
σT + σe
e と f は無相関確率変数
σU2
RX = 2
σU + σf2
βU,T は U と V の傾き
傾きの推定量は Z の測定信頼度だけに依存する．
Z の信頼度が下がると，傾きは小さくなる．
線形回帰分析 – p.14/81
2
e の分散 σ の不偏推定量
誤差分散の推定量は以下のように与えられ，回帰
パラメータの推測に用いられる．
P
∗
2
(X
−
a
−
bZ
)
i
i
2
s =
n−2
P
(Xi − X̄ − a∗ − b(Zi − Z̄))2
=
n−2
不偏推定量 (E[s2 ] = σ 2 ) にするために n − 2 で割る．
線形回帰分析 – p.15/81
傾きに関する推測 (1)
β の 100(1 − α)% 信頼区間
b − tn−2,α/2 se(b) ≤ β ≤ b + tn−2,α/2 se(b)
ただし，
se(b) = pP
s
(Zi − Z̄)2
s ：誤差分散の推定量の平方根
tν,p ：自由度 ν の t 分布上側 p% 点．
線形回帰分析 – p.16/81
傾きに関する推測 (2)
H0 : β = 0 の検定 (有意水準 α%)
b
t=
se(b)
t > tn−2,α/2
(H0
(検定統計量)
=⇒
H0
reject
reject ⇐⇒ 信頼区間が 0 を含まない)
線形回帰分析 – p.17/81
傾きの検定と相関係数の検定
検定統計量 t に対し，
√
b
r n−2
t=
= √
se(b)
1 − r2
が成り立つ．
故に
H0 : β = 0 ⇐⇒ H0 : ρ = 0
ただし，ρ は真の相関係数を表す．
線形回帰分析 – p.18/81
切片に関する推測 (1)
α∗ の 100(1 − α)% 信頼区間
a∗ − tn−2,α/2 se(a∗ ) ≤ α∗ ≤ a∗ + tn−2,α/2 se(a∗ )
ただし，
s
se(a∗ ) = s
1
Z̄ 2
+P
n
(Zi − Z̄)2
線形回帰分析 – p.19/81
切片に関する推測 (2)
H0 : α∗ = 0 の検定 (有意水準 α%)
a∗
t=
se(a∗ )
t > tn−2,α/2
=⇒
H0
reject
線形回帰分析 – p.20/81
R コード (傾き・切片に関する検定)
組み込み関数の使用
coefficients(summary(lm(x˜z)))
アウトプット
Estimate
Std.Error
t value
Pr(>|t|)
(Intercept)
-0.028
0.017
-1.668
0.146
z
0.066
0.002
31.063 7.39e-08
切片は 0 と有意に異なると言えないが，
傾きは 0 と有意に異なる．
線形回帰分析 – p.21/81
R コード (傾き・切片の信頼区間)
組み込み関数の使用
confint(lm(x˜z))
アウトプット
2.5 %
97.5 %
(Intercept) -0.069
0.013
z
0.061
0.071
切片の信頼区間は 0 を含むが，
傾きのは 0 を含まない．
線形回帰分析 – p.22/81
SAS コード (傾き・切片の信頼区間)
オプション (CLB) をつける．
proc reg data=beckmann;
model x=z / CLB;
run;
結果
線形回帰分析 – p.23/81
原点を通ると仮定された直線 (1)
直線が原点を通ると仮定があるときの統計モデル
X = βZ + e
β の最小二乗推定量
P
Zi Xi
β̂ = P 2
Zi
線形回帰分析 – p.24/81
原点を通ると仮定された直線 (2)
β の 100(1 − α)% 信頼区間
β̂ − tn−2,α/2 se(β̂) ≤ β ≤ β̂ + tn−2,α/2 se(β̂)
ただし，
v
se(β) = pP
Zi2
1 X
v =
(Xi − β̂Zi )2
n−1
線形回帰分析 – p.25/81
R コード (原点を通る直線)
summary(lm(x˜z-1))
confint(lm(x˜z-1))
アウトプット（一部分）
Estimate
Std.Error
t value
Pr(>|t|)
0.063
0.001
49.11
<0.001
z
信頼区間
0.060 ≤ β ≤ 0.066
線形回帰分析 – p.26/81
SAS コード (原点を通る直線)
オプション (NOINT) を加える．
proc reg data=beckmann;
model x=z /NOINT CLB;
run;
アウトプット（一部分）
線形回帰分析 – p.27/81
回帰直線に関する推測
ある Z = Z0 に対する
α∗ + βZ0 の 100(1 − α)% 信頼区間
s
2
1
(Z
−
Z̄)
0
a∗ + bZ0 ± tn−2,α/2 × s ×
+P
n
(Zi − Z̄)2
Z0 = 0 のときは，切片の信頼区間と一致する．
線形回帰分析 – p.28/81
回帰直線全体に関する推測 (考え方 1)
Z = Zi (i = 1, 2, . . . ) で同時に信頼区間を作りたい．
=⇒ 多重性の調整 (Bonferroni など)
または，

∗
∗
2
{a
+
bZ
−
(α
+
βZ)}
≤ Mα ,
P
2
(Z−Z̄)
1
2
P
s n + (Z −Z̄)2
i
なる Mα を求める．

∀Z  = 1−α
線形回帰分析 – p.29/81
回帰直線全体に関する推測 (考え方 2)
同値なことで，


∗
∗
2
{a
+
bZ
−
(α
+
βZ)}


P max
≤ Mα  = 1 − α
2
Z̄)
Z
s2 n1 + P(Z−
(Zi −Z̄)2
なる Mα を求める．
=
(X̄−E[X̄])2
σ 2 /n
2
(b−β)
+ σ2 / P(Zi −Z̄)2
s2 /σ 2
∼ 2F2,n−2
となる．
線形回帰分析 – p.30/81
回帰直線全体に関する推測
α∗ + βZ の 100(1 − α)% 信頼区間
L(Z ∗ ) < α∗ + βZ < U (Z ∗ )
ただし，Z ∗ はある範囲の任意の Z の値．
L(Z ∗ ) = X̄ + b(Z ∗ − Z̄) − A(Z ∗ )
U (Z ∗ ) = X̄ + b(Z ∗ − Z̄) + A(Z ∗ )
s
∗ − Z̄)2
1
(Z
2F2,n−2,α s2
A(Z ∗ ) =
+P
n
(Zi − Z̄)2
線形回帰分析 – p.31/81
R コード (問題 2.8) その 1
n<-length(z)
new.z<- seq(2, 12, by=0.02)
ones<-rep(1, length(new.z))
pre.z<-rbind(ones, new.z)
a<-lm(x˜z)$coefficients
plot(new.z, a%*%pre.z, ylim=c(0,1),
type="l", lwd=2, xlab="z", ylab="x")
線形回帰分析 – p.32/81
R コード (問題 2.8) その 2
sz<-sum((z-mean(z))ˆ2)
s2<-sum((lm(x˜z)$residuals)ˆ2)/(n-2)
s1<-1/n+(new.z-mean(z))ˆ2/sz
U<-a%*%pre.z+sqrt(2*qf(0.95,2,6)*s2*s1)
L<-a%*%pre.z-sqrt(2*qf(0.95,2,6)*s2*s1)
points(new.z, U, type="l")
points(new.z, L, type="l")
points(z, x)
線形回帰分析 – p.33/81
0.0
0.2
0.4
x
0.6
0.8
1.0
問題 2.8 結果
2
4
6
8
10
12
z
線形回帰分析 – p.34/81
R コード (点ごとの信頼区間)
p.U<-a%*%pre.z+sqrt(qt(0.95,2,6)*s2*s1)
p.L<-a%*%pre.z-sqrt(qt(0.95,2,6)*s2*s1)
points(new.z, p.U, type="l",lty=3)
points(new.z, p.L, type="l",lty=3)
線形回帰分析 – p.35/81
0.0
0.2
0.4
x
0.6
0.8
1.0
点ごとの信頼区間を重ねて図示
2
4
6
8
10
12
z
線形回帰分析 – p.36/81
出力値からの入力値の推定 (逆補間)
直線のあてはめ際，あらかじめ決められる変数は
入力値にする．(表 2.1 データではアンチモンの量)
観測された X(出力値) から未知の Z(入力値) を推
定する．
X − X̄
Ẑ = Z̄ +
b
Z の信頼区間 =⇒ Fieller の定理
線形回帰分析 – p.37/81
Fieller の定理
U ∼ N (µ, sa1 ), V ∼ N (ν, sa2 ),
s : 推定分散の平方根
g=
t2m,α/2 s2 a22
V2
と定義，
g < 1 と仮定 (V が 0 と有意に異なることと同値)．
µ/ν の 100(1 − α)% 信頼区間
!
r
2 a2
tm,α/2 s
1
U
U
2
(1 − g)a21 +
±
1−g V
V
V2
線形回帰分析 – p.38/81
Fieller の定理の解説
θ = µ/ν とするとき，
E[U − θV ] = 0,
V̂ [U − θV ] = s2 (a21 + θ2 a22 )
θ の 100(1 − α)% 信頼区間
2
(U − θV )
2
θ: 2 2
≤
t
m,α/2
2
2
s (a1 + θ a2 )
この二次不等式は，tm,α/2 s2 a22 /V 2 < 1 であれば，
両側の限界値が求まる．
線形回帰分析 – p.39/81
Z の信頼区間 (Fieller の定理の適用)
U = X − X̄ ，V = b とする．
q
se(U ) = s n+1
se(V ) = √P
n
⇓
a21 = (n + 1)/n
=⇒
g=
⇓
P
2
a2 = 1/ (Zi − Z̄)2
t2m,α/2 s2 a22
V2
s
(Zi −Z̄)2
t2n−2,α/2 s2
= 2P
b
(Zi − Z̄)2
線形回帰分析 – p.40/81
Z の信頼区間 (Fieller の定理の適用)
Fieller の定理に適用すると，
E[U ]
β(Z−Z̄)
= (Z − Z̄) の 100(1 − α)% 信頼区間
E[V ] =
β
X − X̄ tn−2,α/2 s
±
b(1 − g) b(1 − g)
を得る．
s
(1 − g)(n + 1)
(X − X̄)2
+ 2P
n
b
(Zi − Z̄)2
Z̄ を移項すると，Z の 100(1 − α)% 信頼区間になる．
線形回帰分析 – p.41/81
R コード (Z の信頼区間)
new.x<-0.577
hat.z<-mean(z)+(new.x-mean(x))/a[2]
g1<-s2*(qt(0.975, (length(z)-2)))ˆ2
g2<-a[2]ˆ2*sum((z-mean(z))ˆ2)
g<-g1/g2
m1<-mean(z)+(new.x-mean(x))/(a[2]*(1-g))
m2<-qt(0.975, (n-2))*sqrt(s2)/(a[2]*(1-g))
m31<-((1-g)*(n+1))/(n)
m32<-(new.x-mean(x))ˆ2/(a[2]ˆ2*sz)
m1+m2*sqrt(m31+m32)
m1-m2*sqrt(m31+m32)
線形回帰分析 – p.42/81
Z の信頼区間 (結果)
Z の予測値
> hat.z
9.199561
Z の 95% 信頼区間
> m1+m2*sqrt(m31+m32)
10.24142
> m1-m2*sqrt(m31+m32)
8.189849
8.19 ≤ Z ≤ 10.24
線形回帰分析 – p.43/81
応用例：薬物用量と AUC
0
5
10
15
20
time
Subject 1
25
12
0
2
4
6
concentrations
8
10
12
10
8
0
2
4
6
concentrations
6
4
2
0
concentrations
8
10
12
AUC(Area Under the bood - level Curve)
・
・
・血中濃度曲線下面積
0
5
10
15
20
time
Subject 2
25
0
5
10
15
20
25
time
Subject 3
線形回帰分析 – p.44/81
データ
SAS の NLMIXED プロシジャの Example
(Pinheiro and Bates, 1995)
time(測定時間) conc(血中濃度) dose(薬物用量)
12 人が 25 時間で 11 回計測．
AUC は，R の PK パッケージの AUC 関数を用
いて計算．
線形回帰分析 – p.45/81
R コード (AUC の計算)
data<-read.table("AUC_data.txt", header=TRUE)
library(PK)
auc<-matrix(0, 12, 3)
for(i in 1:12)
{auc[i,]<AUC(conc=data$CONC[data$SUBJECT==i],
time=data$TIME[data$SUBJECT==i])$AUC}
colnames(auc) <- c("obs","inter","inf")
auc<-data.frame(auc)
線形回帰分析 – p.46/81
モデル
線形単回帰モデル
log(AU C) = α + β log(用量) + ε
冪モデル
AU C = eα (用量)β eε
(
β=
6 1 =⇒ AU C と用量は非線形関係
β = 1 =⇒ AU C と用量は線形関係
線形回帰分析 – p.47/81
H0 : β = 1 の検定
反応変数を log(AU C) − log(用量) にかえて，
モデル
∗
log(AU C) − log(用量) = α + β log(用量) + ε
より，β ∗ を推定して，
H0 : β ∗ = 0
の検定を行う．
線形回帰分析 – p.48/81
R コード (直線あてはめ)
dose<-numeric(0)
for(i in 1:12)
{dose[i]<-data$DOSE[data$SUBJECT==i]}
ones<-rep(1, length(dose))
predictor<-rbind(ones, log(dose))
plot(log(dose), log(auc$inf),
xlab="log(DOSE)",ylab="log(AUC)",ylim=c(4.4, 5.4))
points(log(dose),
lm(log(auc$inf)˜log(dose))$coe%*%predictor,type="l")
線形回帰分析 – p.49/81
4.4
4.6
4.8
log(AUC)
5.0
5.2
5.4
直線あてはめ結果
1.2
1.3
1.4
1.5
1.6
1.7
log(DOSE)
log(AU C) = 4.13 + 0.42 log(用量)
線形回帰分析 – p.50/81
R コード (H0 : β = 1 の検定)
coefficients(summary(
lm((log(auc$inf)-log(dose))˜log(dose))))
Estimate Std. Error t value Pr(>|t|)
(Intercept)
4.13
0.75
5.47 <0.001
log(dose)
−0.58
0.49 −1.17
0.27
従って，傾きは 1 と有意に異なると言えない．
線形回帰分析 – p.51/81
R コード (β の信頼区間)
confint(lm((log(auc$inf))˜log(dose)))
2.5 % 97.5 %
(Intercept)
2.44
5.81
log(dose) −0.68
1.52
線形回帰分析 – p.52/81
R コード (回帰診断)
plot(lm(log(auc$inf)˜log(dose)))
残差プロット・
・
・データに異常値が含まれてい
るか，誤差の等分散性が満足しているかを調
べる．
Q-Q プロット・
・
・誤差について正規性を満足し
ているか調べる．
Cook の距離・
・
・そのデータの回帰パラメータ
への影響を調べる．
線形回帰分析 – p.53/81
回帰診断
Scale-Location plot
Cook’s distance plot
0.5
0.4
1
0.2
0.5
0.3
10
Cook’s distance
1.0
11
0.1
10
6
0.0
0.0
Standardized residuals
1.5
1
4.60
4.70
4.80
Fitted values
2
4
6
8
10
12
Obs. number
線形回帰分析 – p.54/81
4.4
4.6
4.8
log(AUC)
5.0
5.2
5.4
直線あてはめ結果 (取り除き後)
1.2
1.3
1.4
1.5
1.6
1.7
log(DOSE)
log(AU C) = 3.59 + 0.73 log(用量)
線形回帰分析 – p.55/81
H0 : β = 1 の検定 (取り除き後)
coefficients(summary(
lm((log(auc$inf)[-1]-log(dose)[-1])
˜log(dose)[-1])))
Estimate Std. Error t value Pr(>|t|)
(Intercept)
3.59
0.48
7.44 <0.001
log(dose)[-1]
−0.27
0.31 −0.86
0.41
傾きは 1 と有意に異なると言えない．
線形回帰分析 – p.56/81
β の信頼区間 (取り除き後)
confint(
lm((log(auc$inf)[-1])˜log(dose)[-1]))
2.5 % 97.5 %
(Intercept) 2.50
4.69
log(dose)[-1] 0.02
1.44
線形回帰分析 – p.57/81
結果のまとめ
log(AU C) = 4.13 + 0.42 log(用量)
H0 : β = 1 は棄却されない．
信頼区間 −0.68 ≤ β ≤ 1.52
AUC と用量は非線形関係とは言えない．
(Subject 1 を取り除く)
log(AU C) = 3.59 + 0.73 log(用量)
H0 : β = 1 は棄却されない．
信頼区間 0.02 ≤ β ≤ 1.44
AUC と用量は非線形関係とは言えない．
線形回帰分析 – p.58/81
回帰診断 (Subject1 を取り除いた後)
1.4
Scale-Location plot
Cook’s distance plot
8
2.5
9
8
2.0
1.5
1.0
Cook’s distance
1.0
0.8
0.6
0.5
0.4
0.2
9
0.0
10
0.0
Standardized residuals
1.2
10
4.5
4.6
4.7
Fitted values
4.8
4.9
2
4
6
8
10
Obs. number
線形回帰分析 – p.59/81
4.4 節 量的実験要因
表 4.6 データ
X:ラットの病の回復度合い
Z:ビタミン D の投与用量
群
用量 (Zi )
1
3.5(0)
2
7.0(1)
8 3.31
1.63
3
14.0(2)
13 3.42
1.02
ni
X̄i
10 1.8
si
1.09
Zi = log2 (用量/3.5)
si : 各群における X の標準偏差
線形回帰分析 – p.60/81
解析項目
モデルの次数に関する検定
Z 2 の係数が有意に 0 かどうかの検定
Z の係数が有意に 0 かどうかの検定
結論
線形回帰分析 – p.61/81
モデルの次数に関する推測
X 0 i : Zi での X の予測値．(i = 1, 2, . . . , g)
H0 : 真の式は次数が 2 次以下
P
ni (X̄i − X 0 i )2
F =
(g − 3)s2
F > Fg−3,ν,α/3
=⇒
H0 reject
表 4.6 データでは群が 3 つしかないので，
この検定は行わない．
線形回帰分析 – p.62/81
二次式のあてはめ
X と Z が曲線関係と仮定したときの統計モデル
X̄ = α + βZ + γZ 2 + ε
誤差分散 (プールした群内の分散) の推定量 s
P
2
(n
−
1)s
i
i
s2 = P
(ni − 1)
線形回帰分析 – p.63/81
最小二乗推定量
最小二乗推定量 a，b，c
 P
 

ni X̄i
a
 

−1  P
ni X̄i Zi 
 b =M 
P
c
ni X̄i Zi2
ただし，

 P
P
P
ni Zi2
ni Zi
ni
P
P

 P
2
3
M =
ni Zi
ni Zi 
ni Zi
P
P
P
2
3
ni Zi
ni Zi
ni Zi4
線形回帰分析 – p.64/81
R コード (最小二乗推定量)
z<-c(0, 1, 2)
n<-c(10, 8, 13)
x<-c(1.80, 3.31, 3.42)
M<-c(sum(n), sum(n*z), sum(n*zˆ2),
sum(n*z), sum(n*zˆ2), sum(n*zˆ3),
sum(n*zˆ2), sum(n*zˆ3), sum(n*zˆ4))
M<-matrix(M, ncol=3)
M2<-c(sum(n*x), sum(n*x*z), sum(n*x*zˆ2))
a<-solve(M)%*%M2
線形回帰分析 – p.65/81
R コード (最小二乗推定量:結果)
> a
[,1]
[1,] 1.80
[2,] 2.21
[3,] -0.70
a = 1.80，b = 2.21，c = −0.70
線形回帰分析 – p.66/81
2 次の項の係数に関する推測
H0 : Z 2 の係数 γ = 0 有意水準 (α/2%)
c
Lc =
se(c)
ただし，√
se(c) = s m(3,3) ，
m(3,3) は M −1 の 3 行 3 列成分
|Lc | > tn−3,α/4
=⇒
H0 reject
線形回帰分析 – p.67/81
R コード (2 次の項の係数に関する推測)
s2<-((n-1)%*%sˆ2)/sum(n-1)
se.c<-sqrt(s2*solve(M)[3,3])
L.c<-a[3]/se.c
p.c<-1-pt(abs(L.c), 28, 0.0125)#ｐ値
アウトプット
> p.c
0.089
H0 : γ = 0 は有意水準 0.025% で棄却できない．
線形回帰分析 – p.68/81
回帰曲線全体に関する推測 (1)
X 0 Z : ある値 Z における反応の期待値
 
a

0
2
0
X Z = a + bZ + cZ = Z  b 
c
ただし，Z0 = (1, Z, Z 2 )
線形回帰分析 – p.69/81
回帰曲線全体に関する推測 (2)
2 次式の場合の
Working-Hotelling 信頼帯 (信頼度 100(1 − α)%)
p
X 0 Z ± 3F3,n.−g,α se(X 0 z )
ただし，
n. =
X
ni
√
se(X 0 Z ) = s Z0 M−1 Z
線形回帰分析 – p.70/81
R コード (回帰曲線全体に関する推測)
U<-new.Z%*%a+sqrt(3*qf(0.95, 3, sum(n)-3))*se.X
L<-new.Z%*%a-sqrt(3*qf(0.95, 3, sum(n)-3))*se.X
points(new.z, U, type="l")
points(new.z, L, type="l")
points(z, x)
new.z<-seq(0, 3, 0.2)
ones<-rep(1, length(new.z))
new.Z<-cbind(ones, new.z, new.zˆ2)
se.X<-sqrt(s2*diag(new.Z%*%solve(M)%*%t(new.Z)))
plot(new.z, new.Z%*%a, ylim=c(0,6),
type="l", xlab="z", ylab="x")
線形回帰分析 – p.71/81
3
2
1
0
x
4
5
6
回帰曲線全体に関する推測：結果
0.0
0.5
1.0
1.5
2.0
2.5
3.0
z
線形回帰分析 – p.72/81
一次式のあてはめ
モデル
X̄ = α∗ + βZ + ε
最小二乗推定量
P
ni (Zi − Z̄)(Xi − X̄)
P
β̂ =
ni (Zi − Z̄)2
P
P
ni X̄i − β̂ ni Zi
∗
P
α̂ =
ni
線形回帰分析 – p.73/81
一次式の最小二乗推定値
alpha<-((sum(n*x)-beta*sum(n*z))/(sum(n)))
beta<-((sum(n)*sum(n*x*z)-(sum(n*x))*(sum(n*z)))
/(sum(n)*sum(n*zˆ2)-sum(n*z)ˆ2))
> beta
[1] 0.786
> alpha
[1] 2.007
β̂ = 0.786, α̂∗ = 2.007
線形回帰分析 – p.74/81
一次式の傾きに関する推測
H0 : β = 0 の検定 (有意水準 α/2%)
Lβ̂ =
β̂
se(β̂)
ただし，
pP
s
ni
se(β̂) = pP
ni (Zi − Z̄)2
|Lβ̂ | > tn−3,α/4
=⇒
H0 reject
線形回帰分析 – p.75/81
R コード (一次式の傾きに関する推測)
se.beta<-(sqrt((s2*sum(n))/
(sum(n)*sum(n*zˆ2)-sum(n*z)ˆ2)))
L.beta<-beta/se.beta
p.beta<-1-pt(L.beta, 28, 0.0125)
線形回帰分析 – p.76/81
一次式の傾きに関する推測:結果
アウトプット
> L.beta
3.067
> p.beta
0.002
H0 : β = 0 は有意水準 0.025% で棄却できる．
線形回帰分析 – p.77/81
R コード (一次式と二次式の図示)
a1<-coefficients(lm(x˜z, weight=n))
ones<-rep(1, length(new.z))
new.Z1<-cbind(ones, new.z)
points(new.z, new.Z1%*%a1, type="l")
線形回帰分析 – p.78/81
3
2
1
0
x
4
5
6
一次式と二次式の図示：結果
0.0
0.5
1.0
1.5
2.0
2.5
3.0
z
線形回帰分析 – p.79/81
表 4.6 データの結果のまとめ
二次式では Z 2 の係数は有意に 0 と異なるとは
言えなかった．
一次式では Z の係数は有意に 0 と異なると言
えた．
一次式の Z の係数は正となり，ビタミン D の用
量を増やすと回復具合が増す．
用量を増やし続けると飽和状態になる可能性も
含め，外挿データ部分については議論は出来
ない．
線形回帰分析 – p.80/81
参考文献
J.L. フライス (KR 研究会 訳) (2004)，臨床試験のデザイ
ンと解析，アーム.
RjpWiki，
http://www.okada.jp.org/RWiki/index.php?RjpWiki
Pinheiro, J.C. and Bates, D.M. (1995), "Approximations
to the Log-likelihood Function in the Nonlinear
Mixed-effects Model," Journal of Computational and
Graphical Statistics, 4, 12-35.
線形回帰分析 – p.81/81