model1 <- lm(X ~ age + group)として，(ageは年齢で連続変数，groupはコントロール群と疾患群のfactorialです) ← factorial ではなく factor ですね
model2 <- lm(X ~ group + age)
> summary(model1)summary(model2)は，当然，summary(model1)と一致します．
Call:
lm(formula = X ~ age + group, data = ASL.ROI.affected.df)
Residuals:
Min 1Q Median 3Q Max
-6.8523 -3.1781 -0.6081 1.9931 8.7069
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.685257 4.527610 5.894 8.74e-07 ***
age -0.005839 0.063215 -0.092 0.927
groupnormal 3.458099 2.076886 1.665 0.104
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.285 on 37 degrees of freedom
Multiple R-squared: 0.1607, Adjusted R-squared: 0.1153
F-statistic: 3.541 on 2 and 37 DF, p-value: 0.03917
> anova(model1)つ まり，偏回帰係数表では有意な説明変数はなく，分散分析表では，モデルに加える順番によって，ageまたはgroupのいずれか一方（先に加えた方）が有 意になると言うことです．このような場合，どう解釈すべきなのでしょうか？都合の良いのを取ると言うのもありかもしれませんが（笑），この場合，どのパ ターンでも都合が悪くないので，余計に困ってます．groupのp値が小さいのでこれを取るべきなのかなとも考えています．
Analysis of Variance Table
Response: X
Df Sum Sq Mean Sq F value Pr(>F)
age 1 79.11 79.11 4.3096 0.0449 *
group 1 50.89 50.89 2.7724 0.1044
Residuals 37 679.24 18.36
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(model2)
Analysis of Variance Table
Response: X
Df Sum Sq Mean Sq F value Pr(>F)
group 1 129.85 129.85 7.0734 0.01149 *
age 1 0.16 0.16 0.0085 0.92690
Residuals 37 679.24 18.36
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(lm(X ~ age))こうしてみると，groupの方が，p値は小さいし，anovaでも残差平方和が小さいので，こちらを選ぶべきなのかなと言う気はしますが，よく分かりません．一応両方とも有意になってしまうので困ってしまいます．
Analysis of Variance Table
Response: X
Df Sum Sq Mean Sq F value Pr(>F)
age 1 79.11 79.11 4.1176 0.04949 *
Residuals 38 730.13 19.21
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(lm(X ~ group)
Analysis of Variance Table
Response: X
Df Sum Sq Mean Sq F value Pr(>F)
group 1 129.85 129.85 7.2629 0.01043 *
Residuals 38 679.39 17.88
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(lm(X ~ age), lm(X ~ group))
Analysis of Variance Table
Model 1: X ~ age
Model 2: X ~ group
Res.Df RSS Df Sum of Sq F Pr(>F)
1 38 730.13
2 38 679.39 0 50.74
No.10992 Re: lmでのanovaの結果の齟齬？ 【青木繁伸】 2009/10/03(Sat) 07:14
2変数だけで質問のようなことが生じているのなら，話は若干簡単。従属変数との単相関係数の大きい方を使えばよいでしょう。
他 でもコメントしたとおり，anova は変数編入順序によって結果が左右されるということ。anova の根底にあるのは，事前にある程度，変数の重要度の順序についての仮定があるということでしょう。anova を使ってこの変数の重要度を決めようとしているから，鶏が先か卵が先かという話になっているのではないでしょうか。重要な順に変数を選択しようとしている のですから，最初は変数を1個だけ含むモデルを作って，一番優れたモデル（重要な変数をからなるモデル）を決定し，次の段階で，その変数と別の1つの変数 を含むモデルを作って，一番優れたモデルを作って。あと，3つの変数を含むモデルを作って，，，，とやっていくのと根本的には同じですね。このようにドン ドン変数を追加していけば，最終的に重要な変数だけを含むモデルができそうですが，そうではありません。前の方で取り込まれた変数が，後の方になると役に 立たない変数になってしまうこともあります。このような方法に従えばまだましですが，単変量モデルで有意だった変数だけをあつめて多変量モデルをこうせい するという，誤った方策を採るということも，多く見られますが，多変量解析の本質を無視したやり方です。
別のスレッドの結論にも書きましたが，変数の重要度を決めるのは，データを分析した結果は参考にするのは当然ですが，それに依存してしまってはいけないということもあります。完全にデータに依存してよいと言うことなら，総当たりモデルで決定すれば宜しいかと思います。
No.10994 Re: lmでのanovaの結果の齟齬？ 【taipapa】 2009/10/03(Sat) 11:22
青木先生，コメント有り難うございます．
X ~ ageとX ~ groupのモデルの相関係数をそれぞれ見てみると，Estimate Std. Error t value Pr(>|t|)以上より，単相関係数の自乗は，ageが0.09776，groupが0.1605となって，groupを選ぶべきと言うことですね．なるほど．
(Intercept) 33.00527 2.52507 13.071 1.23e-15 ***
age -0.08561 0.04219 -2.029 0.0495 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.383 on 38 degrees of freedom
Multiple R-squared: 0.09776, Adjusted R-squared: 0.07402
F-statistic: 4.118 on 1 and 38 DF, p-value: 0.04949
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.2765 0.9455 27.792 <2e-16 ***
groupnormal 3.6035 1.3371 2.695 0.0104 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.228 on 38 degrees of freedom
Multiple R-squared: 0.1605, Adjusted R-squared: 0.1384
F-statistic: 7.263 on 1 and 38 DF, p-value: 0.01043
これって，論文に書く時，「年齢補正をしてもコントロール群のXと疾患群のXとの間には有意な差が存在した」としても良いのでしょうか？
と いうのは，この場合，（当然のことながら）上記のモデルでも，t.testでも，t値は2.695と同じになります．上記のようなモデルの検討の過程をへ ていれば，結果としてただのWelchのt検定と同じになった場合でも「年齢補正をした」と言っても良いのかどうかご教示いただければ幸いです．
＃ あっ，「← factorial ではなく factor ですね」，そのとおりです．m(_0_)m
No.10999 Re: lmでのanovaの結果の齟齬？ 【青木繁伸】 2009/10/03(Sat) 21:21
[年齢補正をしても]って，年齢補正はしていないでしょう
独立変数が2つなので，優劣を付けるのは単相関係数を見れば十分と言うだけです。
年齢補正したと言うことなら，年齢とgroupを同時に使った結果を使うべきでしょう。
独立変数が0/1変数1つだけの場合はt検定の結果，回帰直線の傾きの検定結果，回帰分析の分散分析の結果，独立変数と従属変数の相関係数の検定はすべて同じになります（同じP値になります）。
> 単相関係数の自乗は，ageが0.09776，groupが0.1605となって
単純に相関係数を取ればよいのに。。。
No.11002 Re: lmでのanovaの結果の齟齬？ 【taipapa】 2009/10/04(Sun) 00:34
そうすると，lm(X ~ age + group)の結果を用いて，「年齢補正をしてもコントロール群のXと疾患群のXとの間には有意な差が存在しなかった」と言うべきですね．別の方法による 結果も解析してみたところ，差がありませんので，ここでも差がないモデルの方が総合的に良いように思えて来ました．
重ねてアドバイス有り難うございました．
単相関係数は，cor.testやったら，groupが2群のみのせいか（？）エラーになるんですよね (ToT) でも，ageとgroupだとできてしまうから不思議？？？
No.11005 Re: lmでのanovaの結果の齟齬？ 【青木繁伸】 2009/10/04(Sun) 07:14
> 単相関係数は，cor.testやったら，groupが2群のみのせいか（？）エラーにな
どういうエラーが出るんでしょうか。少なくとも group が2群のみと言うのは原因ではないですね。単に相関係数を計算するなら，cor もあります。> X <- rnorm(10)
> g <- factor(sample(2, 10, replace=TRUE))
> cor.test(X, g)
Pearson's product-moment correlation
data: X and g
t = 1.3481, df = 8, p-value = 0.2145
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2734436 0.8339629
sample estimates:
cor
0.4302593
> cor(X, g)
[1] 0.4302593
No.11007 Re: lmでのanovaの結果の齟齬？ 【taipapa】 2009/10/04(Sun) 13:23
青木先生，何度も有り難うございます．
やり直してみると，順序によってできたりできなかったりすることが分かりました．
> cor.test(X, group)
ピアソンの積率相関係数
データ: X と group
t値 = 2.695, 自由度 = 38, P値 = 0.01043
対立仮説: 母相関は，0ではない
95 パーセント信頼区間: 0.1017649 0.6330865
標本推定値:
相関係数
0.4005757
上記は大丈夫ですが，順序を入れ替えると以下のようにエラーになります．
> cor.test(group, X)
ピアソンの積率相関係数
データ: group と X
t値 = NA, 自由度 = 38, P値 = NA
対立仮説: 母相関は，0ではない
95 パーセント信頼区間: NA NA
標本推定値:
相関係数
NA
Warning message:
In cor(x, y) : 強制変換により NA が生成されました
corも同じです．
help をみても，x, y numeric vectors of data values. x and y must have the same length.とあるだけです．そもそも回帰じゃなくて相関なので，順番は関係ないと思うのですが．当たり前かもしれませんが，この順番でも nonparametricだとエラーになりません．
> cor.test(group, X, method="s")
スピアマンの順位相関係数
データ: group と AX
S = 6550.551, P値 = 0.01403
対立仮説: 母相関（ρ）は，0ではない
標本推定値:
ρ
0.3855018
> cor.test(group, X, method="k")
ケンドールの順位相関係数
データ: group と X
z = 2.4075, P値 = 0.01606
対立仮説: 母相関（τ）は，0ではない
標本推定値:
τ
0.3186711
不思議です？？？
No.11008 Re: lmでのanovaの結果の齟齬？ 【知ったかぶり】 2009/10/04(Sun) 17:49
groupがnumericではないためでしょう．
cor.test(as.numeric(group),X)
で問題なく計算できるはずです．
No.11010 Re: lmでのanovaの結果の齟齬？ 【青木繁伸】 2009/10/04(Sun) 20:07
話が横道にそれてしまいましたが。。。
私の手元にある R システムのうち，cor.test が導入された最古のバージョンは Version 1.5.1 (2002-06-17) ですが，以下の通り，順序に関係なく動きます（もちろん，最新版 R-2.9.2 でも同じくちゃんと動きます（メッセージの日本語化のための章制作の print.htest を使っておられるようですが，それを使ってもちゃんと動きます）。
また，(g <- as.character(sample(2, 10, replace=TRUE)))でもちゃんと動きました（class() は当然 "character" です）。(g <- as.logixal(sample(0:1, 10, replace=TRUE))) も，同じく。
なお，この頃と今じゃ set.seed の初期化が同じではないようで，R-2.9.2 では，同じ set.seed(1) でも，異なったベクトルが得られる。R : Copyright 2002, The R Development Core Team
Version 1.5.1 (2002-06-17)
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.
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.
[Previously saved workspace restored]
> set.seed(1)
> (x <- rnorm(10))
[1] -0.97462574 -0.56075523 -2.18392278 0.08951152 0.96951953 -0.44478838
[7] -0.02287879 0.09960695 -0.35515387 0.54290910
> class(x)
NULL なんと，この頃のバージョンは rnorm() が生成するベクトルのクラスは NULL !!（今は numeric）
> (g <- factor(sample(2, 10, replace=TRUE)))
[1] 1 1 1 1 1 1 2 2 2 2
Levels: 1 2
> class(g)
[1] "factor" g は，まちがいなく factor
> cor(x, g)
[1] 0.3466235
> cor(g, x)
[1] 0.3466235 順序の如何に関わらず同じ答え
> cor.test(x, g)
Pearson's product-moment correlation
data: x and g
t = 1.0452, df = 8, p-value = 0.3265
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.3620087 0.8013587
sample estimates:
cor
0.3466235
> cor.test(g, x)
Pearson's product-moment correlation
data: g and x
t = 1.0452, df = 8, p-value = 0.3265
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.3620087 0.8013587
sample estimates:
cor
0.3466235 この場合も順序に関係なく動く
No.11013 Re: lmでのanovaの結果の齟齬？ 【知ったかぶり】 2009/10/04(Sun) 22:35
横道にそれてしまって申し訳ないのですが，ちょっと不可解なので．
g1 <- sample(c("1","2"), 10, replace=TRUE) の場合は，なんの問題もなし．
g2 <- sample(c("a","b"), 10, replace=TRUE) とすると，cor(X,g2)，cor(g2,X)ともにエラー．
g3 <- factor(sample(c("a","b"), 10, replace=TRUE)) とすると，cor(X,g3)はOK，cor(g3,X)はエラー．
g4 <- as.numeric(factor(sample(c("a","b"), 10, replace=TRUE))) これでcor(X,g4)，cor(g4,X)ともにOK.
ダブルクォートで囲われていても「数字」であるか「文字」であるかで異なるようです．
No.11015 Re: lmでのanovaの結果の齟齬？ 【青木繁伸】 2009/10/04(Sun) 23:20
なるほど。そこまではやってみなかったなあ。質問者が状況を呈示してくれていると，こんな無駄なスレッドは生じなかったのに。g2 <- sample(c("a","b"), 10, replace=TRUE) とすると，cor(X,g2)，cor(g2,X)ともにエラー．なるほど，確かにそのようですね。R も，全ての場合にちゃんと対応するわけには行かないということですね。
g3 <- factor(sample(c("a","b"), 10, replace=TRUE)) とすると，cor(X,g3)はOK，cor(g3,X)はエラー．
あ らゆる場合に対応可能なようにユーザ側で，知ったかぶりさんの言うように，as.numeric(ほげほげ) と自衛策を講じるというのが正解。「そんなことぐらい R でやってよ」というのは，ユーザのわがまま。適当に解釈して実行してしまうと，不適切なことも生じるだろうということかな（どのような不適切なことかはさ ておき）。
No.11016 Re: lmでのanovaの結果の齟齬？ 【taipapa】 2009/10/05(Mon) 00:34
青木先生，知ったかぶりさん，有り難うございました．
質問者が知らない間に長いスレッドになって問題が解決していました．お二人には申し訳ありませんでしたが，私には勉強になりました．(^^;;;
cor.test(as.numeric(group), X)で確かにうまくいきます．
factorかどうかだけでなく，numericかどうかまで注意しないといけないのですね．
groupは2群と書きましたが，0,1ではなく，negative, positiveだったのです．
数字か文字か，それが問題だ...
お二方，重ねて有り難うございました．
No.11022 Re: lmでのanovaの結果の齟齬？ 【にゃんちゅう】 2009/10/06(Tue) 13:46
>上記は大丈夫ですが，順序を入れ替えると以下のようにエラーになります．
>> cor.test(group, X)
>ピアソンの積率相関係数
一応解決したようですが。
私は最初にこれを読んだときに最初の変数は尺度水準のチェックが入っているが，第2変数のチェックが入っていない。と思ったのですが大違いですね。
第2変数は文字変数を数値に変換して処理しているということですね。なんでこんなめんどうなことを第2変数だけ行うことにしたのでしょうね。
No.11023 Re: lmでのanovaの結果の齟齬？ 【takahashi】 2009/10/06(Tue) 15:16
正しくは，最初の変数だけas.vectorされている，ということです。
> as.numeric(factor("a"))
[1] 1
> as.numeric(as.vector(factor("a")))
[1] NA
Warning message:
強制変換により NA が生成されました
この違いが原因です。
やっぱりバグっぽいですね。少なくとも開発者が意図した動作じゃないように思います。
● 「統計学関連なんでもあり」の過去ログ--- 043 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る