제 13강: 비선형 및 평활 회귀 I

데이터과학 입문

Author
Affiliation

원중호

서울대학교 통계학과

Published

June 2024

시작하기 전에

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

# install.packages("MASS")
# install.packages("tidyverse")
# install.packages("boot")
library(MASS)
library(tidyverse)
library(boot)

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] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] boot_1.3-30     lubridate_1.9.3 forcats_1.0.0   stringr_1.5.1  
 [5] dplyr_1.1.4     purrr_1.0.2     readr_2.1.5     tidyr_1.3.1    
 [9] tibble_3.2.1    ggplot2_3.5.1   tidyverse_2.0.0 MASS_7.3-60.2  

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

비선형 회귀

  • 선형 회귀에서 평균면은 표본 공간의 평면, 비선형 회귀에서는 임의의 곡면
  • 다른 모든 측면에서는 모형이 동일 (예: 최소제곱법)
  • 실제 대부분의 비선형 회귀 모형의 평균면은 가능도가 가장 높은 영역에서 평면으로 근사 가능
    • 선형회귀를 기반으로 한 몇 가지 좋은 근사치를 사용 가능
    • 여전히 까다로운 계산 및 추론 문제가 생길 수 있다.
  • Bates and Watts (1998)

평활 회귀

예시: wtloss 자료

wtloss 자료

  • 체중 감량 프로그램을 수행 중인 비만 환자는 지방 조직이 감쇄하는 비율로 감소하는 경향이 있음
  • 변수
    • Days: 프로그램 시작 이후의 시간(일)
    • Weight: 표준 조건에서 측정된 환자의 체중(kg)
  • 48세, 키 193cm (6’ 4”), 체격이 큰 남성 환자에 관한 자료
ggplot(MASS::wtloss, aes(x=Days, y=Weight)) + geom_point() + 
      scale_y_continuous(
        "Weight (kg)", 
        sec.axis = sec_axis(~ . * 2.205, name = "Weight (lb)")
      )
Figure 1: Weight loss for an obese patient.

비선형 회귀

  • 다항식 회귀 모형은 관찰된 범위 내에서 이러한 자료를 매우 잘 설명할 수 있지만 이 범위 밖에서는 크게 실패할 수 있다.

  • 일부 이론적 및 경험적 근거를 가진 보다 유용한 모형: \[ y= \beta_0+\beta_1 2^{-t/\theta}+\epsilon \]

  • 모수에 대해 비선형

    • \(\beta_0\) : 궁극적인 순수 체중 또는 점근선 (선형모수)
    • \(\beta_1\) : 감량할 체중 (선형모수)
    • \(\theta\) : 목표치 절반을 감량하는 데 걸리는 시간(비선형모수)

  • 선형모수: 2계 편도함수가 0

비선형 회귀 모형 적합

비선형 회귀 모형의 일반적 형태

\[ y=\eta(\mathbf{x}, \boldsymbol{\beta}) + \epsilon \]

  • \(\mathbf{x} \in \mathbb{R}^p\) 공변량 벡터
  • \(\boldsymbol{\beta} \in \mathbb{R}^p\) : 회귀 모수 벡터
  • \(\epsilon \sim N(0, \sigma^2)\): 오차
  • wtloss 예제에서 \(\boldsymbol{\beta} = (\beta_0, \beta_1, \theta)^T\)

  • \(\boldsymbol{y}=(y_1, \dotsc, y_n)^\top\)가 표본 반응변수 벡터이면 비선형회귀모형 하에서 \(\boldsymbol{\eta} (\boldsymbol{\beta})=(\eta(\boldsymbol{x}_1, \boldsymbol{\beta}), \dotsc, \eta(\boldsymbol{x}_n, \boldsymbol{\beta}))^\top\)는 \(\boldsymbol{y}\)는 평균 벡터

  • \(\boldsymbol{\beta}\)의 최대가능도추정치는 최소제곱 추정치, 즉 \(\| \boldsymbol{y} - \boldsymbol{\eta} (\boldsymbol{\beta}) \|^2\)를 최소화하는 \(\boldsymbol{\beta}\)

  • 분산 모수 \(\sigma^2\)는 선형회귀에서와 같이 잔차 제곱 평균으로 추정

  • 벡터 \(\boldsymbol{\eta} (\boldsymbol{\beta}) \in \mathbb{R}^n\)의 자취는 \(p\)차원 곡면

    • 선형회귀에서는 평면

  • 모형 적합: 유클리드 거리의 의미에서 표본 반응 벡터 \(\boldsymbol{y}\)에 가장 가까운 자취상의 좌표를 찾음

비선형 모형 적합

  • 체중 감소 모형 적합: nls() 사용

    wtloss_st <- c(b0 = 90, b1 = 95, th = 120)
wtloss_fm <- nls(Weight ~ b0 + b1*2^(-Days/th), 
                data = MASS::wtloss, 
                start = wtloss_st, trace = T)
    67.54349    (8.51e-01): par = (90 95 120)
40.18081    (1.55e-01): par = (82.72629 101.3046 138.7137)
39.24489    (2.20e-03): par = (81.39868 102.6584 141.8586)
39.24470    (4.39e-06): par = (81.37375 102.6842 141.9105)
    wtloss_fm
    Nonlinear regression model
  model: Weight ~ b0 + b1 * 2^(-Days/th)
   data: MASS::wtloss
    b0     b1     th 
 81.37 102.68 141.91 
 residual sum-of-squares: 39.24

Number of iterations to convergence: 3 
Achieved convergence tolerance: 4.389e-06
  • 선형 회귀 모형 적합과의 차이
    1. 추정에 대한 명시적인 공식이 없으므로 반복 절차 및 초기값이 필요
    2. 선형 모형 공식은 모형 행렬만 정의하므로 비선형 회귀 모형을 지정하는 데 적합하지 않음

nls()

  • formula: 비선형 모형 공식. response ~ mean 형식
    • 표준 형식: 모수와 설명변수 모두를 포함하는 일반적인 대수식. 연산자는 일반적인 산술적 의미를 가짐
    • 비어있는 좌변은 0으로 간주
  • data: 데이터 프레임
  • start: 모수의 초기값을 지정해주는 list 또는 수치형 벡터.
  • control: 기본 반복 절차의 파라미터 변경에 사용
  • algorithm: 문자열 인수. 특정 적합 알고리즘을 명시함. default 절차는 Gauss-Newton 알고리즘.
  • trace: 반복 절차 출력 여부
  • wtloss 예제에서 모수의 이름은 b0, b1, th와 같이 명시됨
  • 초기값 90, 95, 120은 Figure 1 로부터 눈대중으로 선택
  • trace의 결과로부터 절차가 3번 반복만에 수렴한 것으로 보인다.

가중치

  • nls()에는 weights 인수가 없음!

  • 가중치 W를 알면 비선형 회귀 공식을 y ~ M 대신 ~sqrt(W) * (y-M)을 써서 처리

  • W(의 일부)도 추정해야 한다면 보다 일반적인 최적화 기법을 이용해 적합해야 함

Gauss-Newton 알고리즘

  • \(\boldsymbol{\eta}(\boldsymbol{\beta}) = (\eta_1(\boldsymbol{\beta})), \dotsc, \eta_n(\boldsymbol{\beta})\)로 쓰자.

  • 초기값 \(\boldsymbol{\beta}^{(0)}\)에서 \(\eta_k\)의 1차 테일러 전개: \[ \eta_k(\boldsymbol{\beta}) \approx \eta_k \left(\boldsymbol{\beta}^{(0)}\right) +\underset{j=1}{\stackrel{p}{\sum}} \left(\beta_j-\beta_j^{(0)}\right) \frac{\partial \eta_k}{\partial \beta_j} \Bigg| _{\boldsymbol{\beta}=\boldsymbol{\beta}^{(0)} } \] 벡터 형태: \(\boldsymbol{\eta}(\boldsymbol{\beta}) \approx \boldsymbol{\omega}^{(0)}+\boldsymbol{Z}^{(0)} \boldsymbol{\beta}\), \[ Z_{k j}^{(0)}=\left.\frac{\partial \eta_k}{\partial \beta_j}\right|_{\boldsymbol{\beta}{=} \boldsymbol{\beta}^{(0)}} , \quad \omega_k^{(0)}=\eta_k\left(\boldsymbol{\beta}^{(0)}\right)-\sum_{j=1}^p \beta_j^{(0)} Z_{k j}^{(0)} \]

  • 표본 반응 벡터 \(\boldsymbol{y}\)를 오프셋 \(\boldsymbol{\omega}^{(0)}\)와 가모형행렬 \(Z^{(0)}\)로 정의되는 접평면으로 회귀하여 새로운 \(\boldsymbol{\beta}\)의 근사값 계산: \[ \boldsymbol{\beta}^{(1)}=\left((\boldsymbol{Z}^{(0)})^\top Z^{(0)}\right)^{-1} (\boldsymbol{Z}^{(0)})^{\top}\left(\boldsymbol{y}-\boldsymbol{\omega}^{(0)}\right) \]
  • 반복
  • 선형 회귀의 경우 오프셋이 \(\boldsymbol{0}\)이고 \(\boldsymbol{Z}^{(0)} = \boldsymbol{X}\) (모형 행렬)이므로 이 절차가 한 단계만에 수렴

도함수 정보 이용

  • 도함수 정보를 nls()에 알려주지 않으면 \(Z\) 행렬은 수치적 근사로 계산
  • 도함수 정보가 있으면 대체로 수렴성이 향상
  • 평균 벡터 \(\boldsymbol{\eta}\)를 계산하는 함수를 작성하고 \(Z\) 행렬을 gradient 속성으로 포함

wtloss 예제

  • 도함수 계산 \[ \frac{\partial \eta}{\partial \beta_0}=1, \quad \frac{\partial \eta}{\partial \beta_1}=2^{-x / \theta}, \quad \frac{\partial \eta}{\partial \theta}=\frac{\log (2) \beta_1 x 2^{-x / \theta}}{\theta^2} \]

    expn <- function(b0, b1, th, x){ 
    temp  <-  2^(-x/th)
    model_func <- b0  +  b1 * temp
    Z <- cbind(1, temp, (b1 * x  * temp  * log(2))/th^2)
    dimnames(Z) <- list(NULL, c("b0", "b1", "th")) 
    attr(model_func, "gradient") <- Z
    model_func
}

  • gradient 행렬의 열 이름은 대응하는 모수와 같아야 함

  • 도함수 정보를 이용한 모형 적합

    wtloss_gr <- nls(Weight ~  expn(b0, b1, th, Days), 
                 data = MASS::wtloss, 
                 start = wtloss_st, trace = T)
    67.54349    (8.51e-01): par = (90 95 120)
40.18081    (1.55e-01): par = (82.72629 101.3046 138.7137)
39.24489    (2.20e-03): par = (81.39868 102.6584 141.8586)
39.24470    (4.39e-06): par = (81.37375 102.6842 141.9105)

  • 수렴 속도에 차이가 없는 것처럼 보이지만 함수 expn()을 추적하면 도함수 정보를 줄 때는 6번의 평가만 필요한 반면, 그렇지 않을 때는 21번의 평가가 필요하다는 것을 알 수 있다.

자동 미분

  • expn()과 같은 함수는 기호 미분 함수 deriv()를 사용하여 자동으로 생성될 수 있는 경우가 많음

  • wtloss 예제

    (expn1 <- deriv(y ~ b0 + bl * 2^(-x/th), c("b0", "b1", "th"), 
                function(b0, b1, th, x){}) )
    function (b0, b1, th, x) 
{
    .expr3 <- 2^(-x/th)
    .value <- b0 + bl * .expr3
    .grad <- array(0, c(length(.value), 3L), list(NULL, c("b0", 
        "b1", "th")))
    .grad[, "b0"] <- 1
    .grad[, "b1"] <- 0
    .grad[, "th"] <- bl * (.expr3 * (log(2) * (x/th^2)))
    attr(.value, "gradient") <- .grad
    .value
}
  • deriv() 인수
    1. 모형 공식 (좌변 생략 가능)
    2. 모수의 이름 벡터
    3. 반환값이 될, 인수 리스트가 지정된 빈 함수

비선형 적합 모형 객체와 메소드 함수

nls 객체

  • nls() 호출의 결과는 클래스 nls 객체로, 표준 메소드 함수들을 이용할 수 있다.

  • summary.nls()

    summary(wtloss_gr, cor = T)
    
Formula: Weight ~ expn(b0, b1, th, Days)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
b0   81.374      2.269   35.86   <2e-16 ***
b1  102.684      2.083   49.30   <2e-16 ***
th  141.911      5.295   26.80   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8949 on 49 degrees of freedom

Correlation of Parameter Estimates:
   b0    b1   
b1 -0.99      
th -0.99  0.96

Number of iterations to convergence: 3 
Achieved convergence tolerance: 4.391e-06

  • 이탈도 함수

    deviance.nls <- function(object) sum(summary(object)$residuals^2)  
deviance(wtloss_gr)
    [1] 39.2447

  • vcov.nls(): 평균 모수의 공분산 행렬 추정량

    vcov(wtloss_gr)
               b0        b1        th
b0   5.148412 -4.674526 -11.84142
b1  -4.674526  4.337916  10.54321
th -11.841415 10.543211  28.03203

plinear: 부분 선형 최소제곱 알고리즘 (Golub and Pereyra, 1973)

  • 모든 비선형 모수가 알려진 경우 모형은 선형이 되며 표준 선형 회귀 방법을 사용할 수 있다.
  • 선형 모수를 활용하지 않는 방법보다 훨씬 더 안정적일 수 있고 초기 값이 더 적게 필요하며 다른 절차가 실패하는 시작점에서도 수렴되는 경우가 많다.

MASS::muscle 자료 (Linder, Chakravarti and Vuagnat, 1964)

as.tibble(MASS::muscle)
# A tibble: 60 × 3
   Strip  Conc Length
   <fct> <dbl>  <dbl>
 1 S01    1      15.8
 2 S01    2      20.8
 3 S01    3      22.6
 4 S01    4      23.8
 5 S02    1      20.6
 6 S02    2      26.8
 7 S02    3      28.4
 8 S02    4      27  
 9 S03    0.25    7.2
10 S03    0.5    15.4
# ℹ 50 more rows
length(levels(muscle$Strip))
[1] 21
  • 실험 목적: 쥐의 심장 근육 수축에 대한 용액 내 칼슘의 영향 평가
  • 쥐 21마리의 심장에서 좌심방을 분리하고 여러 차례 전기 자극을 가한 후 다양한 농도의 염화칼슘 용액에 담근 후 수축을 측정
  • 변수
    • Strip: 어떤 심장 근육 조각을 사용했나?
    • Conc: 염화칼슘 용액의 농도(2.2mM의 배수)
    • Length: 조각의 길이 단축 정도, mm 단위 (추정)
  • 모형 \[ \log y_{ij}=\alpha_j+\beta \rho^{x_{i j}}+\varepsilon_{ij} \]
    • \(i\): 염화칼슘 농도, \(j\): 근육 조각
    • 1개의 비선형 모수(\(\rho\)), 22개의 선형 모수(\(\alpha_j\), \(\beta\))
    • 부분 선형 알고리즘이 매우 편리한 예
  • nls()에서 algorithm = plinear를 선택하면 모형 공식의 우변은 열이 비선형 모수의 함수인 행렬을 지정한다. 선형 모수는 행렬의 열에 대한 회귀 계수가 된다. 선형 모형의 경우와 달리 내재적인 절편항은 없다.
    • 초기값은 비선형 모수에만 필요: \(\rho^{(0)} = 0.1\)
A <- model.matrix(~ Strip - 1, data = MASS::muscle)
rats.nlsl <-  nls(log(Length) ~ cbind(A, rho^Conc), data = MASS::muscle,  start = c(rho = 0.1), algorithm = "plinear") 
(B  <-  coef(rats.nlsl))
          rho .lin.StripS01 .lin.StripS02 .lin.StripS03 .lin.StripS04 
   0.07776401    3.08304824    3.30137838    3.44562531    2.80464434 
.lin.StripS05 .lin.StripS06 .lin.StripS07 .lin.StripS08 .lin.StripS09 
   2.60835015    3.03357725    3.52301734    3.38711844    3.46709396 
.lin.StripS10 .lin.StripS11 .lin.StripS12 .lin.StripS13 .lin.StripS14 
   3.81438456    3.73878664    3.51332581    3.39741115    3.47088608 
.lin.StripS15 .lin.StripS16 .lin.StripS17 .lin.StripS18 .lin.StripS19 
   3.72895847    3.31863862    3.37938673    2.96452195    3.58468686 
.lin.StripS20 .lin.StripS21        .lin22 
   3.39628029    3.36998872   -2.96015460

  • plinear 적합을 표준 Gauss-Newton 알고리즘의 초기값으로 사용

    st <- list(alpha = B[2:22], beta = B[23], rho = B[1])
rats.nls2 <- nls(log(Length) ~  alpha[Strip] + beta*rho^Conc,
                  data = MASS::muscle, start = st)
rats.nls2
    Nonlinear regression model
  model: log(Length) ~ alpha[Strip] + beta * rho^Conc
   data: MASS::muscle
alpha..lin.StripS01 alpha..lin.StripS02 alpha..lin.StripS03 alpha..lin.StripS04 
            3.08305             3.30138             3.44563             2.80464 
alpha..lin.StripS05 alpha..lin.StripS06 alpha..lin.StripS07 alpha..lin.StripS08 
            2.60835             3.03358             3.52302             3.38712 
alpha..lin.StripS09 alpha..lin.StripS10 alpha..lin.StripS11 alpha..lin.StripS12 
            3.46709             3.81438             3.73879             3.51333 
alpha..lin.StripS13 alpha..lin.StripS14 alpha..lin.StripS15 alpha..lin.StripS16 
            3.39741             3.47089             3.72896             3.31864 
alpha..lin.StripS17 alpha..lin.StripS18 alpha..lin.StripS19 alpha..lin.StripS20 
            3.37939             2.96452             3.58469             3.39628 
alpha..lin.StripS21         beta..lin22             rho.rho 
            3.36999            -2.96015             0.07776 
 residual sum-of-squares: 1.045

Number of iterations to convergence: 0 
Achieved convergence tolerance: 4.918e-06

  • 예상대로 한 번에 수렴

  • plinear 형식이 아닌 표준 형식으로 적합된 모형 객체에는 generic function predict()가 사용 가능

    attach(MASS::muscle)
# all possible predictors
Muscle <- expand.grid(Conc = sort(unique(Conc)), Strip = levels(Strip))
Muscle$Yhat <- predict(rats.nls2, Muscle)  # prediction
Muscle$logLength <- rep(NA, nrow(Muscle))
# observed responses only
ind <- match(paste(Strip, Conc), paste (Muscle$Strip,  Muscle$Conc))
Muscle$logLength[ind] <- log(Length)
detach(MASS::muscle)
    ggplot(Muscle, aes(x = Conc, y = Yhat)) + 
facet_wrap(~ Strip, ncol=6) +
geom_smooth(method = "loess", se = FALSE, na.rm = TRUE) +  
geom_point(aes(x=Conc, y = logLength)) +  
labs(x = "Calcium Chloride concentration (mM)", 
     y = "log(Length in mm)")

  • 모형은 상황을 상당히 잘 설명하는 것처럼 보이지만 일부 목적에서는 비선형 혼합 효과 모형이 더 적합할 수 있다.

모수의 신뢰구간

잉여 제곱합

  • 모수 벡터: \(\boldsymbol{\theta}=(\theta_1, \boldsymbol{\theta}_2)^T\)
  • 귀무가설 \(H_0 : \theta_1 = \theta_{10}\)
  • \(\widehat{\boldsymbol{\theta}}\): 전체 최소제곱 추정치
  • \(\widehat{\boldsymbol{\theta}}_{2|1}\): \(\theta_1 = \theta_{10}\)로 둔 \(\boldsymbol{\theta}_2\)의 조건부 최소제곱 추정치
    • \(\widehat{\boldsymbol{\theta}} (\theta_{10}) = (\theta_{10}, \widehat{\boldsymbol{\theta}}_{2|1})^T\)

  • 잉여 제곱합 기반 검정 통계량 \[ F(\theta_{10}) = [RSS(\widehat{\boldsymbol{\theta}} (\theta_{10}))-RSS(\widehat{\boldsymbol{\theta}})]/{s^2} \]

    • \(H_0\)하에서 근사적으로 \(F_{1, n-p}\) 분포

  • 부호화 제곱근 \[ \tau(\theta_{10}) = sign(\theta_{10}-\widehat{\theta}_1) \sqrt{F(\theta_{10})} \]

    • \(H_0\)하에서 근사적으로 \(t_{n-p}\) 분포

  • \(\theta_1\)에 대한 신뢰구간: \(\{ \theta_1 | -t_{n-p}(1-\alpha/2) < \tau(\theta_1) < t_{n-p}(1-\alpha/2) \}\)

  • MASS::confint()로 계산 가능

wtloss 자료

  • 체중 감소 프로그램에 참여하는 사람에게 중요한 질문

목표 체중에 도달하기까지 시간이 얼마나 걸릴 것인가?

  • 모형 \(y = \beta_0 + \beta_1 2^{-t/\theta}\) 하에서 체중 \(w_0\)을 달성하는 데 걸린 시간 \[ \delta_0=-\theta \log_2 \{ (w_0 - \beta_0)/\beta_1 \} \]

  • 모형을 \(\delta_0\)에 대해 다시 쓰면 \[ y=\beta_0 + \beta_1 \left( \frac{w_0-\beta_0}{\beta_1} \right) ^{x/\delta_0} + \epsilon \]

  • R 표현

    expn2 <- deriv(~b0 + b1*((w0 - b0)/b1)^(x/d0),
              c("b0","b1","d0"), function(b0, b1, d0, x, w0) {})

  • \(\delta_0\)를 모수로 하는 적합

    wtloss.init <- function(obj, w0) {
    p <- coef(obj)
    d0 <-  - log((w0 - p["b0"])/p["b1"])/log(2) * p["th"]
    c(p[c("b0", "b1")], d0 = as.vector(d0))
}
out <- NULL
w0s <- c(110, 100, 90)
for(w0 in w0s) {
    fm <- nls(Weight ~ expn2(b0, b1, d0, Days, w0),
             wtloss, start = wtloss.init(wtloss_gr, w0))
    out <- rbind(out, c(coef(fm)["d0"], confint(fm, "d0")))
}
dimnames(out)[[1]] <- paste(w0s,"kg:")
out
              d0     2.5%    97.5%
110 kg: 261.5132 256.2303 267.5009
100 kg: 349.4979 334.7293 368.0151
90 kg:  507.0941 457.2637 594.8745
  • 체중이 감소함에 따라 신뢰 구간은 훨씬 넓어지고, 비대칭이 커진다. 추정된 점근 체중에 더 가까운 체중의 경우 프로파일링 절차가 실패할 수 있다.

신뢰영역: stormer 자료

  • Stormer 점도계는 내부의 실린더가 작동 중량에 반응하여 정해진 회전수를 수행하는 데 걸리는 시간을 측정하여 유체의 점도를 측정

  • 점도계 보정은 장치가 점도를 정확히 알고 있는 유체에 정지 상태로 있는 동안 다양한 중량에 대해 소요되는 시간을 측정하여 이루어짐

Weight Viscosity (poise)
(grams) 14.7 27.5 42.0 75.7 89.7 146.6 158.3 161.1 298.3
20 35.6 54.3 75.6 121.2 150.8 229.0 270.0
50 17.6 24.3 31.4 47.2 58.3 85.6 101.1 92.2 187.2
100 24.6 30.0 41.7 50.3 45.1 \(89.0,86.5\)

  • 비선형 회귀모형 \[ T=\frac{\beta_1 v}{w - \beta_2}+\epsilon \]

    • \(T\): 소요시간, \(w\): 중량, \(v\) 점도
    • \(\beta_1\): 선형 모수, \(\beta_2\): 비선형 모수

  • OLS 초기값 (Williams, 1959) \[ \omega T = \beta_1 v + \beta_2 T + (\omega-\beta_2) \epsilon \]

    fm0 <- lm(Wt*Time ~ Viscosity + Time -  1, data = MASS::stormer)
b0 <- coef(fm0); names(b0) <- c("b1", "b2"); b0
           b1        b2 
28.875541  2.843728

  • NLS 회귀

    storm_fm <- nls(Time ~ b1*Viscosity/(Wt-b2), data = stormer,
                start = b0,  trace = T)
    885.3645    (2.71e-01): par = (28.87554 2.843728)
825.1098    (8.42e-03): par = (29.39346 2.233276)
825.0514    (1.69e-05): par = (29.40133 2.218226)
825.0514    (5.24e-08): par = (29.40126 2.218274)

회귀 모수의 신뢰 영역

bc <- coef(storm_fm)
se <- sqrt(diag(vcov(storm_fm)))
dv <- deviance(storm_fm)

\[ \small RSS\left(\beta_1, \beta_2\right)=\sum_{i=1}^{23}\left(T_i-\frac{\beta_1 v_i}{w_i-\beta_2}\right)^2 \]

  • dv == \(d_0 = RSS(\widehat{\beta}_1, \widehat{\beta}_2)\)

  • 잉여 제곱합 통계량 \[ \small F(\beta_1^*, \beta_2^*) = \frac{(RSS(\beta_1^*, \beta_2^*)-d_0)/2}{d_0/21} \stackrel{\cdot}{\sim} F_{2, ~ 21} \]

  • 신뢰영역: \(F(\beta_1, \beta_2)\)가 \(F_{2, ~ 21}\) 분포의 95백분위수보다 작은 \((\beta_1, \beta_2)\)의 집합

par(pty = "s")  # square plotting region
gridsize = 51
b1 <- bc[1] + seq(-3*se[1] , 3*se[1] , length = gridsize)  # 3 * std
b2 <- bc[2] + seq(-3*se[2] , 3*se[2] , length = gridsize) 
bv <- expand.grid(b1, b2)
attach(stormer)
cc <- matrix(Time - rep(bv[,1],rep(nrow(stormer), gridsize*gridsize)) * 
              Viscosity/(Wt - rep(bv[,2], rep(nrow(stormer), gridsize*gridsize))), 
            23)
dbetas <- matrix(drop(rep(1, 23) %*% cc^2), gridsize)   # RSS
fstat <-  matrix( ((dbetas - dv)/2) / (dv/21), gridsize, gridsize)  # F-statistic
detach(stormer)
plot(b1, b2, type = "n")
lev <-  c(1,  2,  5,  7,  10,  15,  20)
contour(b1, b2,  fstat, levels = lev,  labex = 0.75,  lty = 2,  add = T)
# 95%-ile
contour(b1, b2,  fstat, levels = qf(0.95,2,21), add = T,  labex = 0)
text(31.6, 0.3, labels = "95% CR", adj = 0, cex = 0.75)
points(bc[1], bc[2], pch = 3, mkh = 0.1)

  • 가능도 함수가 \(F\)-통계량과 동일한 등고선을 갖기 때문에 타원형에 가까운 등고선의 모양은 정규 선형 회귀에 기반한 근사가 잘 맞을 것임을 시사

  • 등고선의 길쭉한 모양은 추정치 \(\widehat{\beta}_1\), \(\widehat{\beta}_2\)가 높은 (음의) 상관 관계가 있음을 보여줌

  • 다른 여러 모수가 존재하는 경우 두 개의 회귀 모수에 대한 이변량 신뢰 영역을 구하는 것은 적어도 원칙적으로는 어려운 계산일 수 있다. 왜냐하면 등고선을 그리는 \(F\)-통계 표면의 각 지점은 다른 모수에 대해 최적화하여 계산해야 하기 때문이다.

부트스트랩

  • stormer 자료는 설계된 실험 결과이므로 라이브러리 boot의 함수를 사용하여 모형 기반 부트스트랩을 행하여 모수 추정치의 분포를 탐색할 수 있다.

  • 비선형 모형이므로 잔차의 평균이 \(0\)이 아닐 수 있어 이를 제거하고 사용

summary(storm_fm)$parameters
    Estimate Std. Error   t value     Pr(>|t|)
b1 29.401257  0.9155335 32.113795 2.457598e-19
b2  2.218274  0.6655216  3.333136 3.155740e-03
  • 신뢰영역 도표로부터 가능도 기반 분석이 백분위수 구간을 뒷받침할 것임을 시사
st <-  cbind(stormer, fit = fitted(storm_fm))
storm_bf <-  function(rs, i){
    st$Time <- st$fit + rs[i]
    coef(nls(Time ~ b1 * Viscosity/(Wt - b2), st, start = coef(storm_fm)))
}
rs <- scale(resid(storm_fm), scale = F) # remove the mean
(storm.boot <- boot::boot(rs, storm_bf, R = 9999)) 

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot::boot(data = rs, statistic = storm_bf, R = 9999)


Bootstrap Statistics :
     original     bias    std. error
t1* 28.715851  0.6791065   0.8348826
t2*  2.479781 -0.2610913   0.6075570
boot.ci(storm.boot, index = 1,
        type = c("norm", "basic", "pere", "bca"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 9999 bootstrap replicates

CALL : 
boot.ci(boot.out = storm.boot, type = c("norm", "basic", "pere", 
    "bca"), index = 1)

Intervals : 
Level      Normal              Basic                BCa          
95%   (26.40, 29.67 )   (26.45, 29.73 )   (26.12, 29.66 )  
Calculations and Intervals on Original Scale
Warning : BCa Intervals used Extreme Quantiles
Some BCa intervals may be unstable
boot.ci(storm.boot,  index = 2,
        type =  c("norm",  "basic",  "pere",  "bca"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 9999 bootstrap replicates

CALL : 
boot.ci(boot.out = storm.boot, type = c("norm", "basic", "pere", 
    "bca"), index = 2)

Intervals : 
Level      Normal              Basic                BCa          
95%   ( 1.550,  3.932 )   ( 1.556,  3.948 )   ( 1.567,  3.910 )  
Calculations and Intervals on Original Scale

프로파일

RSS의 저차원 프로파일

  • 선형 근사가 정확한 좌표 방향의 경우, \(\beta_i\)에 대한 비선형 t 통계량 \(\tau(\beta_i)\)의 그래프는 최대가능도 추정치 근방에서 직선이 되어야 함

  • 직선에서 크게 벗어난다면 선형 근사가 이 방향으로는 맞지 않음을 경고

plot(profile(wtloss_gr), absVal=FALSE)