![[PukiWiki]](image/tomarufront.png "[PukiWiki]")
某友人から電話があり、複数の回帰直線の比較の方法を尋ねられた。
その返事のメールを少々加工したものを以下に
複数の回帰直線の比較は以下のようにするといいようです。
多変量解析の一種。ANCOVAとしての扱いにもなります。
ウェブ上でデータを入力するようにしているので、データの入力がややこしい のが難点で、ごめんなさい。
まず、以下のようなものを用意してください
x<-c(1,2,3,4,5,6,7,8,9,10,2,3,4,5,5.5,6,7,8,9,10,2,2,3,4,5,5,7,8,10,10)
y<-c(10,22,30,39,44,60,70,89,90,110,22,20,30,40,50,66,66,89,99,109,12,20,30,40,55,55,70,80,99,105)
Cat<-c("A","A","A","A","A","A","A","A","A","A","B","B","B","B","B","B","B","B","B","B","C","C","C","C","C","C","C","C","C","C")
mydata<-data.frame(x,y,Cat)
mydata
fitall<-lm(y~x,mydata)
fitint<-lm(y~Cat+x,mydata)
fitdif<-lm(y~Cat*x,mydata)
summary(fitall)
summary(fitint)
summary(fitdif)
anova(fitint,fitdif)
anova(fitall,fitint)
http://bayes.math.montana.edu/Rweb/Rweb.general.htlm
に行き、入力ボックスに上のコードを貼り付け、submitボタンを押します。
結果がでてきますので、それを保存してください。
作業は以上です。
今回のデータは、3種類（A、B、C）のカテゴリ（Catとしています）からなります。 xが横軸、yが縦軸です。 必要に応じて直線回帰できるように変換しておいてください。
すべてのカテゴリのデータをx、y、Catの順で入力します。
データセットとしては
1, 2, 3, 4,（略）, 9, 10, 2, 3,（略）, 10, 2,（略）, 10 10, 22, 30, 39,（略）, 90,110, 22, 20,（略）,109, 12,（略）,105 "A","A","A","A",（略）,"A","A","B","B",（略）,"B","C",（略）,"C"
こんな感じのものを作り、一行ずつ入力していくイメージです。
縦がずれないように注意してください。
これで、xを入力
x<-c(1,2,3,4,5,6,7,8,9,10,2,3,4,5,5.5,6,7,8,9,10,2,2,3,4,5,5,7,8,10,10)
これで、yを入力
y<-c(10,22,30,39,44,60,70,89,90,110,22,20,30,40,50,66,66,89,99,109,12,20,30,40,55,55,70,80,99,105)
これで、Catを入力
Cat<-c("A","A","A","A","A","A","A","A","A","A","B","B","B","B","B","B","B","B","B","B","C","C","C","C","C","C","C","C","C","C")
３つのデータをひとつにまとめる
mydata<-data.frame(x,y,Cat)
確認のためデータセットを表示
mydata
すべてのデータをひとつにして回帰分析
fitall<-lm(y~x,mydata)
切片は異なるが、傾きは同じものとして回帰分析
fitint<-lm(y~Cat+x,mydata)
切片も傾きも異なるものとして回帰分析
fitdif<-lm(y~Cat*x,mydata)
すべてのデータをひとつにして行った回帰の結果
summary(fitall)
切片は異なるが、傾きは同じものとして行った回帰の結果
summary(fitint)
切片も傾きも異なるものとして行った回帰の結果
summary(fitdif)
「切片は異なるが、傾きは同じもの」と「切片も傾きも異なるもの」を比較
これで有意でなければ、傾きに有意差がないと言える
anova(fitint,fitdif)
「すべてのデータをひとつにしてもの」と「切片は異なるが、傾きは同じもの」
を比較
傾きが有意でなかった場合で、かつ、こちらも有意でなければ、切片にも有意
差がないと言える
よって、傾きも切片も有意差がないので、プールしてよいといえる。
anova(fitall,fitint)
ところどころに説明を挿入
Results from Rweb You are using Rweb1.03 on the server at bayes.math.montana.edu R : Copyright 2003, The R Foundation for Statistical Computing Version 1.8.1 (2003-11-21), ISBN 3-900051-00-3 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. Rweb:> postscript(file= "/tmp/Rout.27085.%03d.ps", onefile=FALSE) Rweb:>
ここまでは、念仏です
Rweb:> x<-c(1,2,3,4,5,6,7,8,9,10,2,3,4,5,5.5,6,7,8,9,10,2,2,3,4,5,5,7,8,10,10)
Rweb:> y<-c(10,22,30,39,44,60,70,89,90,110,22,20,30,40,50,66,66,89,99,109,12,20,30,40,55,55,70,80,99,105)
Rweb:> Cat<-c("A","A","A","A","A","A","A","A","A","A","B","B","B","B","B","B","B","B","B","B","C","C","C","C","C","C","C","C","C","C")
Rweb:> mydata<-data.frame(x,y,Cat)
下にデータセットが表示されています。
Rweb:> mydata
x y Cat
1 1.0 10 A
2 2.0 22 A
3 3.0 30 A
4 4.0 39 A
5 5.0 44 A
6 6.0 60 A
7 7.0 70 A
8 8.0 89 A
9 9.0 90 A
10 10.0 110 A
11 2.0 22 B
12 3.0 20 B
13 4.0 30 B
14 5.0 40 B
15 5.5 50 B
16 6.0 66 B
17 7.0 66 B
18 8.0 89 B
19 9.0 99 B
20 10.0 109 B
21 2.0 12 C
22 2.0 20 C
23 3.0 30 C
24 4.0 40 C
25 5.0 55 C
26 5.0 55 C
27 7.0 70 C
28 8.0 80 C
29 10.0 99 C
30 10.0 105 C
回帰を行う
Rweb:> fitall<-lm(y~x,mydata) Rweb:> fitint<-lm(y~Cat+x,mydata) Rweb:> fitdif<-lm(y~Cat*x,mydata)
すべてのカテゴリをまとめて回帰したもののサマリ ここに、切片（-5.3089）と傾き（11.2396）。また、それぞれ有意かどうか。 今回はどちらも有意なので、切片は0ではなく、傾きも0ではない。
Rweb:> summary(fitall)
Call:
lm(formula = y ~ x, data = mydata)
Residuals:
Min 1Q Median 3Q Max
-9.831 -4.546 0.697 5.050 6.085
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.3089 2.1550 -2.464 0.0202 *
x 11.0280 0.3414 32.301 <2e-16 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 5.135 on 28 degrees of freedom
Multiple R-Squared: 0.9739, Adjusted R-squared: 0.9729
F-statistic: 1043 on 1 and 28 DF, p-value: < 2.2e-16
切片は異なるが、傾きは同じものとして行った回帰のサマリ
Rweb:> summary(fitint)
Call:
lm(formula = y ~ Cat + x, data = mydata)
Residuals:
Min 1Q Median 3Q Max
-8.601 -4.560 0.428 4.799 7.245
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.3828 2.5340 -1.730 0.0956 .
CatB -2.2731 2.3461 -0.969 0.3415
CatC -0.9051 2.3411 -0.387 0.7022
x 11.0514 0.3489 31.679 <2e-16 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 5.234 on 26 degrees of freedom
Multiple R-Squared: 0.9748, Adjusted R-squared: 0.9719
F-statistic: 335.1 on 3 and 26 DF, p-value: < 2.2e-16
切片も傾きも異なるものとして行った回帰のサマリ
Rweb:> summary(fitdif)
Call:
lm(formula = y ~ Cat * x, data = mydata)
Residuals:
Min 1Q Median 3Q Max
-7.5289 -3.5063 0.3771 3.2058 11.0116
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.9333 3.4268 -0.856 0.4005
CatB -10.4386 5.3899 -1.937 0.0646 .
CatC 0.7877 4.8822 0.161 0.8732
x 10.7879 0.5523 19.533 3.06e-16 ***
CatB:x 1.3923 0.8502 1.638 0.1146
CatC:x -0.2976 0.7813 -0.381 0.7066
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 5.016 on 24 degrees of freedom
Multiple R-Squared: 0.9786, Adjusted R-squared: 0.9742
F-statistic: 219.8 on 5 and 24 DF, p-value: < 2.2e-16
「切片は異なるが、傾きは同じもの」と「切片も傾きも異なるもの」を比較
Rweb:> anova(fitint,fitdif) Analysis of Variance Table Model 1: y ~ Cat + x Model 2: y ~ Cat * x Res.Df RSS Df Sum of Sq F Pr(>F) 1 26 712.33 2 24 603.93 2 108.40 2.154 0.1379
確率（Pr(>F)）が 0.1379 で、有意でない。 よって、傾きが違うとはいえない。
Rweb:> anova(fitall,fitint) Analysis of Variance Table Model 1: y ~ x Model 2: y ~ Cat + x Res.Df RSS Df Sum of Sq F Pr(>F) 1 28 738.40 2 26 712.33 2 26.07 0.4757 0.6267
確率（Pr(>F)）が 0.6267 で、有意でない。
よって、切片が違うとはいえない。
以上により、３つのカテゴリをまとめて問題はない。
なお、
切片 -5.3089
傾き 11.2396
（上のほうででてきています）
| Today:1 | Yesterday:3 | Total:11470 since 01 February 2008 |
