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

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
  }

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 점도계는 내부의 실린더가 작동 중량에 반응하여 정해진 회전수를 수행하는 데 걸리는 시간을 측정하여 유체의 점도를 측정

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

• 비선형 회귀모형 \[ 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)\)의 집합

• 가능도 함수가 \(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

• 신뢰영역 도표로부터 가능도 기반 분석이 백분위수 구간을 뒷받침할 것임을 시사

RSS의 저차원 프로파일

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

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

plot(profile(wtloss_gr), absVal=FALSE)