複数の回帰直線の比較 (01 February 2008)

某友人から電話があり、複数の回帰直線の比較の方法を尋ねられた。
その返事のメールを少々加工したものを以下に


複数の回帰直線の比較は以下のようにするといいようです。
多変量解析の一種。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")

3つのデータをひとつにまとめる

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)

今回の例の結果(実行は2008年1月31日)

ところどころに説明を挿入

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 で、有意でない。
よって、切片が違うとはいえない。

以上により、3つのカテゴリをまとめて問題はない。

なお、
切片 -5.3089
傾き 11.2396
(上のほうででてきています)

参考にしたウェブページ

Today:1Yesterday:1Total:5855 since 01 February 2008

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 01 Feb 2008 (金) 11:48:47 (4374d)