No.13923 GLM等において要因の影響の強さではなく、要因が与える分散の大きさを比較したいときの方法  【Sai】 2010/12/06(Mon) 16:58

いつもお世話になっております。
過去ログ等キーワード検索をしましたが,見つからなかったので質問させてください。

一般的に応答変数をYとして,説明変数をX, Zとするような重回帰分析はY=aX + bZ + cと書けると思います。この際,a, b, cは推定されるパラメータです。

こ の時に,a, bの値によって,「YがXあるいはZから独立に受ける影響の強さ」,は定量化することが出来ると思います。ただし,今私が知りたいのは,係数の値ではなく て,「Yの全平均周りのバラツキを生んだ原因はXとZどちらがより大きいのか」,ということを知りたいと思っています。

具体的には,X, Zがカテゴリカル変数でそれぞれ2水準と3水準を持つとします。つまり二元配置の分散分析のような格好です。このときに,X方向のバラツキとZ方向のバラツキというのは,どのように定量化すれば良いでしょうか?以下にサンプルプログラムを示します。
###水準の設定###
X <- c("X1", "X2")
Z <- c("Z1", "Z2", "Z3")
gri <- expand.grid(X=X, Z=Z)

###繰り返し数###
n <- 50

###0, 1の行列を作ってYを作成する###
gri[(ncol(gri)+1):(length(c(X, Z))+ncol(gri))] <- matrix(0, nrow(gri), length(c(X, Z)))

gri[which(gri[, 1]==X[1]), 3] <- 1
gri[which(gri[, 1]==X[2]), 4] <- 1
gri[which(gri[, 2]==Z[1]), 5] <- 1
gri[which(gri[, 2]==Z[2]), 6] <- 1
gri[which(gri[, 2]==Z[3]), 7] <- 1

###繰り返し数だけ取る###
gri2 <- NULL
for(i in 1 : n) {
gri2 <- rbind(gri2, gri)
}

###真のパラメータの設定###
###XよりもZの方が値も大きく,またばらついている###
x <- c(1, 2)
z <- c(1, 20, 300)

###誤差を入れて0, 1行列を乱す###
gri_2_1 <- matrix(as.matrix(gri2[, 3:7]),
nrow(gri)*n, length(c(X, Z)))*
runif(nrow(gri)*n*length(c(X, Z)), 0.1, 1.9)

###もう一度誤差を入れてYを作る###
gri2[, length(c(X, Z)) + 3] <- gri_2_1 %*% c(x, z)+rnorm(nrow(gri2), 0, 10)

###大きくバラついているように見える###
plot(gri2[, 8])

###回帰###
res <- lm(gri2[, 8] ~ gri2$X + gri2$Z)

###このときにXとZの処理の,どちらがYの(全平均値からの)バラツキをより生み出していたのか知りたい###
summary(res)
全分散=群間分散+群内分散に分かれることを利用して,とも思ったのですが,「Zの影響を取り除いたXのYのバラツキに与える影響」,のような関係を見たいと思っているので,その場合の上手い対処法がわかりませんでした。

何か根本的に勘違いをしているのかもしれません。
どなたか,わかる方がいらっしゃいましたら,ご教授頂けると幸いです。

No.13924 Re: GLM等において要因の影響の強さではなく,要因が与える分散の大きさを比較したいときの方法  【青木繁伸】 2010/12/06(Mon) 18:26

上に引き続いて,lm の結果 (res) に対して anova(res) をやるのとは違います?
> anova(res)
Analysis of Variance Table

Response: gri2[, 8]
Df Sum Sq Mean Sq F value Pr(>F)
gri2$X 1 551 551 0.062 0.8035
gri2$Z 2 5762409 2881204 323.977 <2e-16 ***
Residuals 296 2632397 8893
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
交互作用項でもなんでもモデルに加えておけば,それぞれの平方和(または平均平方)が表示されますよね?

No.13925 Re: GLM等において要因の影響の強さではなく,要因が与える分散の大きさを比較したいときの方法  【Sai】 2010/12/06(Mon) 19:10

早速の返信ありがとうございました。
ただ,まだ疑問が解決しません。anovaで出力されるSum sqやそれを自由度で割ったものはZの影響を取り除いたXのYのバラツキに与える影響を定量的に表したものなのでしょうか?

青木先生がせっかく追試してくださったのに恐縮ですが,再現性のためにset.seed(1)として再度上記プログラムを走らせますと,以下のような結果になります。
> anova(res)
Analysis of Variance Table

Response: gri2[, 8]
Df Sum Sq Mean Sq F value Pr(>F)
gri2$X 1 2554 2554 0.3014 0.5834
gri2$Z 2 5905512 2952756 348.4554 <2e-16 ***
Residuals 296 2508257 8474
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘
こ のとき,Xで説明した(逐次)平方和は2554,Zは5905512。従って,Zの方がXのそれよりもYの平均値までの距離の総和が大きいことを意味して いる。そのため,Zの処理ほうがXの処理よりもYに与えるバラツキの程度は大きい。・・・と解釈して良いのでしょうか?

疑問点は2つです。

1. Zの影響はXの影響を取り除いたものとなっているか?
2. ZがXと比べYに与えるバラツキは「このくらい大きい」と定量的に言うことは出来るか?出来るとしたら,どこの値を統計量として用いれば良いのか?

個人的には1.は調整平方和を使えば良いのかと考えました。つまり,
> library(car)
> Anova(res)
Anova Table (Type II tests)

Response: gri2[, 8]
Sum Sq Df F value Pr(>F)
gri2$X 2554 1 0.3014 0.5834
gri2$Z 5905512 2 348.4554 <2e-16 ***
Residuals 2508257 296
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1
です。偶然(?)anovaの結果と一緒なのは若干不思議な気がしますが,投入の順番では値が変わらないということでしょうね。ここで出力されたSum sqで議論すれば,Zの影響を取り除いたXという風になっているでしょうか?

2.に関しては,単に5905512-2554で良いのか,と考えたもののそれは良くない気がします。今回はサンプル数は等しいですが,もしXとZでサンプル数に違いがあったらそのような単純な引き算ではダメかな,と。

まだ何か勘違いをしているかもしれません。
ご教示してくださると幸いです。

No.13927 Re: GLM等において要因の影響の強さではなく,要因が与える分散の大きさを比較したいときの方法  【青木繁伸】 2010/12/06(Mon) 20:54

> Zの影響はXの影響を取り除いたものとなっているか?

平方和については,全体の平方和は総ての平方和の和になっているので,それぞれの和は独立というか,他の平方和に影響されないと言うことでしょう?それが,分散分析というものでは?よくわかりませんけど。

No.13928 Re: GLM等において要因の影響の強さではなく,要因が与える分散の大きさを比較したいときの方法  【Sai】 2010/12/06(Mon) 21:31

返信および修正ありがとうございました。

>それぞれの和は独立というか,他の平方和に影響されないと言うこと

この部分ですが独立変数間に相関があった場合,Sum sqを鵜呑みには出来ません。以下にサンプルコードを示します。カテゴリカルの例が作れれば良かったのですが上手く出来ませんでした。連続変数でご容赦ください。

> set.seed(1)
> library(DAAG)
> x <- seq(1, 100)
> z <- x*1.1 + rnorm(100)
> a <- 1
> b <- 2
> c <- 3
> y <- a*x + b*z + c + rnorm(length(x), 0, 10)
>
> res <- lm(y ~ x + z)
> anova(res)
Analysis of Variance Table

Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x 1 857112 857112 9158.3481 < 2e-16 ***
z 1 317 317 3.3899 0.06866 .
Residuals 97 9078 94
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> Anova(res)
Anova Table (Type II tests)

Response: y
Sum Sq Df F value Pr(>F)
x 68.1 1 0.7273 0.39585
z 317.2 1 3.3899 0.06866 .
Residuals 9078.0 97
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> vif(res)
x z
1262.583 1262.583
こういう場合がありえるので,普通にやっても独立で計算はしてくれません。また,anovaではなく,Anova()を使えば良いのかもしれませんが,そこから計算される調整平方和でYの分散への寄与度をどのように議論するのかの理論も知らないので,困っております。

まだ何か勘違いしているのかもしれません。
その点,ご教示いただければ幸いに思います。

No.13929 Re: GLM等において要因の影響の強さではなく,要因が与える分散の大きさを比較したいときの方法  【青木繁伸】 2010/12/06(Mon) 21:36

> スペースを入れても上手く表を直すことが出来ませんでした。申し訳ありません。

自分の投稿がどのように編集されているかを見れば,どのように編集すればよいかはわかるかも知れませんね(html ファイルの基本的文法--タグの使い方です)。

No.13931 Re: GLM等において要因の影響の強さではなく,要因が与える分散の大きさを比較したいときの方法  【Sai】 2010/12/06(Mon) 22:01

なるほど,<pre>というのを使うのですか。以後気を付けます。

No.13937 Re: GLM等において要因の影響の強さではなく,要因が与える分散の大きさを比較したいときの方法  【kai】 2010/12/07(Tue) 15:15

Rでどう解析するかは分かりませんが,要因の効果を分散として定量化する解析,すなわち,変量効果(Random effect)として解析すれば,Yのバラツキのうち,因子X,因子Zでそれぞれどれぐらい寄与しているのか,分散の値で把握できます.

二元配置の分散分析程度であれば,分散分析表の平均平方から各因子の分散成分が計算できます.(平均平方は分散成分の和になっています)

No.13939 Re: GLM等において要因の影響の強さではなく,要因が与える分散の大きさを比較したいときの方法  【Sai】 2010/12/07(Tue) 20:02

kaiさん,返信ありがとうございました。
つまり平均平方同士の直接比較で良いということでしょうか?

ちなみにRではglmmMLやlmerというパッケージを用いれば出来るようです。

No.13945 Re: GLM等において要因の影響の強さではなく,要因が与える分散の大きさを比較したいときの方法  【kai】 2010/12/09(Thu) 12:19

平均平方は複数の分散成分が混ざった状態なので,純粋な比較で言えば分散成分どうしを比較した方がよいと思います.

No.13946 Re: GLM等において要因の影響の強さではなく,要因が与える分散の大きさを比較したいときの方法  【Sai】 2010/12/09(Thu) 21:29

kaiさん返信ありがとうございました。
glmmやlmerだと上手く出来ませんでしたので,自力で平均 平方から残差を引いて,繰り返し数で割る,という作業で分散成分を抽出しました。一応,「生物統計学入門」山田作太郎・北田修一さんの本の変量モデルの項 を参考に組みましたが,勘違いをしているような気もしています(0以下の分散の解釈など)。もし何かお気づきの点があれば,ご指摘頂けると幸いです。以下 がサンプルプログラムです。
> ###説明されなかった群内分散を説明された平均平方から引く。そのときに繰り返し数で割る###
> sigma <- (anova(res)[, 3]-anova(res)[length(anova(res)[, 3]), 3])/n

> ###分散がマイナスということは帰無仮説を棄却出来なかったということなので,正のものだけ取り出す###
> Sigma <- sigma[sigma>0]

> ###今回は一つだけ残ったが,本来はここに残ったものの大きさの比較を行う###
> Sigma
[1] 57603.91

● 「統計学関連なんでもあり」の過去ログ--- 044 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る