Glmnet Vignette
Trevor Hastie and Junyang Qian
Stanford September 13, 2016
Introduction
Installation
Quick Start
Linear Regression
Logistic Regression
Poisson Models
Cox Models
Sparse Matrices
Appendix 0: Convergence Criteria
Appendix 1: Internal Parameters
Appendix 2: Comparison with other Packages
Installation
Quick Start
Linear Regression
Logistic Regression
Poisson Models
Cox Models
Sparse Matrices
Appendix 0: Convergence Criteria
Appendix 1: Internal Parameters
Appendix 2: Comparison with other Packages
선형 회귀 (Linear Regression)
선형 회귀(Linear regression)로 여기서 2개의 모델 군을 언급한다. 하나는 가우시안(gaussian) 이다, 가우시안 군(Gaussian family), 그리고 다른 하나는 mgaussian 이다, 다중 응답 가우시안 군. 우리는 먼저 보낼의 가우시안을 논의하고 그 후에 다중 응답 하나를 논의.
Gaussian Family
gaussian은 함수 glmnet에서 기본 family 옵션이다. 우리가 관측 xi ∈ Rp 가지고 그리고 응답 yi ∈ Rp, i = 1, ... , N 을 가정하자. 가우시안 군의 목적 함수는
λ ≥ 0 은 복잡도 변수 그리고 0 ≤ α ≤ 1은 ridge ( α = 0) 그리고 lasso (α = 1) 간의 절충(compromise)인 경우.
좌표 강하(Coordinate descent)은 문제를 해결하기 위해 적용되었다. 특히, 우리가 현재 추정 β0~ 그리고 βl~ ∀j ∈ 1, ], ..., p를 가정하자. βj = βj~ 에서 기울기 계산으로 그리고 단순 계산, 갱신은
, 그리고 S(z, γ)는 sign(z)(|z| - γ)+를 가지는 soft-thresholding 연산자인 경우.
위 공식은 x 변수가 단위 변동을 가지기 위해 표준화 되었을 때 적용한다(기본); 그것은 그들이 표준화가 아닐 때 약간 더 복잡하다. "family=gaussian"을 위해, glmnet은 그것의 람다 연속(lambda sequence)를 계산하기 전에 단위 변동을 가지게 y를 표준화함을 주목하라(그리고 그 다음 결과 계수를 비표준화한다); 만약 당신이 다른 소프트웨어를 가지고 결과를 재생산/비교를 원한다면, 먼저 표준화된 y를 제공하는 것이 더 좋다 ("1/N" 변동 식을 사용).
glmnet은 사용자가 적합을 정의하기 위한 다양한 옵션을 제공한다. 우리는 여서거 약간 일반적으로 사용된 옵션을 소개한다 그리고 그들은 glmnet 함수에서 지정할 수 있다.
- alpha는 탄력적-순(elastic-net) 혼합 변수 α 이다, 범위 α ∈ [0, 1]를 가짐. α = 1은 lasso(기본) 이다 그리고 lasso = 0은 ridge 이다.
- weights는 관측 가중치이다. 기본은 각 관측에서 1 이다. (주의: glmnet은 N의 합으로 가중치를 재척도한다, N은 표본 크기.)
- nlambda는 연속으로 λ 값의 수이다. 기본값 100.
- lambda를 제공할 수 있다, 그러나 전형적이지 않다 그러나 프로그램은 연속으로 생성한다. 자동으로 생성될 때, λ 연속은 lambda.max 그리고 lambda.min.ratio 으로 결정된다. 나중에 lambda.max에서 생성된 λ 연속(λ sequence)의 가장 작은 값(lambda.min이라 말함)의 비이다. 프로그램은 그러면 lambda.max에서 아래로 lambda.min까지의 로그 척도에 생성된 nlambda 값 선형. lambda.max는 주어지지 않는다, 그러나 입력 x와 y로 쉽게 구할 수 있다; 그것은 모든 계수가 0인 lambda에서 가장 작은 값이다. alpha=0 (ridge) 에서 lambda.max는 ∞가 될 것이다; 그러므로 이 경우에 우리는 0에 가까운 alpha의 가장 작은 갓에 상응하는 값을 선택한다.)
- standardize는 x 변수 표준화를 위한 논리 표시이다, 모델을 연속으로 적합하기에 앞서, 계수는 항상 원본 척도를 반환한다. 기본은 standardize=TRUE 이다.
더 많은 정보를 위해, help(glmnet) 또는 간단히 ?glmnet를 입력하라.
예와 같이, 우리는 α = 0.2로 설정했다 (ridge 회귀를 더 닮게), 그리고 관측의 나중의 반에 2배의 가중치를 준다. 여기서 너무 긴 표현을 피하기 위해, 우리는 nlambda를 20으로 설정했다. 연습에서, 어쨌든, λ의 값의 수는 100(기본) 또는 더 크게 추천한다. 대부분의 경우에, 그것은 추가 비용을 유발하지 않는다 왜냐하면 알고리즘에서 사용된 warm-starts 때문, 그리고 비선형 모델에서 더 나은 수렴 특징으로 이끈다.
fit = glmnet(x, y, alpha = 0.2, weights = c(rep(1,50),rep(2,50)), nlambda = 20)
우리는 그러면 glmnet 객체를 출력할 수 있다.
print(fit) ## ## Call: glmnet(x = x, y = y, weights = c(rep(1, 50), rep(2, 50)), alpha = 0.2, ## nlambda = 20) ## ## Df %Dev Lambda ## [1,] 0 0.0000 7.939000 ## [2,] 4 0.1789 4.889000 ## [3,] 7 0.4445 3.011000 ## [4,] 7 0.6567 1.854000 ## [5,] 8 0.7850 1.142000 ## [6,] 9 0.8539 0.703300 ## [7,] 10 0.8867 0.433100 ## [8,] 11 0.9025 0.266700 ## [9,] 14 0.9101 0.164300 ## [10,] 17 0.9138 0.101200 ## [11,] 17 0.9154 0.062300 ## [12,] 17 0.9160 0.038370 ## [13,] 19 0.9163 0.023630 ## [14,] 20 0.9164 0.014550 ## [15,] 20 0.9164 0.008962 ## [16,] 20 0.9165 0.005519 ## [17,] 20 0.9165 0.003399
이것은 객체 fit으로 생성된 호출을 표시한다 그리고 열 Df (0이 아닌 계수의 수), %dev (백분율 편차 설명) 그리고 람다 (해당 계수 값)을 가지는 3개의 열 행렬.
(digits 옵션은 출력에서 의미있는 숫자를 지정하기 위해 사용할 수 있음을 주목하라.)
여기서 λ의 실제 수는 여기서 호출에서 지정된 수 보다 작다. 이유는 알고리즘의 중지 기준에 있다. 기본 입력 설정에 따라, 계산은 멈춘다 만약 편차의 아래 방향에서 부분 변량가 10-5 보다 적거나 또는 설명 변량의 부분이 0.999에 도달한다면. 마지막 몇 라인에서, 우리는 변량의 부분이 변화가 없음을 보고 따라서 계산은 중지 기준을 만났을 때 긑난다. 우리는 그런 내부 변수를 변경할 수 있다. 자세한 사항은, Appendix 절을 보거나 또는 help(glmnet.control)를 입력하라.
우리는 이전 절에서처럼 적합된 객체를 그릴 수 있다. plot 함수에 더 많은 옵션이 있다.
사용자는 X-축이 무엇인지 결정할 수 있다. xvar는 3개의 측정을 허용한다: 계수의 ℓ1-norm을 위한 "norm"(기본), 로그-람다 값을 위한 "lambda" 그리고 %편차 설명을 위한 "dev".
사용자는 또한 단순히 label = TRUE 설정으로 변수 연속 번호를 가지는 곡선에 표기할 수 있다.
로그-람다 값과 표기된 각각의 곡선에 대해 "fit"을 플롯해라.
plot(fit, xvar = "lambda", label = TRUE)
지금 우리가 %deviance에 대해서 플롯할 대, 우리는 매우 다른 그림을 얻는다. 이것은 훈련 데이터에서 백분률 편차 설명(percent deviance explained) 이다. 우리가 여기서 보는 무엇은 경로의 끝을 향해서 이 값은 많이 변화하지 않는다, 그러나 계수는 약간 크게 된다. 이것은 우리가 적합의 이 문제 부분에 관심을 가지게 한다. 이것은 특히 다른 모델에서도 참이 될 것이다, 로지스틱 회귀에서도 같음.
plot(fit, xvar = "dev", label = TRUE)
우리는 계수를 추출하고 그리고 λ의 어떤 값에서 예측할 수 있다. 2개의 일반적으로 사용되는 옵션:
- s는 추출될 어떤 값에서 λ의 값을 지정한다.
- exact는 계수의 추출 값을 요구하는지 아닌지를 가리킨다, 그것은, 만약 exact = TRUE 이면, 그리고 원본 적합에 포함되지 않는 s의 값에서 만들어져야 한다. s의 이 값은 object$lambda를 가지고 합성된다, 그리고 모델은 예측이 생성되기 전에 재적합한다. 만약 exact=FALSE 이면(기본), 그러면 예측 함수는 적합 알고리즘에 사용되는 람다를 가지고 일치하지 않는 s의 값에서 예측을 만들기 위해 선형 보간(linear interpolation)을 사용한다.
단순한 예;
any(fit$lambda == 0.5) ## [1] FALSE coef.exact = coef(fit, s = 0.5, exact = TRUE) coef.apprx = coef(fit, s = 0.5, exact = FALSE) cbind2(coef.exact, coef.apprx) ## 21 x 2 sparse Matrix of class "dgCMatrix" ## 1 1 ## (Intercept) 0.19657091 0.199098747 ## V1 1.17495906 1.174650452 ## V2 . . ## V3 0.52934375 0.531934651 ## V4 . . ## V5 -0.76126245 -0.760959480 ## V6 0.46627405 0.468209413 ## V7 0.06148079 0.061926756 ## V8 0.38048793 0.380301491 ## V9 . . ## V10 . . ## V11 0.14213614 0.143260991 ## V12 . . ## V13 . . ## V14 -0.91090226 -0.911207368 ## V15 . . ## V16 . . ## V17 . . ## V18 . 0.009196628 ## V19 . . ## V20 -0.86099392 -0.863117051
왼쪽 열은 exact = TRUE 그리고 FALSE를 위한 오른쪽이다. 우리는 0.01이 열(sequence)에 있지 않음을 위에서 확인한다 그리고 그것은 그러므로 약간 다르다, 비록 크지 않지만. 선형 보간은(Linear interpolation) 대부분 충분한다 만약 특별한 요구가 없으면.
사용자는 적합된 객체로 예측을 만들 수 있다. coef에서 옵션 외에도, 주요 인수 newx 이다, x의 새로운 값의 행렬. type 옵션은 사용자에게 예측 타입을 선택하도록 허용한다: * “link”가 적합된 값을 제공한다
- "response"는 "gaussian" family를 위한 "link"와 같음.
- "coefficients"는 s의 값에서 계수를 계산한다.
- "nonzero"는 s의 각 값에서 0이 아닌 계수의 색인 목록을 반환한다.
예를 들어,
predict(fit, newx = x[1:5,], type = "response", s = 0.05) ## 1 ## [1,] -0.9802591 ## [2,] 2.2992453 ## [3,] 0.6010886 ## [4,] 2.3572668 ## [5,] 1.7520421
λ = 0.05에서 처음 5개의 관측의 적합 값을 제공한다, 만약 s의 다중 값이 제공된다면, 예측 행렬이 생성된다.
사용자가 K-fold 교차-검증을 정의할 수 있다. 모든 glmnet 변수 외에도, cv.glmnet은 nfolds (배 수), foldid (사용자-제공 folds), type.measure(교차-검증을 위해 사용되는 손실)를 포함하는 그것의 특별한 변수를 가진다: * “deviance”또는 “mse”은 제곱 손실을 사용한다
- "mae"는 평균 절대 에러를 사용한다
예를 들어,
cvfit = cv.glmnet(x, y, type.measure = "mse", nfolds = 20)
20-배 교차-검증을 행한다, 평균 제곱 에러(mean squared error)를 기준으로 (하지만 기본).
병렬 컴퓨팅은 또한 cv.glmnet에서 지원한다. 이 작업을 만들기 위해, 사용자는 미리 병렬을 등록해야 한다.
우리는 여기서 비교의 단순한 예를 제공한다.
require(doMC) ## Loading required package: doMC ## Loading required package: iterators ## Loading required package: parallel registerDoMC(cores=2) X = matrix(rnorm(1e4 * 200), 1e4, 200) Y = rnorm(1e4) system.time(cv.glmnet(X, Y)) ## user system elapsed ## 3.173 0.097 3.275 system.time(cv.glmnet(X, Y, parallel = TRUE)) ## user system elapsed ## 0.429 0.061 1.475
위 제안에서, 병렬 컴퓨팅은 거대한-크기 문제를 위해서 특히 계산 처리를 현저하게 빨라질 수 있다.
cv.glmnet 객체에 함수 coef 와 predict는 glmnet 객체의 그것들과 비슷하다, 2개의 특별한 문자열이 s(요구 λ의 값)로 또한 지원되는 점을 제외하면: * “lambda.1se”: MSE가 최조 MSE의 1 표준 에러(1 standard error) 안에 있는 가장 큰 λ.
- "lambda.min": 최소 MSE에 도달한 것에서 λ.
cvfit$lambda.min ## [1] 0.07569327 coef(cvfit, s = "lambda.min") ## 21 x 1 sparse Matrix of class "dgCMatrix" ## 1 ## (Intercept) 0.14867414 ## V1 1.33377821 ## V2 . ## V3 0.69787701 ## V4 . ## V5 -0.83726751 ## V6 0.54334327 ## V7 0.02668633 ## V8 0.33741131 ## V9 . ## V10 . ## V11 0.17105029 ## V12 . ## V13 . ## V14 -1.07552680 ## V15 . ## V16 . ## V17 . ## V18 . ## V19 . ## V20 -1.05278699 predict(cvfit, newx = x[1:5,], s = "lambda.min") ## 1 ## [1,] -1.3638848 ## [2,] 2.5713428 ## [3,] 0.5729785 ## [4,] 1.9881422 ## [5,] 1.5179882
사용자는 사용될 fold를 제어할 수 있다. 여기서 우리는 같은 fold를 사용한다 그래서 우리는 α를 위한 값을 선택할 수 있다.
foldid=sample(1:10,size=length(y),replace=TRUE) cv1=cv.glmnet(x,y,foldid=foldid,alpha=1) cv.5=cv.glmnet(x,y,foldid=foldid,alpha=.5) cv0=cv.glmnet(x,y,foldid=foldid,alpha=0)
그들 모두 같은 plot에 놓기 위한 내장된 plot 함수가 없다, 그래서 우리는 여기서 우리의 것을 보여준다.
par(mfrow=c(2,2))
plot(cv1);plot(cv.5);plot(cv0)
plot(log(cv1$lambda),cv1$cvm,pch=19,col="red",xlab="log(Lambda)",ylab=cv1$name)
points(log(cv.5$lambda),cv.5$cvm,pch=19,col="grey")
points(log(cv0$lambda),cv0$cvm,pch=19,col="blue")
legend("topleft",legend=c("alpha= 1","alpha= .5","alpha 0"),
pch=19,col=c("red","grey","blue"))
우리는 lasso (alpha=1)이 여기서 최고임을 확인한다, 우리는 또한 사용된 람다의 범위는 알파에 따라 다르다는 것을 확인한다.
계수의 상한 그리고 하한 경계 (Coefficient upper and lower bounds)
최근에 모델의 범위를 향상시키는 특징이 추가되었다. 우리가 자신의 모델 적합을 원한다고 가정하자, 그러나 -0.7 보다 크기 0.5 보다 작은 계수로 제한한다. 이것은 upper.limits 과 lower.limits 인수를 통해서 쉽게 달성된다:
tfit=glmnet(x,y,lower=-.7,upper=.5) plot(tfit)
이것들은 오히려 임의의 제한이다: 종종 우리는 양수인 계수를 원한다, 그래서 우리는 단지 lower.limit을 0 이 되게 설정할 수 있다. (주의, 하한은 0 보다 크지 않아야 한다, 그리고 상한은 0 보다 작지 않아야 한다.) 이런 경계는 벡터가 될 수 있다, 각 계수에서 다른 값을 가짐. 만약 scalar로 주어지면, 같은 수가 모두에서 반복된다.
Penalty factors
이 인수는 사용자가 각 계수에 분리된 패널티 요소를 적용하도록 허용한다. 그것의 기본은 각 변수에서 1 이다, 그러나 다른 값을 지정할 수 있다. 특히, penalty.factor가 0을 가지는 어떤 변수는 전혀 패널티가 없다! νj를 j번째 변수를 위한 패털티 요소로 나타내자, 패널티 조건은
이 된다.
패널티 조건은 nvars의 합으로 내부적으로 재척도 됨을 주의하라.
이것은 사람들이 사전 지식 또는 변수 전체에 선호를 가질 때 매우 유용하다. 많은 경우에, 어떤 변수는 너무 중요해서 그들을 언제나 유지되도록 원할 수 있을 것이다, 패널티 요소를 0에 상응하게 설정함으로써 이루어 질 수 있다:
p.fac = rep(1, 20) p.fac[c(5, 10, 15)] = 0 pfit = glmnet(x, y, penalty.factor = p.fac) plot(pfit, label = TRUE)
우리는 0 패널티 요소를 가지는 3개의 변수모델에서 항상 놓여 있는 딱지(label)에서 확인한다, 다른 것은 전형적인 경로를 따른다 그리고 0으로 축소된다.
약간의 다른 유용한 인수. exclude는 하나를 모델이 있지만 어떤 변수를 완전하게 차단한다. 물론, 하나는 이것을 x에서 제외된 단순한 하위집합일 수 있다, 그러나 때로 exclude는 더 유용하다, 그것은 계수의 전체 벡터를 반환하므로, 단지 제외된 것들을 0으로 설정함으로써. 또한 기본값 TRUE인 intercept 인수가 있다; 만약 FALSE이면 교차점(intercept)는 0이 되도록 강제한다.
사용자 정의 플롯 (Customizing plots)
때로는, 특히 변수의 수가 작을 때, 우리는 플롯에 변수 딱지(labels) 추가를 원한다. glmnet은 기본적으로 많은 데이터를 의도했다. 이것은 plot.glmnet에 지원되지 않는다. 어쨌든, 그것을 하는 것은 쉽다, 따라오는 작은 장난스런 예제를 보여준다.
우리는 먼저 약간의 데이터를 생성한다, 10개의 변수를 가짐, 그리고 상상력의 부족 그리고 우리가 그들에게 단순한 글자 이름을 주기 쉽게. 우리는 그러면 glmnet 모델을 적합한다. 그리고 표준 플롯을 생성한다.
set.seed(101)
x=matrix(rnorm(1000),100,10)
y=rnorm(100)
vn=paste("var",1:10)
fit=glmnet(x,y)
plot(fit)
우리는 변수 이름을 가지는 곡선을 표시하기를 원한다. 여기서 이것을 하는 단순한 방법은, R에서 axis 명령을 사용함(그리고 그것을 사용자 정의 방법에 약간의 연구). 우리는 경로의 끝에 계수의 위치를 얻을 필요가 있다, 그리고 우리는 par 명령어를 사용하야 약간의 공간을 만들 필요가 있다, 그래서 우리의 딱지(labels)는 안에 적합될 것이다. 이것은 당신의 딱지가 얼마나 긴지 아는 것이 요구된다, 그러나 여기서 그들은 모두 상당히 짧다.
par(mar=c(4.5,4.5,1,4))
plot(fit)
vnat=coef(fit)
vnat=vnat[-1,ncol(vnat)] # remove the intercept,
# and get the coefficients at the end of the path
axis(4, at=vnat,line=-.5,label=vn,las=1,tick=FALSE, cex.axis=0.5)
우리는 딱지(labels)의 겹침을 피하기 위해 여기서 어떤것도 행하지 않았다, 결국에는 그들은 서로 근접해 있다. 이것은 약간 더 많은 작업이 필요하다, 그러나 아마도 최고는 스스로 해야 한다, 아무튼.
다중 응답 가우시안 군 (Multiresponse Gaussian Family)
다중응답 가우시안 군은 glmnet에서 family = "mgaussian"을 사용하여 얻을 수 있다. 그것은 위의 단일-응답 경우와 매우 비슷하다. 이것은 많은 (correlated, 상관 관계) 응답이 있는 때 유용하다 - 소위 "다중-작업 학습(multi-task learning)" 문제. 여기서 공유 변수로 어느 변수 선택 될지 포함한다, 그러므로 변수를 선택할 때, 계수는 각 응답을 적합한다. 대부분 옵션은 같다, 그래서 우리는 여기서 단일 응답 모델과 차이에 집중한다.
분명히, 이름에서 알 수 있듯이, y는 벡터가 아니다, 그러나 이 절에서 수치 응답의 행렬. 람다의 각 값에서 계수는 또한 결과로는 행렬이다.
여기서 우리는 따라오는 문제를 해결한다:
여기서 βj는 p × K 계수 행렬 β의 j번째 행이다, 그리고 우리는 각 단일 계수에 absolute penalty(절대 패널티)를 단일 예측자 xj를 위한 각 계수 K-벡터 βj에 group-lasso penalty로 대체한다.
우리는 설명하기 위해 미리 생선된 데이터 집합을 사용한다.
data(MultiGaussianExample)
우리는 데이터를 적합한다, 객체 mfit가 반환됨.
mfit = glmnet(x, y, family = "mgaussian")
다중응답 가우시안에서, glmnet에서 옵션은 대부분 단일-응답 경우와 같다, alpha, weights, nlambda, standardize 처럼. 주목해야 할 예외는 standardize.response는 오로지 mgaussian 군에만 있다. 기본값은 FALSE 이다. 만약 standardize.response = TRUE 이면, 그것은 응답 변수를 표준화 한다.
계수를 시각화하기 위해, 우리는 plot 함수를 사용한다.
plot(mfit, xvar = "lambda", label = TRUE, type.coef = "2norm")
우리가 type.coef = "2norm"으로 설정함을 주목하라. 이 설정 하에서, 하나의 곡선은 변수 마다 그려진다, 값이 l2 노름과 같음. 기본 설정은 type.coef = "coef" 이다, 계수 플롯이 각 응답에서 생성됨(다중 그림).
일반적인 그래픽 변수 이외에 xvar 와 label 2개가 있다. 그들은 단일-응답 경우와 같다.
우리는 함수 coef를 사용하여 λ의 요구 값에서 계수를 추출할 수 있다 그리고 predict로 예측을 만들 수 있다. 사용은 비슷하다 그리고 우리는 단지 여기서 predict의 예제를 제공한다.
predict(mfit, newx = x[1:5,], s = c(0.1, 0.01)) ## , , 1 ## ## y1 y2 y3 y4 ## [1,] -4.7106263 -1.1634574 0.6027634 3.740989 ## [2,] 4.1301735 -3.0507968 -1.2122630 4.970141 ## [3,] 3.1595229 -0.5759621 0.2607981 2.053976 ## [4,] 0.6459242 2.1205605 -0.2252050 3.146286 ## [5,] -1.1791890 0.1056262 -7.3352965 3.248370 ## ## , , 2 ## ## y1 y2 y3 y4 ## [1,] -4.6415158 -1.2290282 0.6118289 3.779521 ## [2,] 4.4712843 -3.2529658 -1.2572583 5.266039 ## [3,] 3.4735228 -0.6929231 0.4684037 2.055574 ## [4,] 0.7353311 2.2965083 -0.2190297 2.989371 ## [5,] -1.2759930 0.2892536 -7.8259206 3.205211
예측 결과는 각 응답 변수에서 예측 행렬인 처음 2 차원 그리고 응답 변수를 가리키는 세번째를 가지는 3-차원 배열에 저장된다
cvmfit = cv.glmnet(x, y, family = "mgaussian")
우리는 결과 cv.glmnet 객체 "cvmfit"를 플롯한다.
plot(cvmfit)
λ의 선택 옵션 값을 명확하게 보기 위해, 입력해라
cvmfit$lambda.min ## [1] 0.04731812 cvmfit$lambda.1se ## [1] 0.1316655
이전에서, 첫번째 하나는 최소 평균 제곱 에러가 도달한 것에서 값이다 그리고 두번째는 평균 제곱 에러가 최소의 1 표준 에러 안에 있는 가장 좋은 정규화 모델이다.
cv.glmnet 객체에서 예측은 거의 glmnet 객체에서 같은 작업이다. 우리는 여기서 상세를 뺐다.
댓글 없음:
댓글 쓰기