重回帰分析で標準化回帰係数βを出力するR関数 が紹介されている
ちょっと，見にくいところと余分なことがあるので，修正してみた（代入記号に = を使うのは，不本意だが，許して欲しい）
lm.Beta2 = function(res) {
ysd = sd(res$model[, 1]) # 結果変数のSD
idv = res$model[, -1, drop = FALSE]
N = ncol(idv)
res.beta = sd = res$coefficients[-1]
for (j in 1:N) {
xxx = idv[, j]
if (class(xxx) == "factor") {
for (i in 2:nlevels(xxx)) {
lab = paste(colnames(idv)[j], levels(xxx)[i], sep = "") # 対象の変数
dummy = as.integer(xxx) == i # 1/0 の変数化
sd[lab] = sd(dummy)
}
} else { # 数値だったら簡単にβを計算できる
lab = colnames(idv)[j] # 対象の変数
sd[lab] = sd(xxx)
}
}
res.beta * sd / ysd
}
R の他の関数のように，下請けの関数を組み合わせて目的を達成するという方針でやるならば，以下のように簡単に書くことができる（1行にすることもできるが，そこまでは...）
lm 関数が返すオブジェクト object に対して，
lm.Beta3 = function(object) {
d = model.matrix(object$terms, eval(object$model, parent.frame()))
object$coefficients[-1] * apply(d[, -1, drop=FALSE], 2, sd) / sd(object$model[,1])
}
> lm.Beta3(res)
x zB zC
0.5238095 0.2436580 0.5685352
楽しい，ハッキングだった(^_^)
ありがとうございました
