今回のRに関する記事はマジメです。役に立つ人には多少は役に立つかもしれません。
Rのデフォルト関数で、標準化回帰係数を出力してくれるものはないよう。
独立変数がすべて量的変数ならたとえば、scale() などを介すことで簡単に計算できます。
しかし問題になるのは、名義尺度などの変数です。Rの lm(), glm() では名義尺度が含まれていた場合、自動的にダミー変数化して計算してくれますが、標準化回帰係数を計算するうえでは問題が起こります。
たとえば、標準化回帰係数を出力するとうたう lm.beta()関数が {QuantPsyc}パッケージにありますが*1、僕の理解が間違っていなければ、プログラムにミスがあります。これを名義尺度が含まれる回帰分析に適用すると、おそらく計算がおかしくなるはず（理由は後述）。
lm.Beta()
とまれ修正版をつくってみました。
lm.Beta <- function(res){
b <- res$coefficients #非標準化回帰係数
ysd <- sd(res$model[,1]) #結果変数のSD
idv <- res$model[,-1]
if( class(idv) == "data.frame"){
N <- ncol(idv) }else{ N <- 1}
res.beta <- b[-1]
for(j in 1 : N ){
if( N > 1){xxx <- idv[,j]
}else{xxx <- idv}
if( class(xxx)!= "numeric"){ #もし数値じゃなければ以下をかます
for(i in 2:length(levels(xxx)) ){
dummy <- as.numeric(xxx == levels(xxx)[i]) # 1/0 の変数化
lab <- paste(colnames(idv)[j],levels(xxx)[i],sep="") #対象の変数
beta <- b[lab] * ( sd(dummy) / ysd ) #ベータの計算方法：2変数のSDの商に回帰係数をかける
res.beta[lab] <- beta
}
}else{ #数値だったら簡単にβを計算できる
lab <- paste(colnames(idv)[j],sep="") #対象の変数
beta <- b[lab] * ( sd(xxx) / ysd )
res.beta[lab] <- beta
}
}
res.beta
}
使い方
lm <- (y ~ x1 + x2 + ...) の結果を lm.Beta() の中に入れてください。
使用例
> #### y = x + z ... の重回帰分析
> y <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
> x <- c(2, 1, 4, 3, 6, 5, 8, 7, 10, 9)
> z <- factor(c("A","A","B","B","B","B","C","C","C","C"))
>
> #y, x, z をそれぞれスケール
> Y <- scale(y)
> X <- scale(x)
> Z.B <- scale(as.numeric(z=="B"))
> Z.C <- scale(as.numeric(z=="C"))
>
> #ふつうの重回帰
> res <- lm(y ~ x + z)
> #スケール化した重回帰
> res.scale <- lm(Y ~ X + Z.B + Z.C)
>
> res
Call:
lm(formula = y ~ x + z)
Coefficients:
(Intercept) x zB zC
0.7143 0.5238 1.4286 3.3333
> #lm.Beta()をつかって標準化回帰係数
> lm.Beta(res)
x zB zC
0.5238095 0.2436580 0.5685352
>
> #スケール化したデータの回帰係数
> res.scale$coef[2:4]
X Z.B Z.C
0.5238095 0.2436580 0.5685352
最後の数行に注目。
lm.Beta()による係数と、スケール化したデータの回帰係数が一致していることがわかります。
lm.beta() が間違う場合があるワケ
結論からいえば、名義尺度の変数を sd() 関数に直接ぶち込んでいるところ。
sd() 関数は、見た目からもわかるとおり標準偏差を出力する関数ですが、名義尺度（factor形式）が入れると勝手に数値に読み替えてしまいます。
たとえば、あ／い／う／え／お という5水準の変数も、以下のように
> aiueo <- factor(c("あ","い","う","え","お"))
> sd(aiueo)
[1] 1.581139なぜか計算できてしまいます。
この理由。sd()関数の中身を見てもらえれば話は早いんですが、sd()関数は、名義尺度がきたとき、自動的に 1, 2, ..., n というように番号を振って、それを計算してしまうからです。以下がそのイメージ
> hoge <- c(1,2,3,4,5) > sd(hoge) [1] 1.581139
