No.10265 Re: 2つの標準偏回帰係数の有意差の検定 【青木繁伸】 2009/07/04(Sat) 18:33
以下のような，プログラムを書いてみました。
lm 関数の戻り値と，偏回帰係数の差の検定をしたい2個の従属変数名を渡すと，結果を返してくれます。プログラムがどのように作られているか，読んでみてください。func <- function(a, i, j)
{
x <- summary(a)$coefficients[-1,1:2]
z <- abs(x[i, 1]-x[j, 1])/sqrt(x[i, 2]^2+x[j, 2]^2)
p <- pnorm(z, lower.tail=FALSE)*2
return(list(coeff=x[c(i, j),], test=c("Z-value"=z, "P-value"=p)))
}
set.seed(123)
x <- data.frame(matrix(rnorm(100), 20))
colnames(x) <- c(paste("X", 1:4, sep=""), "Y")
a <- lm(Y~., x)
summary(a)
func(a, "X1", "X3")
func(a, "X3", "X4")
func(a, "X3", "X2")
実行例は以下の通り> func <- function(a, i, j)
+ {
+ x <- summary(a)$coefficients[-1,1:2]
+ z <- abs(x[i, 1]-x[j, 1])/sqrt(x[i, 2]^2+x[j, 2]^2)
+ p <- pnorm(z, lower.tail=FALSE)*2
+ return(list(coeff=x[c(i, j),], test=c("Z-value"=z, "P-value"=p)))
+ }
>
> set.seed(123)
> x <- data.frame(matrix(rnorm(100), 20))
> colnames(x) <- c(paste("X", 1:4, sep=""), "Y")
> a <- lm(Y~., x)
> summary(a)
Call:
lm(formula = Y ~ ., data = x)
Residuals:
Min 1Q Median 3Q Max
-1.3224 -0.3175 -0.1559 0.5211 1.1749
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4380 0.1651 2.654 0.0181 *
X1 -0.1783 0.1726 -1.033 0.3180
X2 0.4545 0.2071 2.194 0.0444 *
X3 -0.3975 0.1764 -2.254 0.0396 *
X4 -0.2330 0.1782 -1.308 0.2106
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7177 on 15 degrees of freedom
Multiple R-squared: 0.4082, Adjusted R-squared: 0.2504
F-statistic: 2.587 on 4 and 15 DF, p-value: 0.07943
> func(a, "X1", "X3")
$coeff
Estimate Std. Error
X1 -0.1783128 0.1726378
X3 -0.3974719 0.1763580
$test
Z-value P-value
0.8880343 0.3745223
> func(a, "X3", "X4")
$coeff
Estimate Std. Error
X3 -0.3974719 0.1763580
X4 -0.2330306 0.1781542
$test
Z-value P-value
0.6559777 0.5118384
> func(a, "X3", "X2")
$coeff
Estimate Std. Error
X3 -0.3974719 0.1763580
X2 0.4544601 0.2070967
$test
Z-value P-value
3.131951375 0.001736486
No.10421 Re: 2つの標準偏回帰係数の有意差の検定 【青木繁伸】 2009/07/21(Tue) 11:59
> 最後の質問にできればどなたかお答願えないでしょうか？
普通すぎて，何処にそれが書いてあったかなどは忘却の彼方。探してみれば幾つも（と言うほどでもないが）見つかる。
例
えば，R の asympTest パッケージ。これは，添付図に示すようなことを目的とするもの。まあ，答えを言えば，"Asymptotic
parametoric test" 漸近的パラメトリック検定（Asymptotic test 漸近検定でよいと思う）ということ。
な
お，本来この漸近検定は二つの統計量が独立であることが要求されます。質問の場合のような二つの偏回帰係数も独立ではないかもしれませんが，まあ，漸近近
似ということで許されるようでもあります（いくつかの時点における生存率の差の検定にも使われるので）。しかし，偏回帰係数の検定において注意すべきこと
は，二つの変数が同じ土俵で比較できるものでなければならないと言うことです。その1として，意味的に比較できるものであること。たとえれば，数学と英語
の試験の平均値に差があるかどうか調べても意味がないだろうというようなこと。その2として，測定単位が同じであること（さらに，分散がほぼ同じであるこ
と）。有るデータセットで，二つの偏回帰係数の漸近検定を行った結果と，一方の変数を100倍（あるいは100分の1倍）したデータセットに対して同じ検
定を行うと，結果として得られるZ値が異なることがわかるでしょう。当然，その結果としてP値も異なり，検定結果（有意か有意でないか）も変わることがあ
るでしょう。
そういう点に注意して，この検定を行ってください。
以下は，Rを使っての例です。> library(asympTest)
> set.seed(123)
> gendat1 <- function(n, mu=0, sigma=1) # 指定された平均値と標準偏差を持つデータベクトルを生成
+ {
+ x <- rnorm(n) # n 個の正規乱数を発生
+ return((x-mean(x))/sd(x)*sigma+mu) # 基準化する
+ }
> sim <- function(n1, m1, sd1, n2, m2, sd2)
+ {
+ x1 <- gendat1(n1, m1, sd1)
+ x2 <- gendat1(n2, m2, sd2)
+ print(asymp.test(x1, x2)) # 漸近検定（statistic が z 統計量）
+ z <- (mean(x1)-mean(x2))/seDMean(x1, x2) # 実際にやっていること（z 値の計算）
+ cat(sprintf("z value = %f statistic と同じになる\n", z))
+ cat(sprintf("seDMean = %f\n", seDMean(x1, x2))) # 推定される標準偏差
+ # > asympTest:::seDMean.default # seDMean のやっていること
+ # function (x, y, rho = 1, ...)
+ # {
+ # sqrt(var(x)/length(x) + rho^2 * var(y)/length(y))
+ # }
+ #
+ # つまり，デフォルトの rho=1 のときは，以下と同じ
+ se1 <- sd1/sqrt(n1)
+ se2 <- sd2/sqrt(n2)
+ se.est <- sqrt(se1^2+se2^2) # この値と，seDMean の値が同じになる
+ cat(sprintf("se.est = %f seDMean と同じになる\n", se.est))
+ # ちなみに，等分散を仮定する t.test の結果は
+ print(t.test(x1, x2, var.equal=TRUE))
+ }
>
> n1 <- 50 # 第1群のサンプルサイズ，平均値，標準偏差
> m1 <- 50
> sd1 <- 10
> n2 <- 70 # 第2群のサンプルサイズ，平均値，標準偏差
> m2 <- 52
> sd2 <- 11
> sim(n1, m1, sd1, n2, m2, sd2)
Two-sample asymptotic difference of means test
data: x1 and x2
statistic = -1.0358, p-value = 0.3003
alternative hypothesis: true difference of means is not equal to 0
95 percent confidence interval:
-5.784594 1.784594
sample estimates:
difference of means
-2
z value = -1.035759 statistic と同じになる
seDMean = 1.930951
se.est = 1.930951 seDMean と同じになる
Two Sample t-test
data: x1 and x2
t = -1.0193, df = 118, p-value = 0.3101
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-5.885367 1.885367
sample estimates:
mean of x mean of y
50 52