重回帰分析で標準化回帰係数βを出力するR関数

今回の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

SPSSではどう計算しているか教えてください

ちなみにけっこう気になるのは、SPSSがこれをどう計算しているのかという点です。

私はSPSSを持っていないので、上記の計算結果が一致するかどうかわかりません。

どなたか優しい方、試してはもらえませんでしょうか。よろしくお願いいたします。
各変数はこんな感じです。

Y

1, 2, 3, 4, 5, 6, 7, 8, 9, 10

x

2, 1, 4, 3, 6, 5, 8, 7, 10, 9

z

A, A, B, B, B, B, C, C, C, C

*1:http://www.inside-r.org/packages/cran/QuantPsyc/docs/lm.beta