glmnet에 의한 ridge 회귀의 복원

Ridge 계산에 관한 주석에 따라 glmnet에 의한 ridge 회귀 결과를 직접 계산하여 확인해 본다. 데이터는 여기 참고.

주석 내용을 반복하면 다음과 같다.

glmnet에 standardize = TRUE 옵션(디폴트)을 주고 계산하는 ridge 추정값은 다음과 같이 구할 수 있다. ① x̃_ij^* = (x_ij−x̄_j)/σ̂_j로 표준화한다. 여기서 x̄_j는 x_ij의 표본평균, σ̂_j는 x_ij의 표본평균편차(표본분산 계산 시 n으로 나눔). 표본평균을 빼는 것은 절편을 처리하기 위함이며, σ̂_j로 나누는 것은 변수 표준화. ② β̂^* = (n⁻¹X̃^*′X̃^*+λ^*I)⁻¹n⁻¹X̃^*′y를 구한다(X̃^*는 x̃_ij^*의 행렬, λ^* = λ/σ̂_y). 이것은 ‘베타계수’ β_jσ̂_j들의 추정값 벡터이다. ③ 표준화를 되돌리기 위해 β̂_j = β̂_j^*/σ̂_j를 계산하면 이것이 최종적인 기울기 계수의 ridge 추정값. ④ 절편 추정값은 앞과 동일하게 α̂ = ȳ − β̂₁x̄₁ − ⋯ − β̂_px̄_p로 계산.

`glmnet`에 의한 ridge 회귀

다음 glmnet에 의한 결과를 보라.

ridge <- glmnet(X, Y, alpha = 0)
(lamb <- ridge$lambda[50])  # an arbitrary choice
# [1] 3485.563
coef(ridge, s = lamb)
# 20 x 1 sparse Matrix of class "dgCMatrix"
#                         1
# (Intercept) 1349.96474879
# grdp          -1.46733791
# regpop       -70.45896296
# popgrowth     -2.60572536
# eq5d          -2.24328091
# deaths        -0.01593032
# drink         -2.97176050
# hdrink         1.03296288
# smoke          1.90460780
# aged           2.90703405
# divorce       -8.98595269
# medrate        0.43610941
# gcomp         -0.90484486
# vehipc        97.75098402
# accpv         -1.13114311
# dumppc        -4.58445915
# stratio       -3.16084850
# deathrate      0.07065207
# pctmale       -3.94300918
# accpc          1.88459406

이 주어진 λ (lamb)을 이용한 ridge 추정값을 행렬연산에 의하여 똑같이 구해 보는 것이 목적이다.

행렬 연산에 의한 ridge 구현

먼저, glmnet에서 말하는 “표준편차”들은 n − 1이 아닌 n으로 나누어 구한다. 이를 위해 sd0라는 함수를 만든다.

sd0 <- function(x) sqrt(mean((x-mean(x))^2))

## X, Y defined in data preparation page
sig.Y <- sd0(Y)
sig.X <- apply(X, 2, sd0)

준비: λ의 변환

(lamb.star <- lamb/sig.Y)
# [1] 10.35528

glmnet에서 사용한 λ 값 3485.563은 행렬 연산을 통한 수동 연산에서 10.35528에 해당한다.

① X변수들의 표준화

Xs <- apply(X, 2, function(x) (x-mean(x))/sd0(x))

② β̂^*의 계산

n <- nrow(Xs)
bhat.star <- solve(1/n*crossprod(Xs) + lamb.star*diag(ncol(Xs)), 1/n*crossprod(Xs, Y))

③ β̂의 복원

(bhat <- bhat.star / sig.X)
#                   [,1]
# grdp       -1.46731962
# regpop    -70.45811386
# popgrowth  -2.60571609
# eq5d       -2.24326913
# deaths     -0.01593022
# drink      -2.97175109
# hdrink      1.03295937
# smoke       1.90460738
# aged        2.90703561
# divorce    -8.98595524
# medrate     0.43611098
# gcomp      -0.90484667
# vehipc     97.75111247
# accpv      -1.13114176
# dumppc     -4.58445394
# stratio    -3.16085391
# deathrate   0.07065214
# pctmale    -3.94300686
# accpc       1.88459496

이 결과는 앞에서 glmnet으로 구한 결과와 동일함을 확인할 수 있다.

④ 절편 추정값 확인

mean(Y) - sum(colMeans(X)*bhat)
# [1] 1349.963

대안적인 계산 방법

주석에 다음과 같이 쓰여 있기도 하다.

참고로, 위 ①∼③의 β̂은 D = diag(σ̂₁²,…,σ̂_p²)라 할 때 (n⁻¹X̃′X̃+λ^*D)⁻¹n⁻¹X̃′y으로 계산한 것과 동일하다.

이를 확인해 보자.

Xd <- apply(X, 2, function(x) x - mean(x))
solve(1/n*crossprod(Xd) + lamb.star*diag(sig.X^2), 1/n*crossprod(Xd,Y))
#                   [,1]
# grdp       -1.46731962
# regpop    -70.45811386
# popgrowth  -2.60571609
# eq5d       -2.24326913
# deaths     -0.01593022
# drink      -2.97175109
# hdrink      1.03295937
# smoke       1.90460738
# aged        2.90703561
# divorce    -8.98595524
# medrate     0.43611098
# gcomp      -0.90484667
# vehipc     97.75111247
# accpv      -1.13114176
# dumppc     -4.58445394
# stratio    -3.16085391
# deathrate   0.07065214
# pctmale    -3.94300686
# accpc       1.88459496

결과는 glmnet과 동일하다.

[Ridge 계산에 관한 주석으로 돌아가기]

glmnet에 의한 ridge 회귀