前回の
glmnetでロジスティック回帰の標準回帰係数を計算する - Qiita
を元に、今度はポワソン回帰をした時の標準回帰係数について調べてみました。
前回同様、データとしてbirthwtを使います。(1列目は削除)
R
data(birthwt)
df<-birthwt[,-1]
同様に説明変数に、生のデータと標準化したデータの2つを用意します。
R
# 目的変数
dy <- df[, 9]
# 生データ
dx <- data.matrix(df[, -9])
# 標準化したデータ
dx_std<-data.matrix(scale(df[,-9]))
目的変数としてbwtをつかっていいのか?という議論はココでは無視しますw
まずは、生のデータを使い、glmnet()の中でデータを標準化standardize=TRUEさせた場合。
R
require(glmnet)
set.seed(123)
cvfit <- cv.glmnet(x = dx, y = dy, family = "poisson", standardize = TRUE)
coef(cvfit, s = "lambda.min")
9 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) 8.0390164625
age .
lwt 0.0008950326
race -0.0520953345
smoke -0.0999303805
ptl -0.0098804539
ht -0.1635103496
ui -0.1671385398
ftv .
次に、説明変数を事前に標準化して、glmnet()の中でデータを標準化させないstandardize=FALSE場合
R
set.seed(123)
cvfit_std <- cv.glmnet(x = dx_std, y = dy, family = "poisson", standardize=FALSE)
coef(cvfit, s = "lambda.min")
9 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) 7.982804430
age .
lwt 0.027369544
race -0.047841346
smoke -0.048904913
ptl -0.004874442
ht -0.039977228
ui -0.059533132
ftv .
とロジスティック回帰の時と同じような結果に。
同じglmnet()なので、係数もロジスティック回帰と同様にagresti変換しているんじゃないかな～？、ということで検証。
R
set.seed(123)
cvfit <- cv.glmnet(x = dx, y = dy, family = "poisson", standardize = TRUE)
# 標準化されていない回帰係数
cf <- coef(cvfit, s = "lambda.min")[,1]
# 各変数の標準偏差
sds <- sapply(as.data.frame(dx), sd)
# 各変数の平均
mus <- sapply(as.data.frame(dx), mean)
# 回帰係数を標準化する
cf_std <- cf * c(1, c(sds[names(cf)][-1]))
# 切片を標準化する
cf_std[1] <- cf[1] + sum(cf[-1] * mus[names(cf)][-1])
as.data.frame(cf_std)
cf_std
(Intercept) 7.982804430
age 0.000000000
lwt 0.027369544
race -0.047841346
smoke -0.048904913
ptl -0.004874442
ht -0.039977228
ui -0.059533132
ftv 0.000000000
と、読みが大当たり！標準化したデータによる回帰係数と同じ値になりました。