![[PukiWiki]](image/tomarufront.png "[PukiWiki]")
*複数の回帰直線の比較 （01 February 2008） [#s3652589] 某友人から電話があり、複数の回帰直線の比較の方法を尋ねられた。~ その返事のメールを少々加工したものを以下に ---- 複数の回帰直線の比較は以下のようにするといいようです。~ 多変量解析の一種。ANCOVAとしての扱いにもなります。 ウェブ上でデータを入力するようにしているので、データの入力がややこしい のが難点で、ごめんなさい。 **用意するもの [#a6247fc1] まず、以下のようなものを用意してください 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) 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ボタンを押します。~ 結果がでてきますので、それを保存してください。~ 作業は以上です。~ **コードの説明 [#j58cb272] 今回のデータは、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) **今回の例の結果（実行は2008年1月31日） [#j77991d1] ところどころに説明を挿入 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~ （上のほうででてきています） **参考にしたウェブページ [#fc6d3c44] -http://hosho.ees.hokudai.ac.jp/~kato/unix/R.html -http://hosho.ees.hokudai.ac.jp/~kato/seminar/020909/index.html -http://www.tuat.ac.jp/~ethology/cgi-bin/wiki.cgi?page=R%A4%C7%C3%B1%B2%F3%B5%A2 -http://www001.upp.so-net.ne.jp/ito-hi/stat/R1.html |Today:&counter(today); |Yesterday:&counter(yesterday); |Total:&counter(); since 01 February 2008|
