某友人から電話があり、複数の回帰直線の比較の方法を尋ねられた。
その返事のメールを少々加工したものを以下に
複数の回帰直線の比較は以下のようにするといいようです。
多変量解析の一種。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)
ところどころに説明を挿入
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:2 | Yesterday:7 | Total:11487 since 01 February 2008 |