제 14강: 비선형 및 평활 회귀 II

데이터과학 입문

원중호

서울대학교 통계학과

June 2024

시작하기 전에

다음의 패키지가 설치되어 있지 않으면 설치한다.

# install.packages("MASS")
# install.packages("tidyverse")
# install.packages("boot")
# install.packages("KernSmooth")
# install.packages("gam")
# install.packages("mars")
# install.packages("nnet")
library(MASS)
library(tidyverse)
library(boot)
library(KernSmooth)
library(gam)
library(mda)
library(nnet)

sessionInfo()
R version 4.4.0 (2024-04-24)
Platform: aarch64-apple-darwin20
Running under: macOS Sonoma 14.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Asia/Seoul
tzcode source: internal

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] nnet_7.3-19        mda_0.5-4          class_7.3-22       gam_1.22-3        
 [5] foreach_1.5.2      KernSmooth_2.23-22 boot_1.3-30        lubridate_1.9.3   
 [9] forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2       
[13] readr_2.1.5        tidyr_1.3.1        tibble_3.2.1       ggplot2_3.5.1     
[17] tidyverse_2.0.0    MASS_7.3-60.2     

loaded via a namespace (and not attached):
 [1] utf8_1.2.4        generics_0.1.3    stringi_1.8.3     hms_1.1.3        
 [5] digest_0.6.35     magrittr_2.0.3    evaluate_0.23     grid_4.4.0       
 [9] timechange_0.3.0  iterators_1.0.14  fastmap_1.1.1     jsonlite_1.8.8   
[13] fansi_1.0.6       scales_1.3.0      codetools_0.2-20  cli_3.6.2        
[17] rlang_1.1.3       munsell_0.5.1     withr_3.0.0       yaml_2.3.8       
[21] tools_4.4.0       tzdb_0.4.0        colorspace_2.1-0  vctrs_0.6.5      
[25] R6_2.5.1          lifecycle_1.0.4   pkgconfig_2.0.3   pillar_1.9.0     
[29] gtable_0.3.5      glue_1.7.0        xfun_0.43         tidyselect_1.2.1 
[33] rstudioapi_0.16.0 knitr_1.46        htmltools_0.5.8.1 rmarkdown_2.26   
[37] compiler_4.4.0

평활회귀 – 1변수 함수 적합

산점도 평활화

  • 목적: 연속형 반응변수 \(y \in \mathbb{R}\) vs 연속형 설명변수 \(x \in \mathbb{R}\)의 산점도를 평활화한 곡선 얻기
  • 곡선의 함수적 형태에 대해 거의 알지 못함
    • 1부와의 큰 차이!
  • 고전적 방법
    • 제한된 범위의 변환 (예: 2차 및 3차 항을 모형에 추가)
    • 변수의 범위를 분할하고 조각별 상수 또는 조각별 함수를 사용

GAGurine 자료

as_tibble(GAGurine)
# A tibble: 314 × 2
     Age   GAG
   <dbl> <dbl>
 1  0     23  
 2  0     23.8
 3  0     16.9
 4  0     18.6
 5  0.01  17.9
 6  0.01  25.9
 7  0.01  16.5
 8  0.01  26.3
 9  0.01  26.9
10  0.01  17.9
# ℹ 304 more rows
  • 0–17세 어린이 314명의 소변 내 화학 물질 GAG 농도에 대한 자료
  • 전진선택법은 6차 다항식을 제안
GAG.lm <- lm(GAG ~ Age + I(Age^2) + I(Age^3) + I(Age^4) + I(Age^5) + 
              I(Age^6) + I(Age^7) + I(Age^8), data=MASS::GAGurine)
anova(GAG.lm)
Analysis of Variance Table

Response: GAG
           Df  Sum Sq Mean Sq  F value    Pr(>F)    
Age         1 12589.8 12589.8 593.5762 < 2.2e-16 ***
I(Age^2)    1  3750.8  3750.8 176.8409 < 2.2e-16 ***
I(Age^3)    1  1491.5  1491.5  70.3202 1.891e-15 ***
I(Age^4)    1   449.3   449.3  21.1814 6.137e-06 ***
I(Age^5)    1   174.3   174.3   8.2156 0.0044421 ** 
I(Age^6)    1   286.0   286.0  13.4831 0.0002842 ***
I(Age^7)    1    57.2    57.2   2.6979 0.1015115    
I(Age^8)    1    44.9    44.9   2.1173 0.1466717    
Residuals 305  6469.1    21.2                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
GAG.lm2 <- lm(GAG ~ Age + I(Age^2) + I(Age^3) + I(Age^4) + I(Age^5) + 
              I(Age^6), data=GAGurine)
mypred <- data.frame(Age = seq(0, 17, len = 200))
mypred$GAG <- predict(GAG.lm2, newdata = mypred)
ggplot(GAGurine, aes(x = Age, y = GAG)) +
  geom_point(alpha=0.7) +  # Scatterplot of points
  geom_line(data = mypred, aes(x = Age, y = GAG), color = "blue") +  # Polynomial regression line
  labs(title = "Degree 6 polynomial")

스플라인

  • 수직선을 매듭(knot)이라 부르는 점의 순서집합 \(\{ z_i \}\)으로 분할
  • 구간 \([z_i, z_{i+1}]\)에서 3차 스플라인은 3차 다항식이며 연속인 1차 및 2차 도함수를 가짐
    • 각 매듭은 3개의 제약식을 만족
  • 매듭이 \(n\)개인 경우 3차 스플라인을 나타내려면 \(n+4\)개의 모수가 필요
    • \(4(n + 1)\)개 자유도에서 \(3n\) 연속 조건을 뺌
  • B 스플라인 (basis splines): 임의의 스플라인 함수는 B 스플라인의 선형결함으로 표현 가능
  • R 함수 splines::bs()는 B 스플라인 기저 행렬을 생성하여 선형 회귀 적합에 포함될 수 있음 (회귀 스플라인)

자연 3차 스플라인 (natural cubic splines)

  • 제한된 형태의 \(B-\)스플라인은 경계 매듭(boundary knots) 바깥에서 선형이 되도록 B 스플라인을 제약

    • 내부 매듭 (internal knots) \(n\)개, 경계 매듭 2개: \(n+2\)개 모수

  • R 함수 splines::ns()는 자유도 df를 지정하여 설명변수의 경험분포의 분위수로서 내부 매듭을, 최대최소값으로 경계 매듭을 선택할 수 있음

  • 주의: 스플라인 모형으로부터의 예측은 함수의 기저가 설명변수의 관측치에 따라 달라지므로 주의가 필요.

    • predict.lm()을 사용하면 새로운 기저함수 집합이 형성되고 적합된 계수가 잘못 적용됨

평활 스플라인 (smoothing splines)

  • \(n\)개의 순서쌍 \((x_i, y_i)\)에 대해 \[ \sum w_i\left[y_i-f\left(x_i\right)\right]^2+\lambda \int\left(f^{\prime \prime}(x)\right)^2 \mathrm{~d} x \] 을 최소화하는 두 번 미분가능한 함수 \(f\)를 적합

  • \(x_i\)에 매듭이 있는 3차 스플라인

    • 자료점을 보간하지 않음(\(\lambda > 0\))
    • 자유도: \(\operatorname{tr}(S_{\lambda})\), \(S_{\lambda}\)는 \(\widehat{y}=S_{\lambda} y\)인 \(n \times n\) 행렬 (Green and Silverman, 1994, pp. 37-–8, Hastie, Tibshirani, and Friedman, 2009, Ch. 5)

  • R 함수 smooth.spline()은 \(\lambda\) 또는 자유도를 입력으로 받거나 그렇지 않으면 교차 검증을 통해 자동으로 \(\lambda\)를 선택

attach(GAGurine)
GAG_ns5  <- lm(GAG ~ splines::ns(Age, df = 5))
GAG_ns10 <- lm(GAG ~ splines::ns(Age, df = 10))
GAG_ns20 <- lm(GAG ~ splines::ns(Age, df = 20))
GAG_sms  <- smooth.spline(Age, GAG)
xx <- seq(min(Age), max(Age), length.out = 200)
detach(GAGurine)
predicted05 <- predict(GAG_ns5,  newdata = data.frame(Age=xx))
predicted10 <- predict(GAG_ns10, newdata = data.frame(Age=xx))
predicted20 <- predict(GAG_ns20, newdata = data.frame(Age=xx))
predicted_sms <- predict(GAG_sms,  newdata = data.frame(Age=xx))
predicted_sms <- rename(as_tibble(predicted_sms), Age = x, GAG = y)
predicted_sms$df <- "Smoothing"
mypred <- data.frame(Age = xx, GAG = predicted05, df = "5")
mypred <- rbind(mypred, data.frame(Age = xx, GAG = predicted10, df = "10"))
mypred <- rbind(mypred, data.frame(Age = xx, GAG = predicted20, df = "20"))
mypred <- rbind(mypred, predicted_sms)
ggplot(GAGurine) + 
  geom_point(aes(x=Age, y=GAG), alpha=0.5, size=0.5) + 
  geom_line(data=mypred, aes(x=Age, y=GAG, color=df, linetype=df)) +
  labs(title = "Splines")

국소 회귀

  • 주어진 값 \(x\) 근방의 자료점에 선형 또는 다항식 회귀모형을 적합한 다음 \(x\)에서의 예측치를 평활화된 값으로 사용하여 작동하는 평활화 기법

lowess (locally weighted scatterplot smoothing, Cleveland, (1979), (1981))

  • \(x\) 주위에 창(window)를 설정, \(x\)에 가까운 점이 가장 큰 가중치를 갖도록 창 내부의 자료점에 가중치 적용
  • 강인(robust) 가중 회귀로 \(x\)에서 값을 예측
  • 인수 f는 창에 포함되는 자료의 비율. 기본값 f = 2/3.

loess (locally estimated scatterplot smoothing, Cleveland, Grosse and Shyu, 1992)

  • 유사한 방식으로 1차원, 2차원 또는 그 이상의 차원에서 작동하는 lowess 개념의 확장
attach(GAGurine)
mypred <- as.data.frame(loess.smooth(Age, GAG))
names(mypred) <- c("Age", "GAG")
detach(GAGurine)
ggplot(GAGurine, aes(x=Age, y=GAG)) + 
  geom_point(alpha=0.75) + 
  # default: predict(loess(y ~ x))
  geom_smooth(se=FALSE, color="magenta") + 
  geom_line(data=mypred, aes(x=Age, y=GAG), color="blue") +
  labs(title = "loess")

커널 평활화

\[ \widehat{y}_i=\sum_{j=1}^n y_i K\left(\frac{x_i-x_j}{b}\right) / \sum_{j=1}^n K\left(\frac{x_i-x_j}{b}\right) \] - \(b\): 대역폭, \(K\): 커널 함수 (예: Gaussian)

attach(GAGurine)
mypred1 <- as.data.frame(ksmooth(Age, GAG, "normal", bandwidth = 1))
names(mypred1) <- c("Age", "GAG")
mypred1$bandwidth <- "1"
mypred2 <- as.data.frame(ksmooth(Age, GAG, "normal", bandwidth = 5))
names(mypred2) <- c("Age", "GAG")
mypred2$bandwidth <- "5"
mypred <- rbind(mypred1, mypred2)
detach(GAGurine)
ggplot(GAGurine, aes(x=Age, y=GAG)) + 
  geom_point(alpha=0.75) + 
  geom_line(data=mypred, aes(x=Age, y=GAG, color=bandwidth)) +
  labs(title = "ksmooth")

locfit

  • 커널 회귀는 국소적으로 상수함수를 적합
  • 다항함수(특히 홀수 차수)의 국소적 적합이 이론적 이점이 있음 (Wand and Jones, 1995)
attach(GAGurine)
h <- KernSmooth::dpill(Age, GAG)  # bandwidth selection
mypred1 <- as.data.frame(KernSmooth::locpoly(Age, GAG, degree = 0, bandwidth = h))
names(mypred1) <- c("Age", "GAG")
mypred1$degree <- "0"
mypred2 <- as.data.frame(KernSmooth::locpoly(Age, GAG, degree = 1, bandwidth = h))
names(mypred2) <- c("Age", "GAG")
mypred2$degree <- "1"
mypred3 <- as.data.frame(KernSmooth::locpoly(Age, GAG, degree = 2, bandwidth = h))
names(mypred3) <- c("Age", "GAG")
mypred3$degree <- "2"
mypred <- rbind(mypred1, mypred2, mypred3)
detach(GAGurine)
ggplot(GAGurine, aes(x=Age, y=GAG)) + 
  geom_point(alpha=0.75) + 
  geom_line(data=mypred, aes(x=Age, y=GAG, color=degree)) +
  labs(title = "locpoly")

가법 모형

일반화 가법모형(generalized additive models, GAM)

  • 선형회귀 모형 \[ Y=\alpha+\sum_{j=1}^p \beta_j X_j+\epsilon \]

  • 가법 모형: \(\beta_j X_j \gets f_j(X_j)\), \(f_j: \mathbb{R}\to\mathbb{R}\): \[ Y=\alpha+\underset{j=1} { \stackrel{p}{\sum} } f_j (X_j)+\epsilon \]

    • 임의의 함수 \(f_j\)를 허용하는 것은 유용하지 않으므로 평활함수로 제약
    • 함수 \(f_j\)가 \(\alpha\)를 포함할 수 있음

역적합(backfitting)

  • 잔차 \[ Y - \alpha - \sum_{k \ne j} f_k (X_k) \] 에 대해 산점도 평활자를 이용, 평활함수 \(f_j\)르 적합
    • 모형의 선형 항들은 최소제곱법으로 적합
    • 수렴될 때까지 반복
  • gam::gam()
    • 모형 공식은 평활 스플라인과 loess를 나타내는 s(x), lo(x) 항을 허용

rock 자료

  • 석유를 함유한 암석 12개 각각의 4개 단면에 대한 측정값(area, peri, shape—암석 단면에 있는 기공의 원형도) 포함
  • 목표: 다른 세 가지 측정값을 통해 투과성 perm(유체 흐름의 특성)을 예측
  • 투과율의 스케일(6.3–1300)이 크게 변하므로 로그 척도 사용
rock.lm <-  lm(log(perm) ~ area + peri + shape, data = rock) 
rock.gam <- gam::gam(log(perm) ~ s(area) + s(peri) + s(shape), 
                     control = gam::gam.control(maxit = 50), 
                     data = rock)
summary(rock.gam)

Call: gam::gam(formula = log(perm) ~ s(area) + s(peri) + s(shape), 
    data = rock, control = gam::gam.control(maxit = 50))
Deviance Residuals:
    Min      1Q  Median      3Q     Max 
-1.6848 -0.4702  0.1235  0.5430  1.2955 

(Dispersion Parameter for gaussian family taken to be 0.7446)

    Null Deviance: 126.9322 on 47 degrees of freedom
Residual Deviance: 26.0597 on 35.0005 degrees of freedom
AIC: 134.8981 

Number of Local Scoring Iterations: NA 

Anova for Parametric Effects
              Df Sum Sq Mean Sq  F value    Pr(>F)    
s(area)    1.000 14.577  14.577  19.5788 8.988e-05 ***
s(peri)    1.000 75.670  75.670 101.6321 6.856e-12 ***
s(shape)   1.000  0.413   0.413   0.5548    0.4613    
Residuals 35.001 26.060   0.745                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova for Nonparametric Effects
            Npar Df  Npar F  Pr(F)
(Intercept)                       
s(area)           3 0.34165 0.7953
s(peri)           3 0.94055 0.4315
s(shape)          3 1.43200 0.2500
anova(rock.lm,  rock.gam) 
Analysis of Variance Table

Model 1: log(perm) ~ area + peri + shape
Model 2: log(perm) ~ s(area) + s(peri) + s(shape)
  Res.Df    RSS     Df Sum of Sq      F Pr(>F)
1 44.000 31.949                               
2 35.001 26.060 8.9995    5.8891 0.8789 0.5528
rock.gam1 <- gam::gam(log(perm) ~ area + peri + s(shape), data = rock) 
anova(rock.lm,  rock.gam1,  rock.gam)
Analysis of Variance Table

Model 1: log(perm) ~ area + peri + shape
Model 2: log(perm) ~ area + peri + s(shape)
Model 3: log(perm) ~ s(area) + s(peri) + s(shape)
  Res.Df    RSS     Df Sum of Sq      F Pr(>F)
1 44.000 31.949                               
2 41.000 28.999 3.0000    2.9495 1.3205 0.2833
3 35.001 26.060 5.9995    2.9396 0.6581 0.6835
  • 결과는 각 일변수함수가 하나의 선형 자유도와 세 개의 비선형 자유도를 가짐을 보여준다.
  • RSS를 31.95(선형 적합)에서 26.06(완전비선형적합)으로 줄이는 것은 추가 자유도 9에서는 유의하지 않음
  • peri에 대한 비선형 항은 유의하지 않아 보임
    • shape의 비선형 항만으로 RSS는 29.00
par(mfrow =  c(2, 3), pty = "s") 
plot(rock.gam,  se = T)
plot(rock.gam1,  se = T)
  • 평활함수 중 하나만 비선형으로 나타나고 이조차 완전히 설득력이 있는 것도 아니다.

mgcv::gam()

  • gam::gam() 사용시 가장 큰 어려움 중 하나는 각 항의 평활화 정도를 정하는 것
    • s()의 경우 기본 자유도는 4, lo()의 경우 f=0.5
  • mgcv::gam()에서는 평활화 정도가 적합의 일부로 추정(Wood, Pya, & Saefken, 2016)
mgcv::gam(log(perm) ~ s(area) + s(peri) + s(shape), data=rock)

Family: gaussian 
Link function: identity 

Formula:
log(perm) ~ s(area) + s(peri) + s(shape)

Estimated degrees of freedom:
1.0 1.0 1.4  total = 4.4 

GCV score: 0.7886465

다중 가법 회귀 스플라인(MARS, Friedman, 1991)

  • 기본값으로 선형 스플라인(조각 연속 1차 함수)을 각 변수에 적합하는 가법 모형
    • 고차 교호작용 허용
  • 조각의 수는 소프트웨어적으로 선택

cpus 자료

head(MASS::cpus)
            name syct mmin  mmax cach chmin chmax perf estperf
1  ADVISOR 32/60  125  256  6000  256    16   128  198     199
2  AMDAHL 470V/7   29 8000 32000   32     8    32  269     253
3  AMDAHL 470/7A   29 8000 32000   32     8    32  220     253
4 AMDAHL 470V/7B   29 8000 32000   32     8    32  172     253
5 AMDAHL 470V/7C   29 8000 16000   32     8    16  132     132
6  AMDAHL 470V/8   26 8000 32000   64     8    32  318     290
  • 미니컴퓨터와 메인프레임의 혼합 벤치마크 성능에 대한 연구(Ein-Dor and Feldmesser, 1987)
  • syst: 사이클 시간(나노초), cach: 캐시 크기(KB), mmin, mmax: 메인 메모리 크기(KB), chmin, chmax: 채널 수
    • 측정값은 IBM 370/158-3을 기준으로 정규화
  • 선형모형(교재 6.9절)
set.seed(123)
cpus2 <- select(cpus, -c(name, estperf)) # excludes names, authors' predictions
cpus.samp <- sample(seq_len(nrow(cpus2)), 100)
cpus.lm <- lm(log10(perf)~., data=cpus2 |> 
                              dplyr::slice(cpus.samp) )
cpus.lm2 <- stepAIC(cpus.lm, trace=FALSE)
test.cpus <- function(fit)
  sqrt(sum((log10(cpus2 |> slice(-cpus.samp) |> select(perf)) -
            predict(fit, cpus2 |> slice(-cpus.samp)))^2)/109)
test.cpus(cpus.lm)
[1] 0.2270138
cpus.lm2 <- stepAIC(cpus.lm, trace=FALSE)
cpus.lm2$anova
Stepwise Model Path 
Analysis of Deviance Table

Initial Model:
log10(perf) ~ syct + mmin + mmax + cach + chmin + chmax

Final Model:
log10(perf) ~ syct + mmin + mmax + cach

     Step Df    Deviance Resid. Df Resid. Dev       AIC
1                               93   2.936109 -338.8085
2 - chmax  1 0.001168625        94   2.937277 -340.7687
3 - chmin  1 0.042380328        95   2.979658 -341.3362
test.cpus(cpus.lm2)
[1] 0.2273272
  • 따라서 더 작은 모델을 선택한다고 해서 예측성능이 향상되지는 않음
  • 선형스플라인 가법모형 (역적합과 다른 방식으로 적합):
Xin <- as.matrix(cpus2 |> dplyr::slice(cpus.samp) |> select(-perf) )
cpus.mars <- mda::mars(Xin, log10(cpus2 |> slice(cpus.samp) |> select(perf) )) 
test2 <- function(fit) {  # test error
    Xp <- as.matrix(cpus2 |> slice(-cpus.samp) |> select(-perf) ) #cpus2[-cpus.samp,  1:6])
    sqrt(sum((log10(cpus2 |> slice(-cpus.samp) |> select(perf) ) - predict(fit, Xp))^2)/109)
}
test2(cpus.mars) 
[1] 0.2064594
showcuts <- function(obj) {
    tmp  <-  obj$cuts[obj$sel, ]
    dimnames(tmp)  <- list(NULL,  colnames(Xin)) 
    tmp
}
showcuts(cpus.mars) 
     syct mmin  mmax cach chmin chmax
[1,]    0    0     0    0     0     0
[2,]    0    0 20970    0     0     0
[3,]    0    0     0   30     0     0
[4,]    0    0     0   30     0     0
[5,]    0 1500     0    0     0     0
[6,]    0    0     0    0     0    38
[7,]    0    0     0    0     0    38
[8,]  124    0     0    0     0     0
# examine the fitted functions
Xp <- matrix(cpus2 |> slice(cpus.samp) |> select(-perf) |> map_dbl(mean), 
             100, 6, byrow = T)
xr <- sapply(cpus2, range) # variable ranges
df_plt <- list()
for (i in seq_len(ncol(Xp))) {
    Xp1 <- Xp  # hold other variables
    xvals <- seq(xr[1, i], xr[2, i], len = 100)
    Xp1[, i] <- xvals
    df_plt[[i]] <- data.frame(x = xvals, y = predict(cpus.mars, Xp1), 
                              variable = names(cpus2)[i])
}
df_long <- dplyr::bind_rows(df_plt)
ggplot(df_long, aes(x = x, y = y, color = variable)) + geom_line() +
  facet_wrap(~ variable, ncol = 3, scales = "free")
## allow pairwise interaction terms
cpus.mars2 <- mars(Xin, log10(cpus2 |> slice(cpus.samp) |> select(perf) ), degree = 2) 
showcuts(cpus.mars2) 
      syct mmin  mmax cach chmin chmax
 [1,]    0    0     0    0     0     0
 [2,]    0    0 20970    0     0     0
 [3,]    0    0 20970    0     0     0
 [4,]    0    0 20970    6     0     0
 [5,]    0    0 20970    6     0     0
 [6,]    0 2620     0    0     0     0
 [7,]    0    0     0    0     5     0
 [8,]  185 2620     0    0     0     0
 [9,]    0    0 20970    0     0    38
test2(cpus.mars2) 
[1] 0.231812
## allow arbitrary interactions
cpus.mars6 <- mars(Xin, log10(cpus2 |> slice(cpus.samp) |> select(perf) ), degree = 6) 
showcuts(cpus.mars6) 
      syct mmin  mmax cach chmin chmax
 [1,]    0    0     0    0     0     0
 [2,]    0    0 20970    0     0     0
 [3,]    0    0 20970    0     0     0
 [4,]    0    0 20970    6     0     0
 [5,]    0 2620     0    0     0     0
 [6,]    0    0     0    0     5     0
 [7,]  185 2620     0    0     0     0
 [8,]  200    0 20970    6     0     0
 [9,]  160    0 20970    6     0     0
[10,]  185 2620  6300    0     0     0
test2(cpus.mars6) 
[1] 0.2198926
  • 조각별 교호작용 항(degree=2)이나 임의의 교호작용(degree=6)을 허용하는 것은 예측의 효율성에 거의 차이가 없다.

사영 추적 회귀

사영 추적 회귀모형 (Friedman and Stuetzle, 1981)

\[ Y = \alpha_0 + \sum_{j=1}^{M} f_j (\boldsymbol{\alpha}_j^T \boldsymbol{X}) + \epsilon , \quad \|\boldsymbol{\alpha}_j\| = 1 \]

  • 1차원으로 사영된 변수들에 가법모형 적용
  • 차원 \(M\) 및 방향 \(\boldsymbol{\alpha}_j\)는 사용자가 선택

  • 충분히 큰 \(M\)에 대해 \(\boldsymbol{X}\)에 대한 임의의 연속 함수를 근사할 수 있다(uniform approximation on a compact set, Diaconis and Shahshahani, 1984).

  • 가법항 \(f_j (\boldsymbol{\alpha}_j^T \boldsymbol{X})\)은 한 방향을 제외한 모든 방향에서 일정하기 때문에 능형 함수(ridge function)로 불림

stats::ppr()

  • 최소제곱법으로 사영 추적 회귀모형 적합
    1. \(M_{\max}\)개 항을 순차적으로 적합
    2. 각 단계에서 가장 효과가 작은 항 제거
    3. 항이 \(M\) 개가 될 때까지 재적합
  • \(M,..., M_{max}\) 항 모형들의 RSS 반환
  • \(f_j\): Friedman’s supersmoother (supsmu), 평활스플라인

rock 자료

rock1 <- rock |> dplyr::mutate(area = area / 10000, peri = peri / 10000)
(rock.ppr <- ppr(log(perm) ~ area + peri + shape, 
                 data = rock1, nterms = 2, max.terms = 5))
Call:
ppr(formula = log(perm) ~ area + peri + shape, data = rock1, 
    nterms = 2, max.terms = 5)

Goodness of fit:
 2 terms  3 terms  4 terms  5 terms 
8.737806 5.289517 4.745799 4.490378 
summary(rock.ppr) 
Call:
ppr(formula = log(perm) ~ area + peri + shape, data = rock1, 
    nterms = 2, max.terms = 5)

Goodness of fit:
 2 terms  3 terms  4 terms  5 terms 
8.737806 5.289517 4.745799 4.490378 

Projection direction vectors ('alpha'):
      term 1      term 2     
area   0.34357179  0.37071027
peri  -0.93781471 -0.61923542
shape  0.04961846  0.69218595

Coefficients of ridge terms ('beta'):
   term 1    term 2 
1.6079271 0.5460971
  • \(\beta_{ij}\): 다중 응답 회귀 모형에 대응 \[ Y_i = \alpha_{i0} + \sum_{j=1}^{M} \beta_{ij} f_j (\boldsymbol{\alpha}_j^T \boldsymbol{X}) + \epsilon \]
    • 능형함수 \(f_j\)는 자료의 사영에 대해 평균 0 및 분산 1이 되도록 스케일링
par(mfrow = c(3,  2))
plot(rock.ppr)
plot(update(rock.ppr,  bass = 5)) # supsmu control
# use smoothing spline with GCV
plot(rock.ppr2 <- update(rock.ppr,  sm.method = "gcvspline", gcvpen = 2))
summary(rock.ppr2) 
Call:
ppr(formula = log(perm) ~ area + peri + shape, data = rock1, 
    nterms = 2, max.terms = 5, sm.method = "gcvspline", gcvpen = 2)

Goodness of fit:
 2 terms  3 terms  4 terms  5 terms 
22.52331 21.56513 21.56516  0.00000 

Projection direction vectors ('alpha'):
      term 1      term 2     
area   0.34852171  0.44220823
peri  -0.93698661 -0.85610648
shape  0.02426324 -0.26745763

Coefficients of ridge terms ('beta'):
   term 1    term 2 
1.4647599 0.2008702 

Equivalent df for ridge terms:
term 1 term 2 
  2.73   2.00
  • 2항 모형이 효과적인 것으로 보이고, area, peri가 지배적

능형함수 시각화

summary(rock1) # to find the ranges of the variables
      area             peri             shape              perm        
 Min.   :0.1016   Min.   :0.03086   Min.   :0.09033   Min.   :   6.30  
 1st Qu.:0.5305   1st Qu.:0.14149   1st Qu.:0.16226   1st Qu.:  76.45  
 Median :0.7487   Median :0.25362   Median :0.19886   Median : 130.50  
 Mean   :0.7188   Mean   :0.26822   Mean   :0.21811   Mean   : 415.45  
 3rd Qu.:0.8870   3rd Qu.:0.39895   3rd Qu.:0.26267   3rd Qu.: 777.50  
 Max.   :1.2212   Max.   :0.48642   Max.   :0.46413   Max.   :1300.00  
Xp  <- expand.grid(area = seq(0.1, 1.2, 0.05), peri = seq(0, 0.5, 0.02), 
                   shape = 0.2)
rock.grid <- cbind(Xp, fit = predict(rock.ppr2, Xp)) 
lattice::wireframe(fit ~ area +  peri, rock.grid,
                   screen = list(z = 160,  x = -60),  
                   aspect = c(1,  0.5), drape = T)

예제: cpus 자료

supsmu

(cpus.ppr <- ppr(log10(perf)~., data = cpus2 |> dplyr::slice(cpus.samp),
                 nterms = 2, max.terms = 10, bass = 5)) 
Call:
ppr(formula = log10(perf) ~ ., data = dplyr::slice(cpus2, cpus.samp), 
    nterms = 2, max.terms = 10, bass = 5)

Goodness of fit:
 2 terms  3 terms  4 terms  5 terms  6 terms  7 terms  8 terms  9 terms 
2.000989 1.729778 1.726409 1.416302 1.295887 1.075994 1.111010 1.040365 
10 terms 
1.178248 
cpus.ppr <- ppr(log10(perf)~., data = cpus2 |> dplyr::slice(cpus.samp),
                nterms = 8, max.terms = 10, bass = 5)  # select 8 terms
test.cpus(cpus.ppr) 
[1] 0.2197131

spline

ppr(log10(perf)~., data = cpus2 |> dplyr::slice(cpus.samp),
    nterms = 2, max.terms = 10, sm.method = "spline")
Call:
ppr(formula = log10(perf) ~ ., data = dplyr::slice(cpus2, cpus.samp), 
    nterms = 2, max.terms = 10, sm.method = "spline")

Goodness of fit:
 2 terms  3 terms  4 terms  5 terms  6 terms  7 terms  8 terms  9 terms 
1.913118 1.863657 1.731443 1.463755 1.311034 1.261891 1.219237 1.207179 
10 terms 
1.161124 
cpus.ppr2 <- ppr(log10(perf)~., data = cpus2 |> dplyr::slice(cpus.samp),
                 nterms = 7,  max.terms = 10, sm.method = "spline")
test.cpus(cpus.ppr2) 
[1] 0.2349891
  • 예측값 비교
res3 <- unlist(log10(cpus2 |> dplyr::slice(-cpus.samp) |> select(perf)) - 
                  predict(cpus.ppr, cpus2 |> dplyr::slice(-cpus.samp)) ) 
res2 <- unlist(log10(cpus2 |> slice(-cpus.samp) |> select(perf)) -
                  predict(cpus.lm2, cpus2 |> slice(-cpus.samp)))  # lm
wilcox.test(res2^2, res3^2, paired = T, alternative = "greater") 

    Wilcoxon signed rank test with continuity correction

data:  res2^2 and res3^2
V = 3313, p-value = 0.1705
alternative hypothesis: true location shift is greater than 0

신경망

인공신경망

  • 순방향(feed-forward) 신경망은 선형 회귀 함수를 유연하게 일반화 (Bishop, 1995; Ripley, 1996; Goodfellow, Bengio, and Courville, 2016).

  • 모수가 너무 많은(overparameterized) 비선형 회귀 모형

  • 단일 은닉층 순방향 신경망: \[ y_k = \phi_o \left( \alpha_k+\sum_h w_{hk}\phi_h \left( \alpha_h+\sum_i w_{ih}x_i \right) \right) \]

    • \(\phi_h\): 은닉층 단위의 활성화 함수. 로지스틱 함수 \(\ell(z) = \frac{e^z}{1+e^z}\) 또는 ReLU \(\ell(z) = (z)_{+}\).
    • \(\phi_o\): 출력 단위의 활성화 함수. 선형, 로지스틱 또는 임계값 함수(\(\phi_o(x)=I(x>0)\)).
      • 선형 출력함수의 경우 사영 추적 회귀의 특수한 경우

  • 층 건너뛰기(skip-layer, skip-connection) \[ y_k = \phi_o \left( \alpha_k+\sum_{i \rightarrow k} w_{ik}x_i + \sum_{j \rightarrow k} w_{jk}\phi_h \left( \alpha_j+\sum_{i \rightarrow j} w_{ij}x_i \right) \right) \]

  • 모수: \(\{\alpha_j\}\) (biases), \(\{w_{jk}\}\) (weights)

  • 선형 출력함수의 경우 은닉층의 크기를 늘리면 옹골집합에서 모든 연속 함수를 균일하게 근사 (Cybenko, 1989; Hornik, Stinchcombe and White, 1989)

  • 근사 결과는 비건설적: 실제로는 적합 기준 (오차제곱합, 소프트맥스 등)을 최소화하도록 가중치를 선택 \[ \small E=\underset{p}{\sum} || t^p - y^p ||^2 , \quad E=\sum_{p}\sum_{k}-t_k^p \log p_k^p, \quad p_k^p=\frac{e^{y_k^p}}{\sum_{c=1}^K e^{y_c^p}} \]

    • \(t^p\)는 타겟, \(y^p\)는 \(p\)번째 예제 패턴의 출력값
  • 조절(regularization): \(E+\lambda C(f)\), \(f\) 회귀함수
    • 모형의 자유도 제약, 과적합 방지
    • Weight decay: \(C(f) = \sum_{i,j} w_{ij}^2\) (입력이 \([0, 1]\) 범위로 조정되는 경우)
    • 참고: 드롭아웃

nnet::nnet()

nnet(formula, data, weights, size, Wts, linout = F, entropy = F, 
    softmax = F, skip = F, rang = 0.7, decay = 0, maxit = 100,
    trace = T)
  • 인수
    • size: 은닉층 단위의 개수
    • Wts: \(w_{ij}\) 초기값
    • linout: 선형 출력함수?
    • entropy: 최소제곱 적합이 아닌 엔트로피 적합?
    • softmax: 소프트맥스 적합 (분류 문제)?
    • skip: 건너뛰기 허용?
    • rang: Wts가 누락된 경우 runif(n,-rang, rang)의 임의의 가중치를 사용한다.
    • decay: 가중치붕괴 모수 \(\lambda\)
    • maxit: 최적화에 대한 최대 반복 횟수
    • Hess: 해의 헤시안 행렬이 반환되어야 하는가?
    • trace: 적합과정 출력?

rock 자료

rock1 <- rock |> dplyr::mutate(area = area / 10000, peri = peri / 10000)
rock.nn <- nnet::nnet(log(perm) ~ area + peri + shape, data = rock1, 
                      size = 3, decay = 1e-3, linout = T,  skip = T,
                      maxit = 1000, Hess = T)
# weights:  19
initial  value 1519.181470 
iter  10 value 32.184361
iter  20 value 21.786682
iter  30 value 16.573606
iter  40 value 14.773713
iter  50 value 14.653865
iter  60 value 14.532256
iter  70 value 14.493436
iter  80 value 14.486145
iter  90 value 14.483477
iter 100 value 14.454209
iter 110 value 14.316777
iter 120 value 13.904518
iter 130 value 13.507059
iter 140 value 13.488607
iter 150 value 13.469334
iter 160 value 13.432777
iter 170 value 13.267998
iter 180 value 13.015840
iter 190 value 12.861436
iter 200 value 12.791917
iter 210 value 12.787941
iter 220 value 12.786843
iter 230 value 12.786574
iter 240 value 12.786264
iter 250 value 12.786038
final  value 12.785771 
converged
summary(rock.nn)
a 3-3-1 network with 19 weights
options were - skip-layer connections  linear output units  decay=0.001
 b->h1 i1->h1 i2->h1 i3->h1 
 -8.94  15.47 -20.84  22.05 
 b->h2 i1->h2 i2->h2 i3->h2 
 -4.29  11.28 -16.38   8.30 
 b->h3 i1->h3 i2->h3 i3->h3 
 -2.12  11.24  -1.60   4.12 
  b->o  h1->o  h2->o  h3->o  i1->o  i2->o  i3->o 
  8.18  -8.36  15.28  -5.78  -2.61  -3.12   0.98 
sum((log(rock$perm) - predict(rock.nn))^2)
[1] 10.49205
eigen(rock.nn$Hessian, T)$values
 [1] 1.351179e+03 1.137257e+02 5.296300e+01 1.775020e+01 5.015666e+00
 [6] 3.973783e+00 1.377261e+00 9.390317e-01 4.443876e-01 2.564340e-01
[11] 1.299196e-01 4.984606e-02 2.623271e-02 1.275362e-02 8.286848e-03
[16] 5.596481e-03 5.304720e-03 4.482875e-03 1.692169e-04

예제 : 자료 cpus

  • 신경망 모형을 효과적으로 사용하려면 변수를 스케일하는 것이 필수적
attach(cpus2)
cpus3 <- data.frame(syct = syct - 2, mmin = mmin - 3, mmax = mmax - 4,
                    cach = cach / 256,  chmin = chmin / 100,
                    chmax = chmax / 100, perf = perf)
detach(cpus2)
test.cpus <-  function(fit) {
    sqrt(sum((log10(cpus3 |> slice(-cpus.samp) |> select(perf)) - 
              predict(fit, cpus3 |> slice(-cpus.samp)))^2)/109)
}
cpus.nn1 <- nnet(log10(perf)~., cpus3 |> slice(cpus.samp), linout = T, 
                 skip = T,  size = 0) # linear model
# weights:  7
initial  value 942410632.692998 
iter  10 value 3.776454
final  value 2.936109 
converged
test.cpus(cpus.nn1)
[1] 0.2270138
  • 비선형성(은닉층) 도입
cpus.nn2 <- nnet(log10(perf)~., cpus3 |> slice(cpus.samp), linout = T, 
                 skip=T, size=4, decay=0.01, maxit = 1000, trace=F)
test.cpus(cpus.nn2) 
[1] 0.2006723
cpus.nn3 <- nnet(log10(perf)~., cpus3 |> slice(cpus.samp), linout = T,
                 skip=T, size=10, decay=0.01, maxit = 1000, trace=F)
test.cpus(cpus.nn3)
[1] 0.2268575
cpus.nn4 <- nnet(log10(perf)~., cpus3 |> slice(cpus.samp), linout = T, 
                 skip=T, size=25, decay=0.01, maxit = 1000, trace=F)
test.cpus(cpus.nn4)  
[1] 0.3481916
  • 적합도가 은닉층 단위의 수(충분한 경우)보다는 가중치 붕괴 정도에 의해 거의 완전히 제어된다는 것을 보여줌

교차검증에 의한 가중치 붕괴 선택

CVnn.cpus <- function(formula, data = cpus3,
    size = c(0, 4, 4, 10, 10),
    lambda = c(0, rep(c(0.003, 0.01), 2)),
    nreps = 5, nifold = 10, ...) {
  CVnn1 <- function(formula, data, nreps=1, ri,  ...) {
    truth <- log10(data$perf)
    res <- numeric(length(truth))
    cat("  fold")
    for (i in sort(unique(ri))) {
      cat(" ", i,  sep="")
      for(rep in 1:nreps) {
        learn <- nnet(formula, data[ri !=i,], trace=FALSE, ...)
        res[ri == i] <- res[ri == i] +
                        predict(learn, data[ri == i,])
      }
    }
    cat("\n")
    sum((truth - res/nreps)^2)
  }
  choice <- numeric(length(lambda))
  ri <- sample(nifold, nrow(data), replace = TRUE)
  for(j in seq(along=lambda)) {
    cat("  size =", size[j], "decay =", lambda[j], "\n")
    choice[j] <- CVnn1(formula, data, nreps=nreps, ri=ri,
                       size=size[j], decay=lambda[j], ...)
    }
  cbind(size=size, decay=lambda, fit=sqrt(choice/100))
}
CVnn.cpus(log10(perf) ~ ., data = cpus3 |> slice(cpus.samp),
          linout = TRUE, skip = TRUE, maxit = 1000)
  size = 0 decay = 0 
  fold 1 2 3 4 5 6 7 8 9 10
  size = 4 decay = 0.003 
  fold 1 2 3 4 5 6 7 8 9 10
  size = 4 decay = 0.01 
  fold 1 2 3 4 5 6 7 8 9 10
  size = 10 decay = 0.003 
  fold 1 2 3 4 5 6 7 8 9 10
  size = 10 decay = 0.01 
  fold 1 2 3 4 5 6 7 8 9 10
     size decay       fit
[1,]    0 0.000 0.1902077
[2,]    4 0.003 0.1757065
[3,]    4 0.010 0.1913253
[4,]   10 0.003 0.1926077
[5,]   10 0.010 0.1999800
  • 교차 검증 결과는 모형 선택에 다소 둔감
  • 비선형성은 타당하지 않은 것으로 보임

결론

  • 선형모형의 확장

    1. 블랙박스: 사영 추적 회귀, 신경망
    2. 사용자 제어권: (모수적) 비선형 회귀, 가법 모형

  • 후자가 해석상 이득을 얻을 수는 있지만 rock 자료의 예에서 보듯 충분히 일반적이지 않을 수 있음

  • 해석이 반드시 쉽지도 않음 (MASS Fig. 8.12)

  • 어떤 모형이 특정 문제에 가장 적합할지는 문제의 목적, 특히 예측이 중요한지 해석이 중요한지에 따라 다름

  • 모형을 뒷받침하는 분포 이론이 거의 없음 — 과적합, 과설명의 위험

    • 최근 연구 결과: conformal prediction, 예: jackknife+