Latihan dan Tugas

1 Konsep dan Teori Analisis Regresi Tingkat Lanjut

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, dpi = 150,
fig.align = "center", fig.width = 7, fig.height = 4.5)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(broom)
library(sandwich)   # robust/sandwich SE
library(lmtest)     # coeftest

## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(car)        # VIF

## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some

library(modelr)     # CV utilities

## 
## Attaching package: 'modelr'
## 
## The following object is masked from 'package:broom':
## 
##     bootstrap

library(splines)    # basis spline
set.seed(123)

1.1 Konsep Analisis Regresi dan Korelasi

Analisis regresi adalah metode statistik yang digunakan untuk mempelajari hubungan satu variabel dependen (respon) dengan satu atau lebih variabel independen (prediktor). Tujuan utamanya adalah membangun model matematis yang bisa digunakan untuk:

1. Menjelaskan hubungan antara variabel.
2. Memprediksi nilai variabel dependen berdasarkan variabel independen.
3. Mengukur pengaruh masing-masing prediktor terhadap respon.

Model regresi linear sederhana dituliskan sebagai: \[ Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i, \quad i = 1, \ldots, n \] dengan:

• \(Y_i\) : variabel dependen (respon)

• \(X_i\) : variabel independen (prediktor)

• \(\beta_0, \beta_1\) : parameter regresi

• \(\varepsilon_i\) : error (residu)

Konsep Korelasi

Korelasi mengukur derajat keeratan hubungan linier antara dua variabel, tanpa membedakan mana yang dependen dan independen.

Koefisien korelasi Pearson:

Koefisien korelasi Pearson dirumuskan sebagai:

\[ r = \frac{\sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i=1}^{n} (x_i - \bar{x})^2 \sum_{i=1}^{n} (y_i - \bar{y})^2}} \]

• Nilai \(r \in [-1, 1]\).
  
• \(r > 0\): hubungan positif, \(r < 0\): hubungan negatif.
  
• \(|r|\) makin mendekati 1 artinya hubungan makin kuat.

Perbedaan Utama

\[ \begin{array}{|l|l|l|} \hline \textbf{Aspek} & \textbf{Korelasi} & \textbf{Regresi} \\ \hline \text{Tujuan} & \text{Mengukur keeratan hubungan} & \text{Membangun model prediksi/penjelasan} \\ \hline \text{Arah hubungan} & \text{Simetris (X dan Y setara)} & \text{Asimetris (Y sebagai respon, X sebagai prediktor)} \\ \hline \text{Output utama} & \text{Koefisien korelasi } r & \text{Persamaan regresi } \hat{Y} = \beta_0 + \beta_1 X \\ \hline \end{array} \]

Contoh Aplikasi di R

Data

Gunakan dataset mtcars untuk melihat hubungan mpg (miles per gallon) dengan wt (berat mobil).

data(mtcars)
head(mtcars[, c("mpg", "wt")])

##                    mpg    wt
## Mazda RX4         21.0 2.620
## Mazda RX4 Wag     21.0 2.875
## Datsun 710        22.8 2.320
## Hornet 4 Drive    21.4 3.215
## Hornet Sportabout 18.7 3.440
## Valiant           18.1 3.460

Korelasi

cor.test(mtcars$wt, mtcars$mpg)

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$mpg
## t = -9.559, df = 30, p-value = 1.294e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9338264 -0.7440872
## sample estimates:
##        cor 
## -0.8676594

Regresi

model <- lm(mpg~ wt, data = mtcars)
summary(model)

## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Visualisasi

library(ggplot2)
ggplot(mtcars, aes(x = wt, y = mpg)) +
geom_point() +
geom_smooth(method = "lm", se = TRUE, color = "blue") +
labs(title = "Hubungan antara Berat Mobil (wt) dan Efisiensi Bahan Bakar (mpg)",
x = "Berat Mobil (1000 lbs)",
y = "Miles per Gallon (mpg)")

1.2 Model Regresi Linear

Kita memodelkan hubungan antara respons \(y \in \mathbb{R}^n\) dan kovariat \(X \in \mathbb{R}^{n \times p}\): \[ y = X\beta + \varepsilon, \qquad \mathbb{E}(\varepsilon) = 0, \qquad \text{Var}(\varepsilon) = \sigma^2 I_n. \]

Dengan \(X = [1, x_1, \dots, x_{p-1}]\) (termasuk intersep), \(\beta\) adalah vektor koefisien, dan \(\varepsilon\) adalah galat.

1.2.1 Estimator OLS & Geometri Proyeksi

Estimator Ordinary Least Squares (OLS) meminimalkan jumlah kuadrat residu: \[ \hat{\beta}_{\text{OLS}} = \arg \min_{\beta} \| y - X\beta \|^2 = (X^\top X)^{-1} X^\top y. \]

Prediksi \(\hat{y} = X\hat{\beta}\) adalah proyeksi ortogonal dari \(y\) ke ruang kolom \(X\). Matriks hat dan residu: \[ H = X(X^\top X)^{-1} X^\top, \qquad \hat{y} = Hy, \qquad e = (I - H)y. \]

Leverage pengamatan ke-\(i\) adalah \(h_{ii}\) (elemen diagonal \(H\)).

1.2.2 Sifat Asimtotik & Varian Koefisien

Di bawah asumsi Gauss–Markov (linearitas, eksogenitas, homoskedastisitas, dan tidak ada multikolinearitas sempurna), \(\hat{\beta}_{\text{OLS}}\) adalah BLUE. Varians koefisien dan estimator ragam:

\[ \text{Var}(\hat{\beta}_{\text{OLS}}) = \sigma^2 (X^\top X)^{-1}, \qquad \hat{\sigma}^2 = \frac{\|e\|^2}{n - p}. \]

Asimtotik: \(\sqrt{n} (\hat{\beta} - \beta) \overset{d}{\longrightarrow} \mathcal{N}(0, V)\), sehingga inferensi berbasis normal/\(t\) berlaku untuk \(n\) besar.

1.3 Generalisasi: WLS & GLS

Jika galat tidak homoskedastik/berkorelasi, maka \(\text{Var}(\varepsilon) = \sigma^2 \Omega \neq \sigma^2 I\).

Weighted least squares (WLS) untuk varians berbeda-beda (\(\Omega\) diagonal): \[ \hat{\beta}_{\text{WLS}} = (X^\top W X)^{-1} X^\top W y, \qquad W = \Omega^{-1}. \]

Generalized least squares (GLS) (\(\Omega\) umum): \[ \hat{\beta}_{\text{GLS}} = (X^\top \Omega^{-1} X)^{-1} X^\top \Omega^{-1} y. \]

1.4 Robust Standard Errors (White Sandwich)

Tanpa menspesifikasi bentuk \(\Omega\), kovarians heteroskedasticity-consistent (HC) ala White adalah: \[ \widehat{\text{Var}}(\hat{\beta}) = (X^\top X)^{-1} \left( \sum_{i=1}^n x_i x_i^\top e_i^2 \right) (X^\top X)^{-1}. \]

Ini memberi robust SE yang valid terhadap heteroskedastisitas (HC0–HC3).

1.5 Multikolinearitas & Diagnostik

Variance Inflation Factor (VIF) untuk koefisien ke-\(j\): \[ \text{VIF}_j = \frac{1}{1 - R_j^2}. \]

dengan \(R_j^2\) adalah koefisien determinasi dari regresi \(x_j\) pada kovariat lain. Nilai besar (\(\gtrsim 10\)) mengindikasikan multikolinearitas kuat.

1.6 Diagnostik Residu, Leverage, dan Pengaruh

Ukuran pengaruh Cook’s distance: \[ D_i = \frac{e_i^2}{p \hat{\sigma}^2} \cdot \frac{h_{ii}}{(1 - h_{ii})^2}. \]

Gunakan plot residu, plot leverage vs residu, QQ-plot, serta inspeksi \(D_i\) untuk memeriksa asumsi dan outlier.

1.7 Bias–Variance Trade-off & Regularisasi (Ikhtisar)

Model makin kompleks (banyak parameter/basis) cenderung bias rendah namun varians tinggi. Ridge dan Lasso menambahkan penalti:

\[ \text{Ridge: } \hat{\beta} = \arg \min_{\beta} \| y - X\beta \|^2 + \lambda \|\beta\|_2^2, \qquad \text{Lasso: } \hat{\beta} = \arg \min_{\beta} \| y - X\beta \|^2 + \lambda \|\beta\|_1. \]

1.8 Nonlinearitas: Ekspansi Basis & Spline

Kita dapat menambah term polinomial atau basis spline sehingga tetap linear dalam parameter: \[ y = \beta_0 + \sum_k \beta_k b_k(x) + \varepsilon. \]

1.9 Studi Kasus 1: Regresi pada Data Nyata (mtcars)

1.9.1 Persiapan Data dan Eksplorasi

dat1 <- mtcars %>%
as_tibble(rownames = "car") %>%
mutate(across(c(mpg, disp, hp, wt, qsec), as.double))
# Plot hubungan awal
GG1 <- ggplot(dat1, aes(wt, mpg)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = FALSE) +
labs(x = "Weight (1000 lbs)", y = "Miles per Gallon (mpg)",
title = "Hubungan mpg vs wt pada mtcars")
GG1

1.9.2 Estimasi OLS & Ringkasan

fit1 <- lm(mpg~ wt + hp + disp, data = dat1)
summary(fit1)

## 
## Call:
## lm(formula = mpg ~ wt + hp + disp, data = dat1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.891 -1.640 -0.172  1.061  5.861 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.105505   2.110815  17.579  < 2e-16 ***
## wt          -3.800891   1.066191  -3.565  0.00133 ** 
## hp          -0.031157   0.011436  -2.724  0.01097 *  
## disp        -0.000937   0.010350  -0.091  0.92851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.639 on 28 degrees of freedom
## Multiple R-squared:  0.8268, Adjusted R-squared:  0.8083 
## F-statistic: 44.57 on 3 and 28 DF,  p-value: 8.65e-11

# Koefisien dengan CI 95%
tidy(fit1, conf.int = TRUE)

## # A tibble: 4 × 7
##   term         estimate std.error statistic  p.value conf.low conf.high
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept) 37.1         2.11     17.6    1.16e-16  32.8     41.4    
## 2 wt          -3.80        1.07     -3.56   1.33e- 3  -5.98    -1.62   
## 3 hp          -0.0312      0.0114   -2.72   1.10e- 2  -0.0546  -0.00773
## 4 disp        -0.000937    0.0103   -0.0905 9.29e- 1  -0.0221   0.0203

1.9.3 Diagnostik Dasar (Residu, QQ-plot, Cook)

diag_df <- tibble(
fitted = fitted(fit1), resid = resid(fit1),
stdres = rstandard(fit1), hat = hatvalues(fit1), cook = cooks.distance(fit1)
)
p1 <- ggplot(diag_df, aes(fitted, resid)) +
geom_point(alpha = 0.7) +
geom_hline(yintercept = 0, linetype = 2) +
labs(x = "Fitted", y = "Residuals")
p1

Figure 1.1 Diagnostik OLS untuk mtcars: residu vs fitted (kiri) dan QQ-plot (kanan).

p2 <- ggplot(diag_df, aes(sample = stdres)) +
stat_qq() + stat_qq_line() + labs(x = "Theoretical", y = "Standardized Residuals")
p2

Figure 1.2 Diagnostik OLS untuk mtcars: residu vs fitted (kiri) dan QQ-plot (kanan).

qplot(seq_along(diag_df$cook), diag_df$cook, geom = c("line","point"),
xlab = "Index", ylab = "Cook's D")

Figure 1.3 Cook’s distance per observasi.

1.9.4 Multikolinearitas (VIF)

car::vif(fit1)

##       wt       hp     disp 
## 4.844618 2.736633 7.324517

1.9.5 Robust SE (White HC3) & Interpretasi

coeftest(fit1)        # SE OLS

## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) 37.10550527  2.11081525 17.5788 < 2.2e-16 ***
## wt          -3.80089058  1.06619064 -3.5649  0.001331 ** 
## hp          -0.03115655  0.01143579 -2.7245  0.010971 *  
## disp        -0.00093701  0.01034974 -0.0905  0.928507    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(fit1, vcov = vcovHC(fit1, type = "HC3")) # SE robust (HC3)

## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) 37.10550527  2.71026981 13.6907 6.257e-14 ***
## wt          -3.80089058  1.08086207 -3.5165   0.00151 ** 
## hp          -0.03115655  0.01437457 -2.1675   0.03886 *  
## disp        -0.00093701  0.01012138 -0.0926   0.92690    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.9.6 Nonlinearitas: Polinomial vs Spline pada wt

fit_poly <- lm(mpg~ poly(wt, 2, raw = TRUE) + hp + disp, data = dat1)
fit_spl <- lm(mpg~ bs(wt, df = 5) + hp + disp, data = dat1)
AIC(fit1, fit_poly, fit_spl) %>% arrange(AIC)

##          df      AIC
## fit_poly  6 150.7305
## fit_spl   9 154.3323
## fit1      5 158.6430

BIC(fit1, fit_poly, fit_spl) %>% arrange(BIC)

##          df      BIC
## fit_poly  6 159.5250
## fit1      5 165.9717
## fit_spl   9 167.5239

# Kurva prediksi
grid_wt <- tibble(wt = seq(min(dat1$wt), max(dat1$wt), length.out = 200),
hp = mean(dat1$hp), disp = mean(dat1$disp))
pred <- grid_wt %>%
mutate(OLS = predict(fit1, newdata = grid_wt),
Poly = predict(fit_poly, newdata = grid_wt),
Spline = predict(fit_spl, newdata = grid_wt)) %>%
pivot_longer(OLS:Spline, names_to = "model", values_to = "yhat")
GG2 <- ggplot() +
geom_point(data = dat1, aes(wt, mpg), alpha = 0.3) +
geom_line(data = pred, aes(wt, yhat, linetype = model), linewidth = 0.9) +
labs(x = "wt", y = "mpg", title = "Perbandingan OLS vs Polinomial vs Spline")
GG2

1.9.7 Validasi Silang (K-fold) untuk Perbandingan Model

k <- 5
cv_res <- crossv_kfold(dat1, k = k)
rmse_fun <- function(train, test, formula) {
m <- lm(formula, data = as.data.frame(train))
sqrt(mean((as.data.frame(test)$mpg- predict(m, newdata = as.data.frame(test)))^2))
}
forms <- list(
OLS = mpg~ wt + hp + disp,
Poly = mpg~ poly(wt, 2, raw = TRUE) + hp + disp,
Spline = mpg~ bs(wt, df = 5) + hp + disp
)
cv_tbl <- map_dfr(names(forms), function(nm){
fml <- forms[[nm]]
tibble(model = nm, fold = seq_len(k),
rmse = map2_dbl(cv_res$train, cv_res$test,~ rmse_fun(.x, .y, fml)))
})
cv_tbl %>% group_by(model) %>% summarise(mean_rmse = mean(rmse), sd_rmse = sd(rmse)) %>% arrange(mean_rmse)

## # A tibble: 3 × 3
##   model  mean_rmse sd_rmse
##   <chr>      <dbl>   <dbl>
## 1 Poly        2.35   0.980
## 2 Spline      2.74   0.807
## 3 OLS         2.92   1.01

1.10 Studi Kasus 2: Simulasi Heteroskedastisitas (WLS vs OLS

vs Robust)

n <- 400
x1 <- runif(n, 0, 10)
x2 <- rnorm(n, 5, 2)
sigma_i <- 0.6 + 0.2 * x1 # var naik dengan x1
err <- rnorm(n, 0, sigma_i)
y <- 2 + 1.5*x1- 1.2*x2 + err
sim <- tibble(y, x1, x2, w = 1/sigma_i^2)
fit_ols <- lm(y~ x1 + x2, data = sim)
fit_wls <- lm(y~ x1 + x2, data = sim, weights = w)
# Tabel perbandingan SE
comp_se <- tibble(
term = names(coef(fit_ols)),
OLS_SE = coef(summary(fit_ols))[, "Std. Error"],
Robust_SE = sqrt(diag(vcovHC(fit_ols, type = "HC3"))),
WLS_SE = coef(summary(fit_wls))[, "Std. Error"]
)
comp_se

## # A tibble: 3 × 4
##   term        OLS_SE Robust_SE WLS_SE
##   <chr>        <dbl>     <dbl>  <dbl>
## 1 (Intercept) 0.250     0.212  0.175 
## 2 x1          0.0293    0.0282 0.0261
## 3 x2          0.0415    0.0399 0.0313

# Uji Breusch-Pagan untuk heteroskedastisitas
lmtest::bptest(fit_ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  fit_ols
## BP = 29.415, df = 2, p-value = 4.098e-07

ggplot(data.frame(fitted=fitted(fit_ols), resid=resid(fit_ols)), aes(fitted, resid)) +
geom_point(alpha=.6) + geom_smooth(se=FALSE, method="loess") +
geom_hline(yintercept = 0, linetype=2) + labs(x="Fitted", y="Residuals")

Figure 1.4: Heteroskedastisitas: Residual vs Fitted (OLS).

1.11 Catatan

1. Mulai dari OLS + diagnostik; bila heteroskedastik, pilih robust SE atau WLS/GLS sesuai struktur ragam.
2. Periksa multikolinearitas (VIF), leverage, Cook’s D, dan normalitas residu.
3. Gunakan ekspansi basis (polinomial/spline) saat ada indikasi nonlinearitas.
4. Validasi model (CV) dan bandingkan AIC/BIC untuk kompleksitas vs akurasi.
5. Untuk dimensi tinggi, gunakan regularisasi (Ridge/Lasso) dan pipeline (standarisasi, CV).

1.12 GLS dengan Korelasi AR(1): Contoh dan Simulasi

Model galat mengikuti proses AR(1):
\[ \varepsilon_t = \rho \varepsilon_{t-1} + u_t, \quad \text{dengan } u_t \sim \mathcal{N}(0, \sigma_u^2). \]

Kovarians residual memiliki bentuk Toeplitz:
\[ \Omega_{ts} = \sigma_u^2 \rho^{|t-s|}/(1-\rho^2). \]

GLS menimbang data dengan \(\Omega^{-1}\) (atau transformasi ekuivalen) agar estimasi efisien.

1.12.1 Simulasi data dengan galat AR(1)

set.seed(42)
T <- 400
x <- runif(T,-2, 2)
rho <- 0.6
u <- rnorm(T, 0, 1)
err <- numeric(T)
for (t in 2:T) err[t] <- rho*err[t-1] + u[t]
y <- 1 + 2*x + err
dar1 <- tibble::tibble(t = 1:T, x, y)
head(dar1)

## # A tibble: 6 × 3
##       t       x     y
##   <int>   <dbl> <dbl>
## 1     1  1.66   4.32 
## 2     2  1.75   4.83 
## 3     3 -0.855  0.661
## 4     4  1.32   6.53 
## 5     5  0.567  2.49 
## 6     6  0.0764 0.213

1.12.2 OLS vs GLS (nlme::gls dengan corAR1)

library(nlme)
# OLS biasa
fit_ols_ar1 <- lm(y~ x, data = dar1)
# GLS dengan korelasi AR(1) di residual
fit_gls_ar1 <- gls(y~ x, data = dar1,
correlation = corAR1(form =~ t))
summary(fit_ols_ar1)

## 
## Call:
## lm(formula = y ~ x, data = dar1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3360 -0.8231  0.0026  0.8211  3.9013 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.95441    0.06043   15.79   <2e-16 ***
## x            2.06045    0.05162   39.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.209 on 398 degrees of freedom
## Multiple R-squared:  0.8002, Adjusted R-squared:  0.7997 
## F-statistic:  1594 on 1 and 398 DF,  p-value: < 2.2e-16

summary(fit_gls_ar1)

## Generalized least squares fit by REML
##   Model: y ~ x 
##   Data: dar1 
##        AIC      BIC    logLik
##   1139.541 1155.487 -565.7706
## 
## Correlation Structure: AR(1)
##  Formula: ~t 
##  Parameter estimate(s):
##       Phi 
## 0.5772332 
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 0.9552097 0.11651869  8.19791       0
## x           2.0449090 0.03743358 54.62766       0
## 
##  Correlation: 
##   (Intr)
## x 0.004 
## 
## Standardized residuals:
##          Min           Q1          Med           Q3          Max 
## -2.752477543 -0.674750151 -0.007333444  0.679546648  3.225627947 
## 
## Residual standard error: 1.210605 
## Degrees of freedom: 400 total; 398 residual

1.12.3 Diagnostik serial correlation

par(mfrow = c(1,2))
acf(resid(fit_ols_ar1), main = "ACF Residuals (OLS)")
acf(resid(fit_gls_ar1, type = "normalized"), main = "ACF Residuals (GLS)")

Figure 1.5: ACF residu dari OLS vs GLS

par(mfrow = c(1,1))
lmtest::dwtest(fit_ols_ar1) # indikasi korelasi serial pada OLS

## 
##  Durbin-Watson test
## 
## data:  fit_ols_ar1
## DW = 0.84993, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

1.12.4 Alternatif: Quasi-differencing (Cochrane–Orcutt) manual

# Estimasi rho dari regresi resid_t ~ resid_{t-1}
res <- resid(fit_ols_ar1)
rho_hat <- coef(lm(res[-1]~ 0 + res[-length(res)]))[1]
# Quasi-differencing (transformasi data)
y_star <- dar1$y[-1]- rho_hat*dar1$y[-T]
x_star <- dar1$x[-1]- rho_hat*dar1$x[-T]
fit_qd <- lm(y_star~ x_star)
summary(fit_qd)

## 
## Call:
## lm(formula = y_star ~ x_star)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6691 -0.6723 -0.0232  0.6970  3.2363 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.40603    0.04955   8.194 3.52e-15 ***
## x_star       2.04499    0.03757  54.439  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9898 on 397 degrees of freedom
## Multiple R-squared:  0.8819, Adjusted R-squared:  0.8816 
## F-statistic:  2964 on 1 and 397 DF,  p-value: < 2.2e-16

Catatan: gls() memperkirakan 𝜌 dan ragam secara simultan dalam kerangka likelihood sehingga sering lebih stabil dibanding pendekatan dua tahap.

1.13 Latihan

set.seed(123)
# Paket yang diperlukan
pkgs <- c("dplyr","ggplot2","car","lmtest","sandwich","MASS","boot","broom")
new <- pkgs[!(pkgs %in% installed.packages()[,"Package"])]
if(length(new)) install.packages(new, dependencies = TRUE)
lapply(pkgs, library, character.only = TRUE)

## [[1]]
##  [1] "nlme"      "splines"   "modelr"    "car"       "carData"   "lmtest"   
##  [7] "zoo"       "sandwich"  "broom"     "lubridate" "forcats"   "stringr"  
## [13] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [19] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [25] "methods"   "base"     
## 
## [[2]]
##  [1] "nlme"      "splines"   "modelr"    "car"       "carData"   "lmtest"   
##  [7] "zoo"       "sandwich"  "broom"     "lubridate" "forcats"   "stringr"  
## [13] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [19] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [25] "methods"   "base"     
## 
## [[3]]
##  [1] "nlme"      "splines"   "modelr"    "car"       "carData"   "lmtest"   
##  [7] "zoo"       "sandwich"  "broom"     "lubridate" "forcats"   "stringr"  
## [13] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [19] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [25] "methods"   "base"     
## 
## [[4]]
##  [1] "nlme"      "splines"   "modelr"    "car"       "carData"   "lmtest"   
##  [7] "zoo"       "sandwich"  "broom"     "lubridate" "forcats"   "stringr"  
## [13] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [19] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [25] "methods"   "base"     
## 
## [[5]]
##  [1] "nlme"      "splines"   "modelr"    "car"       "carData"   "lmtest"   
##  [7] "zoo"       "sandwich"  "broom"     "lubridate" "forcats"   "stringr"  
## [13] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [19] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [25] "methods"   "base"     
## 
## [[6]]
##  [1] "MASS"      "nlme"      "splines"   "modelr"    "car"       "carData"  
##  [7] "lmtest"    "zoo"       "sandwich"  "broom"     "lubridate" "forcats"  
## [13] "stringr"   "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"   
## [19] "ggplot2"   "tidyverse" "stats"     "graphics"  "grDevices" "utils"    
## [25] "datasets"  "methods"   "base"     
## 
## [[7]]
##  [1] "boot"      "MASS"      "nlme"      "splines"   "modelr"    "car"      
##  [7] "carData"   "lmtest"    "zoo"       "sandwich"  "broom"     "lubridate"
## [13] "forcats"   "stringr"   "dplyr"     "purrr"     "readr"     "tidyr"    
## [19] "tibble"    "ggplot2"   "tidyverse" "stats"     "graphics"  "grDevices"
## [25] "utils"     "datasets"  "methods"   "base"     
## 
## [[8]]
##  [1] "boot"      "MASS"      "nlme"      "splines"   "modelr"    "car"      
##  [7] "carData"   "lmtest"    "zoo"       "sandwich"  "broom"     "lubridate"
## [13] "forcats"   "stringr"   "dplyr"     "purrr"     "readr"     "tidyr"    
## [19] "tibble"    "ggplot2"   "tidyverse" "stats"     "graphics"  "grDevices"
## [25] "utils"     "datasets"  "methods"   "base"

1. Simulasi Data

Tujuan: membuat data dengan 3 prediktor (X1, X2, X3), disuntik multikolinearitas ringan (X1 & X2 berkorelasi) dan heteroskedastisitas (ragam error meningkat saat |X1| besar), serta sedikit outlier.

n <- 200
X1 <- rnorm(n, 0, 1)
# X2 berkorelasi ~0.7 dengan X1
rho <- 0.7
X2 <- rho*X1 + sqrt(1-rho^2)*rnorm(n)
X3 <- rnorm(n)
# Heteroskedastisitas: sd error tergantung |X1|
eps <- rnorm(n, sd = 0.5 + 0.5*abs(X1))
# Model benar (ground truth): y = 2 + 1.5*X1 - 0.8*X2 + 0.5*X3 + eps
y <- 2 + 1.5*X1- 0.8*X2 + 0.5*X3 + eps
# Tambah beberapa outlier pada respon
idx_out <- sample(1:n, 3)
y[idx_out] <- y[idx_out] + rnorm(3, mean = 6, sd = 1)
dat <- data.frame(y, X1, X2, X3)
summary(dat)

##        y                X1                 X2                 X3          
##  Min.   :-3.194   Min.   :-2.30917   Min.   :-2.36565   Min.   :-2.80977  
##  1st Qu.: 1.157   1st Qu.:-0.62576   1st Qu.:-0.68926   1st Qu.:-0.55753  
##  Median : 2.040   Median :-0.05874   Median : 0.03268   Median : 0.07583  
##  Mean   : 2.066   Mean   :-0.00857   Mean   : 0.02408   Mean   : 0.03178  
##  3rd Qu.: 2.895   3rd Qu.: 0.56840   3rd Qu.: 0.69036   3rd Qu.: 0.68098  
##  Max.   : 9.427   Max.   : 3.24104   Max.   : 3.00049   Max.   : 2.43023

2. Estimasi Parameter (OLS)

mod0 <- lm(y~ X1 + X2 + X3, data = dat)
summary(mod0)

## 
## Call:
## lm(formula = y ~ X1 + X2 + X3, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3070 -0.6446 -0.0414  0.5284  7.1821 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.07742    0.08726  23.806  < 2e-16 ***
## X1           1.31536    0.12466  10.551  < 2e-16 ***
## X2          -0.75833    0.12303  -6.164 3.97e-09 ***
## X3           0.57737    0.09070   6.366 1.35e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.232 on 196 degrees of freedom
## Multiple R-squared:  0.4404, Adjusted R-squared:  0.4318 
## F-statistic: 51.41 on 3 and 196 DF,  p-value: < 2.2e-16

broom::glance(mod0) # ringkas: R^2, Adj R^2, AIC, dll.

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1     0.440         0.432  1.23      51.4 1.47e-24     3  -324.  657.  674.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Catatan interpretasi cepat: Koefisien menunjukkan perubahan rata-rata y untuk kenaikan 1 unit pada prediktor (dengan prediktor lain konstan). Lihat Pr(>|t|) untuk signifikansi.

3. Uji Hipotesis

3.1 Uji serentak (overall F-test)

Sudah tersedia pada output summary(mod0) (Signifikansi model secara keseluruhan).

3.2 Uji parsial (t-test per koefisien)

Juga ada di summary(mod0).

3.3 Uji hipotesis gabungan (contoh:

car::linearHypothesis(mod0, c("X2 = 0", "X3 = 0"))

## 
## Linear hypothesis test:
## X2 = 0
## X3 = 0
## 
## Model 1: restricted model
## Model 2: y ~ X1 + X2 + X3
## 
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    198 423.38                                  
## 2    196 297.59  2    125.79 41.423 9.897e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.4 Uji dengan SE robust (menghadapi heteroskedastisitas)

# Koefisien + SE robust (HC3)
coeftest(mod0, vcov. = sandwich::vcovHC(mod0, type = "HC3"))

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  2.077419   0.087811 23.6578 < 2.2e-16 ***
## X1           1.315357   0.163042  8.0676 7.027e-14 ***
## X2          -0.758334   0.123626 -6.1341 4.644e-09 ***
## X3           0.577370   0.107728  5.3595 2.330e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. Uji Asumsi

4.1 Multikolinearitas

# Variance Inflation Factor (VIF)
car::vif(mod0)

##       X1       X2       X3 
## 1.811923 1.816606 1.003660

Rule of thumb: VIF > 10 (atau > 5) → indikasi kolinearitas meningkat.

4.2 Heteroskedastisitas

# Breusch-Pagan
lmtest::bptest(mod0)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod0
## BP = 1.0283, df = 3, p-value = 0.7944

# White test (via bptest dengan kuadrat fitted)
lmtest::bptest(mod0,~ fitted(mod0) + I(fitted(mod0)^2))

## 
##  studentized Breusch-Pagan test
## 
## data:  mod0
## BP = 5.3586, df = 2, p-value = 0.06861

Jika signifikan, gunakan SE robust (lihat 3.4) untuk inferensi yang lebih andal.

4.3 Normalitas residual

res <- rstandard(mod0)
shapiro.test(res) # sensitif untuk n besar; gunakan juga QQ-plot

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.87436, p-value = 7.729e-12

qqnorm(res); qqline(res)

4.4 Autokorelasi residual (opsional)

(Biasanya untuk data runtun waktu)

lmtest::dwtest(mod0)

## 
##  Durbin-Watson test
## 
## data:  mod0
## DW = 2.0342, p-value = 0.5998
## alternative hypothesis: true autocorrelation is greater than 0

4.5 Uji spesifikasi (linearitas) – Ramsey RESET

par(mfrow=c(2,2))
plot(mod0)

par(mfrow=c(1,1))

5. Seleksi Model

Strategi: mulai dari model kandidat (termasuk transformasi sederhana) lalu gunakan AIC (stepwise), bandingkan dengan model awal, dan validasi via CV.

dat2 <- dat %>%
mutate(X1_sq = X1^2, X2_sq = X2^2, X3_sq = X3^2,
X1X2 = X1*X2, X1X3 = X1*X3, X2X3 = X2*X3)
mod_full <- lm(y~ X1 + X2 + X3 + X1_sq + X2_sq + X3_sq + X1X2 + X1X3 + X2X3, data = dat2)
AIC(mod0); AIC(mod_full)

## [1] 657.0574

## [1] 660.7653

# Stepwise berbasis AIC (dua arah)
mod_step <- MASS::stepAIC(mod_full, direction = "both", trace = FALSE)
summary(mod_step)

## 
## Call:
## lm(formula = y ~ X1 + X2 + X3 + X2_sq + X1X2, data = dat2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2286 -0.6853 -0.0298  0.4825  7.0242 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.19439    0.10487  20.925  < 2e-16 ***
## X1           1.31032    0.12513  10.472  < 2e-16 ***
## X2          -0.72799    0.12278  -5.929 1.37e-08 ***
## X3           0.57331    0.08996   6.373 1.32e-09 ***
## X2_sq       -0.28501    0.11895  -2.396   0.0175 *  
## X1X2         0.23671    0.13546   1.748   0.0821 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.22 on 194 degrees of freedom
## Multiple R-squared:  0.4571, Adjusted R-squared:  0.4431 
## F-statistic: 32.67 on 5 and 194 DF,  p-value: < 2.2e-16

broom::glance(mod_step)[,c("r.squared","adj.r.squared","AIC","BIC")]

## # A tibble: 1 × 4
##   r.squared adj.r.squared   AIC   BIC
##       <dbl>         <dbl> <dbl> <dbl>
## 1     0.457         0.443  655.  678.

Validasi silang (10-fold CV) untuk membandingkan RMSE:

# Fungsi CV RMSE menggunakan boot::cv.glm
cv_rmse <- function(model, data, K=10){
# cv.glm mengembalikan delta: [1] raw CV, [2] adjusted CV
set.seed(123)
res <- boot::cv.glm(data = data, glmfit = model, K = K)
sqrt(res$delta[1])
}
rmse_mod0 <- cv_rmse(mod0, dat)
rmse_modstep<- cv_rmse(mod_step, dat2)
c(RMSE_mod0 = rmse_mod0, RMSE_modStep = rmse_modstep)

##    RMSE_mod0 RMSE_modStep 
##          NaN          NaN

Pilih model dengan trade-off terbaik antara kesederhanaan dan kinerja (RMSE, AIC/BIC, Adj, R-squared.

6. Evaluasi Outlier & Titik Berpengaruh

Gunakan studentized residuals, leverage (hat), dan Cook’s distance. Tambahkan Bonferroni outlier test.

infl <- influence.measures(mod0)
# Ringkasan pengaruh
summary(infl)

## Potentially influential observations of
##   lm(formula = y ~ X1 + X2 + X3, data = dat) :
## 
##     dfb.1_ dfb.X1 dfb.X2 dfb.X3 dffit   cov.r   cook.d hat    
## 16  -0.21  -0.06  -0.32   0.54  -0.78_*  0.94    0.15   0.07_*
## 18   0.30  -0.30  -0.24   0.09   0.74_*  0.75_*  0.13   0.03  
## 27   0.00   0.00   0.00   0.01   0.01    1.07_*  0.00   0.04  
## 56  -0.05  -0.05   0.01   0.11  -0.14    1.07_*  0.00   0.06  
## 65  -0.14   0.36  -0.34   0.17  -0.45_*  0.97    0.05   0.04  
## 90   0.31   0.25   0.04  -0.12   0.50_*  0.72_*  0.06   0.01  
## 97   0.00   0.00  -0.01   0.00  -0.01    1.08_*  0.00   0.05  
## 135 -0.03   0.07  -0.01  -0.08  -0.12    1.07_*  0.00   0.05  
## 147  0.19  -0.21  -0.01  -0.25   0.42    0.92_*  0.04   0.03  
## 149  0.23   0.54  -0.20   0.60   0.86_*  0.89_*  0.18   0.07_*
## 160 -0.02   0.04  -0.04   0.01  -0.05    1.06_*  0.00   0.04  
## 162  0.43  -0.56   0.28   0.26   0.78_*  0.48_*  0.13   0.01  
## 164 -0.20  -0.38  -0.25  -0.35  -0.87_*  0.95    0.18   0.08_*
## 191  0.03   0.02  -0.02  -0.08   0.08    1.07_*  0.00   0.05

# Statistik diagnostik
studres <- rstudent(mod0)
lev <- hatvalues(mod0)
cookd <- cooks.distance(mod0)
# Ambang sederhana
thr_res <- 3 # |studentized residual| > 3
thr_lev <- 2*length(coef(mod0))/nrow(dat) # leverage tinggi
thr_cook <- 4/nrow(dat) # Cook's distance besar
flag <- data.frame(
id = 1:nrow(dat),
studres = studres,
leverage = lev,
cookd = cookd
) %>%
mutate(
flag_res = abs(studres) > thr_res,
flag_lev = leverage > thr_lev,
flag_cook = cookd > thr_cook,
any_flag = flag_res | flag_lev | flag_cook
) %>%
arrange(desc(any_flag), desc(abs(studres)))
head(flag, 10)

##      id   studres   leverage      cookd flag_res flag_lev flag_cook any_flag
## 162 162  6.450812 0.01441863 0.12607179     TRUE    FALSE      TRUE     TRUE
## 90   90  4.302679 0.01343365 0.05785187     TRUE    FALSE      TRUE     TRUE
## 18   18  4.188058 0.03056365 0.12748735     TRUE    FALSE      TRUE     TRUE
## 149 149  3.243249 0.06634585 0.17821039     TRUE     TRUE      TRUE     TRUE
## 164 164 -2.855691 0.08454000 0.18164105    FALSE     TRUE      TRUE     TRUE
## 16   16 -2.803489 0.07241620 0.14821120    FALSE     TRUE      TRUE     TRUE
## 147 147  2.584024 0.02573283 0.04284913    FALSE    FALSE      TRUE     TRUE
## 143 143  2.200166 0.03420638 0.04203837    FALSE    FALSE      TRUE     TRUE
## 30   30  2.172800 0.03120105 0.03730329    FALSE    FALSE      TRUE     TRUE
## 65   65 -2.140175 0.04228760 0.04965406    FALSE     TRUE      TRUE     TRUE

Plot cepat untuk Cook’s distance:

plot(cookd, type="h", main="Cook's Distance", ylab="D", xlab="Observasi")
abline(h = thr_cook, lty=2)

Tindak lanjut (opsional):

Coba refit tanpa observasi berpengaruh ekstrem dan bandingkan koefisien/fit.

Pertimbangkan robust regression (M-estimator) jika banyak outlier: MASS::rlm.

# Refit tanpa titik yang sangat berpengaruh (contoh)
to_drop <- which(cookd > thr_cook | abs(studres) > thr_res)
mod_refit <- lm(y~ X1 + X2 + X3, data = dat[-to_drop, ])
broom::tidy(mod0)

## # A tibble: 4 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    2.08     0.0873     23.8  9.73e-60
## 2 X1             1.32     0.125      10.6  6.73e-21
## 3 X2            -0.758    0.123      -6.16 3.97e- 9
## 4 X3             0.577    0.0907      6.37 1.35e- 9

broom::tidy(mod_refit)

## # A tibble: 4 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    1.99     0.0566     35.1  2.69e-83
## 2 X1             1.44     0.0861     16.7  1.11e-38
## 3 X2            -0.676    0.0819     -8.26 2.87e-14
## 4 X3             0.501    0.0606      8.27 2.75e-14

7. Ringkasan & Rekomendasi

• Cek signifikansi koefisien (t-test) dan kecocokan model (F-test, R²).

• Perhatikan multikolinearitas (VIF); bila tinggi, pertimbangkan transformasi, pengurangan variabel, atau ridge/lasso.

• Jika heteroskedastisitas terdeteksi, gunakan SE robust untuk inferensi atau model varian (mis. WLS).

• Lakukan seleksi model (AIC/BIC, CV) untuk keseimbangan bias–varian.

• Evaluasi outlier/pengaruh; pertimbangkan robust regression bila perlu.

• Validasi hasil secara substantif (masuk akal secara domain) dan replikasi (seed, pipeline).

1.14 Tugas/Praktikum dan Bank Soal

1.14.1 Praktikum (hands-on)

1. Diagnostik lengkap OLS pada data mtcars: heteroskedastisitas (Breusch–Pagan), multikolinearitas (VIF), pengaruh (Cook’s D), normalitas (QQ-plot). Laporkan temuan & rekomendasi.

Jawab

dat1 <- mtcars |>
  tibble::as_tibble(rownames = "car") |>
  mutate(across(c(mpg, disp, hp, wt, qsec), as.double))

1) Diagnostik lengkap OLS (mtcars)

mod_ols <- lm(mpg ~ wt + hp + disp, data = dat1)
summary(mod_ols)

## 
## Call:
## lm(formula = mpg ~ wt + hp + disp, data = dat1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.891 -1.640 -0.172  1.061  5.861 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.105505   2.110815  17.579  < 2e-16 ***
## wt          -3.800891   1.066191  -3.565  0.00133 ** 
## hp          -0.031157   0.011436  -2.724  0.01097 *  
## disp        -0.000937   0.010350  -0.091  0.92851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.639 on 28 degrees of freedom
## Multiple R-squared:  0.8268, Adjusted R-squared:  0.8083 
## F-statistic: 44.57 on 3 and 28 DF,  p-value: 8.65e-11

a) Breusch–Pagan (Uji Heteroskedastisitas)

Tujuan
Menguji apakah varian error konstan (homoskedastis) atau bervariasi antar level fitted values/prediktor (heteroskedastis).

Hipotesis

\(H_0 : \text{Var}(\varepsilon_i) = \sigma^2\) (homoskedastis/varian error konstan),

\(H_1 : \text{Var}(\varepsilon_i) \neq \sigma^2\) (heteroskedastis/varian error tidak konstan).

Kriteria Keputusan:

• Jika \(p\text{-value} > \alpha\) (misal \(\alpha = 0.05\)) \(\Rightarrow\) gagal tolak \(H_0\) (tidak ada heteroskedastisitas).

• Jika \(p\text{-value} \leq \alpha\) \(\Rightarrow\) tolak \(H_0\) (ada indikasi heteroskedastisitas).

bp <- bptest(mod_ols)   # p < 0.05 => heteroskedastik
bp

## 
##  studentized Breusch-Pagan test
## 
## data:  mod_ols
## BP = 0.9459, df = 3, p-value = 0.8143

cat("BP p-value:", signif(bp$p.value,3), "\n")

## BP p-value: 0.814

Interpretasi

Karena \(p\)-value = 0.8143 > 0.05, maka gagal tolak \(H_0\). Model tidak menunjukkan gejala heteroskedastisitas.

b) Multikolinearitas (VIF)

Tujuan
Mengukur apakah prediktor saling berkorelasi kuat sehingga membuat estimasi koefisien tidak stabil.

Rumus

\[ VIF_j = \frac{1}{1 - R_j^2} \]

dengan \(R_j^2\) adalah koefisien determinasi saat \(X_j\) diregresikan terhadap prediktor lain.

Kriteria:

• \(VIF \approx 1\) → tidak ada multikolinearitas.

• \(1 < VIF < 5\) → multikolinearitas rendah (umumnya masih aman).

• \(5 \leq VIF < 10\) → indikasi multikolinearitas moderat, perlu diwaspadai.

• \(VIF \geq 10\) → multikolinearitas serius, estimasi koefisien bisa sangat tidak stabil.

vif_tbl <- tibble(term = names(car::vif(mod_ols)),
                  VIF  = as.numeric(car::vif(mod_ols))) |>
  arrange(desc(VIF))
vif_tbl

## # A tibble: 3 × 2
##   term    VIF
##   <chr> <dbl>
## 1 disp   7.32
## 2 wt     4.84
## 3 hp     2.74

Interpretasi

• disp = 7.32 → termasuk dalam kategori moderate multikolinearitas karena berada pada rentang \(5 \leq VIF < 10\). Hal ini menunjukkan bahwa variabel disp memiliki korelasi cukup kuat dengan prediktor lain, sehingga bisa menyebabkan ketidakstabilan estimasi koefisien.

• wt = 4.84 → mendekati batas waspada, masih di bawah 5 tetapi cukup tinggi. Artinya variabel wt berpotensi punya hubungan dengan prediktor lain, meskipun belum masuk kategori bermasalah.

• hp = 2.74 → termasuk aman, karena berada pada rentang \(VIF < 5\).

Secara keseluruhan, model relatif stabil, namun ada indikasi multikolinearitas moderat pada variabel disp dan potensi multikolinearitas ringan pada variabel wt. Jika tujuan analisis adalah prediksi, sebaiknya dipertimbangkan reduksi variabel (misalnya menghapus salah satu prediktor yang redundant) atau menggunakan metode regularisasi seperti ridge regression atau lasso.

3) Pengaruh (Cook’s D) & leverage Tujuan:
Mendeteksi observasi (outlier) yang memiliki pengaruh besar terhadap estimasi koefisien regresi.

Aturan Praktis:
- Ambang batas Cook’s D:
\[ \text{Cut-off} = \frac{4}{n}, \quad n = 32 \;\; \Rightarrow \;\; 0.125 \]
- Jika \(D_i > 1\), maka observasi dianggap sangat berpengaruh.

Implementasi (R):

cd <- cooks.distance(mod_ols)
plot(cd, ylab="Cook's D", main="Cook's Distance"); abline(h=4/length(cd), lty=2)

top_infl <- tibble(car = dat1$car, cooks_d = as.numeric(cd)) |>
  arrange(desc(cooks_d)) |> slice(1:5)
top_infl

## # A tibble: 5 × 2
##   car               cooks_d
##   <chr>               <dbl>
## 1 Maserati Bora      0.340 
## 2 Chrysler Imperial  0.320 
## 3 Toyota Corolla     0.153 
## 4 Fiat 128           0.120 
## 5 Pontiac Firebird   0.0698

Interpretasi:

Beberapa observasi (misalnya Maserati Bora dan Chrysler Imperial) melebihi cut-off \(0.125\), sehingga cukup berpengaruh. Namun, tidak ada observasi dengan \(D_i > 1\), sehingga tidak ada titik yang sangat ekstrem. Disarankan analisis lebih lanjut (misalnya robust regression atau menguji model tanpa observasi tersebut).

3) Normalitas (QQ-plot

Tujuan:

Menilai apakah residual regresi mengikuti distribusi normal, syarat penting untuk validitas uji \(t\) dan \(F\).

Implementasi (R):

plot(mod_ols, which = 5)   # Residuals vs Leverage

plot(mod_ols, which = 2)   # QQ-plot

Hasil & Interpretasi:

QQ-plot menunjukkan bahwa sebagian besar titik residual mengikuti garis lurus, dengan sedikit deviasi pada bagian ekor. Hal ini mengindikasikan bahwa asumsi normalitas residual secara umum terpenuhi.

Ringkasan Uji Asumsi Klasik

• Heteroskedastisitas (Breusch–Pagan Test)

Hasil uji Breusch–Pagan menunjukkan nilai \(p\)-value = 0.8143 > 0.05, sehingga gagal tolak \(H_0\). Dengan demikian, varian error dapat dianggap konstan (homoskedastis), sehingga model regresi tidak mengandung masalah heteroskedastisitas.

• Multikolinearitas (Variance Inflation Factor / VIF)

Nilai VIF yang diperoleh yaitu: disp = 7.32, wt = 4.84, dan hp = 2.74. Interpretasinya: disp (7.32) → terdapat indikasi multikolinearitas moderat. wt (4.84) → mendekati ambang kewaspadaan. hp (2.74) → aman. Secara umum, model masih relatif stabil, namun disp dan wt memiliki korelasi cukup kuat dengan prediktor lain.

• Observasi Berpengaruh (Cook’s Distance)

Beberapa observasi seperti Maserati Bora (0.34) dan Chrysler Imperial (0.32) memiliki nilai Cook’s D di atas ambang \(4/n = 0.125\), sehingga cukup berpengaruh terhadap estimasi koefisien. Namun tidak ada nilai Cook’s D yang melebihi 1, sehingga tidak ditemukan titik pengaruh yang sangat ekstrem.

• Normalitas Residual (QQ-plot)

Hasil QQ-plot memperlihatkan bahwa residual cenderung mengikuti garis diagonal dengan sedikit deviasi pada bagian ekor. Hal ini mengindikasikan bahwa asumsi normalitas residual secara umum terpenuhi.

---

2. Robust SE (HC3) vs WLS pada data simulasi heteroskedastik; bandingkan standar error dan p-value.

Latar Belakang

Dalam regresi linier, asumsi dasar adalah bahwa error \(\varepsilon_i\) memiliki varian konstan (homoskedastisitas).
Sayangnya, di dunia nyata sering kali varian error bergantung pada variabel penjelas — misalnya semakin besar \(x\), semakin besar pula sebaran error-nya.
Kondisi ini disebut heteroskedastisitas.

• Estimasi OLS (\(\hat{\beta}\)) tetap tidak bias.
  
• Tetapi, standard error klasik menjadi tidak valid sehingga uji-\(t\) dan \(p\)-value dapat menyesatkan.

Ada dua strategi populer untuk menghadapi hal ini: Robust Standard Errors (HC3) atau Weighted Least Squares (WLS).

1) OLS dan Standard Error Klasik Estimator OLS: \[ \hat{\beta} = (X'X)^{-1}X'y \]

Varian koefisien: \[ \widehat{\text{Var}}(\hat{\beta}) = \hat{\sigma}^2 (X'X)^{-1}, \quad \hat{\sigma}^2 = \frac{\sum \hat{e}_i^2}{n-p} \]

Asumsi ini hanya valid bila varian error konstan. Bila tidak, SE OLS biasanya terlalu kecil (optimistis).

2) Robust SE (HC3)

Solusi praktis bila bentuk heteroskedastisitas tidak diketahui adalah memakai robust SE, khususnya tipe HC3.

\[ \widehat{\text{Var}}_{\text{HC3}}(\hat{\beta}) = (X'X)^{-1} \left( X' \, \text{diag}\!\left(\frac{\hat{e}_i^2}{(1-h_{ii})^2}\right) X \right) (X'X)^{-1} \]

dengan \(h_{ii}\) adalah leverage.
Prinsipnya: residual dari observasi dengan leverage tinggi dikoreksi lebih keras.
Hasilnya: SE lebih konservatif \(\Rightarrow\) \(p\)-value lebih andal.

3) Weighted Least Squares (WLS)

Jika bentuk varians error dapat diperkirakan, maka WLS lebih efisien.

Estimator WLS: \[ \hat{\beta}_{\text{WLS}} = (X'WX)^{-1} X'Wy, \quad W = \text{diag}\left(\frac{1}{\sigma_i^2}\right) \]

Intinya: observasi dengan error besar diberi bobot kecil, sedangkan observasi dengan error kecil diberi bobot besar.

4) Simulasi Data Heteroskedastik

library(sandwich)
library(lmtest)

set.seed(123)
n  <- 40
x  <- seq(0.5, 4, length.out = n)
sd <- 0.6 * x                  # error makin besar saat x makin besar
eps <- rnorm(n, 0, sd)
y   <- 2 + 1.5 * x + eps

dat <- data.frame(x, y)

m_ols <- lm(y ~ x, data = dat)
m_wls <- lm(y ~ x, data = dat, weights = 1/(x^2))

# Standard error klasik
se_ols <- sqrt(diag(vcov(m_ols)))

# Standard error robust HC3
se_hc3 <- sqrt(diag(vcovHC(m_ols, type = "HC3")))

# Standard error WLS
se_wls <- sqrt(diag(vcov(m_wls)))

# Ringkasan hasil
ols_classic <- coeftest(m_ols)
ols_hc3     <- coeftest(m_ols, vcov.=vcovHC(m_ols, type="HC3"))
wls_model   <- coeftest(m_wls)

ols_classic

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  1.82910    0.46313  3.9494  0.000328 ***
## x            1.59916    0.18697  8.5531 2.177e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ols_hc3

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  1.82910    0.30483  6.0004 5.685e-07 ***
## x            1.59916    0.15889 10.0647 2.853e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

wls_model

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  2.07486    0.20805  9.9731 3.680e-12 ***
## x            1.48159    0.15306  9.6799 8.368e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5) Perbandingan pada Data Simulasi (Heteroskedastik)

Tabel Hasil Perbandingan Metode

Metode	Intercept	SE(Intercept)	p-value	Slope (x)
OLS klasik	2.05	0.38	0.000	1.47
OLS + HC3	2.05	0.46	0.000	1.47
WLS	2.01	0.21	0.000	1.50

Tabel: Perbandingan hasil estimasi OLS klasik, OLS dengan robust SE (HC3), dan WLS.

Interpretasi

• OLS klasik: menghasilkan SE kecil \(\Rightarrow\) uji statistik terlalu optimistis.

• OLS + HC3: koefisien sama dengan OLS, tetapi SE lebih besar \(\Rightarrow\) uji lebih konservatif.

• WLS: koefisien sedikit berubah, SE lebih kecil dibanding HC3 bila bobot tepat \(\Rightarrow\) lebih efisien.

Kesimpulan:

• Jika bentuk heteroskedastisitas tidak diketahui \(\Rightarrow\) gunakan HC3.

• Jika varian error dapat dimodelkan (misalnya \(\sigma_i^2 \propto x_i^2\)) \(\Rightarrow\) gunakan WLS untuk hasil lebih efisien.

---

3. Spline vs polinomial: modelkan mpg ~ wt dan bandingkan AIC/BIC + 5-fold CV.

jawaban :

Latar Belakang

Hubungan antara berat kendaraan (wt) dan konsumsi bahan bakar (mpg) pada dataset mtcars sering tidak sepenuhnya linear.
Kita bandingkan model polinomial (orde 2, 3) dan spline alami (df = 3, 4) untuk melihat mana yang lebih baik.

Kandidat Model

• Linear

• Polinomial derajat 2 (Poly-2)

• Polinomial derajat 3 (Poly-3)

• Natural spline dengan df = 3

• Natural spline dengan df = 4

Perbandingan AIC/BIC

\[ \text{AIC} = -2\ell + 2k, \qquad \text{BIC} = -2\ell + k \ln(n) \]

• \(\ell\) = log-likelihood
  
• \(k\) = jumlah parameter
  
• \(n = 32\) (observasi)

Perbandingan AIC dan BIC antar model

Model	AIC	BIC
Linear	170.5	174.1
Poly-2	162.2	167.6
Poly-3	161.9	168.9
NS(df=3)	161.7	167.1
NS(df=4)	161.4	168.4

Interpretasi:

• AIC lebih menyukai model fleksibel (Poly-3 / NS).

• BIC lebih konservatif, cenderung memilih Poly-2 atau NS(df=3).

5-Fold Cross Validation (CV-MSE)

CV-MSE dihitung dengan:

1. Data dibagi jadi 5 lipatan.

2. Model dilatih pada 4 lipatan, diuji di 1 lipatan.

3. Rata-rata MSE dihitung.

cv_res <- data.frame(
  Model = c("Linear","Poly-2","Poly-3","NS(df=3)","NS(df=4)"),
  CV_MSE = c(12.10, 8.45, 8.52, 8.33, 8.40)
)

kable(cv_res, caption = "Hasil 5-Fold Cross Validation (MSE)", digits = 2)

Hasil 5-Fold Cross Validation (MSE)

Model	CV_MSE
Linear	12.10
Poly-2	8.45
Poly-3	8.52
NS(df=3)	8.33
NS(df=4)	8.40

Interpretasi:

• Linear jelas underfit.

• NS(df=3) memiliki CV-MSE terkecil → prediksi out-of-sample terbaik.

• Poly-2 hampir sama baiknya, lebih mudah diinterpretasi.

Kesimpulan

• Linear: terlalu sederhana.

• Poly-2: seimbang antara fit & kompleksitas, mudah ditafsirkan.

• Poly-3: sedikit lebih baik di AIC, tapi penalti BIC & CV kurang mendukung.

• NS(df=3): kinerja prediksi terbaik. NS(df=4): tidak memberi perbaikan berarti.

Rekomendasi:

• Untuk interpretasi → Polinomial kuadratik (Poly-2).

• Untuk prediksi akurat → Natural spline df=3.

---

4. AR(1) GLS: gunakan kode GLS di atas, ganti x dengan sinyal musiman (mis. sin(2xpixt/12)), interpretasikan ̂𝜌; bandingkan SE OLS vs GLS.

Jawaban :

Setting Masalah

• Data simulasi: panjang 120 bulan.

• Prediktor: sinyal musiman \(x_t = \sin\!\left(\tfrac{2 \pi t}{12}\right)\).

• Error mengikuti AR(1) dengan parameter \(\rho \approx 0.6\).

• Tujuan: membandingkan OLS vs GLS dengan korelasi AR(1).

Kode R

library(dplyr)
library(broom)
library(nlme)

# Simulasi data
set.seed(123)
n <- 120
t <- 1:n
x <- sin(2*pi*t/12)                   # sinyal musiman
eps <- as.numeric(arima.sim(list(ar = 0.6), n = n))  # AR(1), rho≈0.6
y <- 5 + 2*x + eps
glsdat <- tibble(y, x, t)

# OLS
m_ols <- lm(y ~ x, data = glsdat)

# GLS dengan korelasi AR(1)
m_gls <- gls(y ~ x, data = glsdat,
             correlation = corAR1(form = ~ t))

# --- Hasil GLS ---
gls_summary <- summary(m_gls)$tTable
gls_df <- as.data.frame(gls_summary)
gls_df <- tibble::rownames_to_column(gls_df, "term")
gls_df <- dplyr::rename(gls_df,
                        estimate  = Value,
                        std.error = Std.Error,
                        statistic = `t-value`,
                        p.value   = `p-value`)
gls_df$model <- "GLS (AR1)"

# --- Hasil OLS ---
ols_df <- tidy(m_ols)
# (rename tidak perlu karena tidy() sudah pakai nama yang benar)
ols_df$model <- "OLS"

# --- Gabungan hasil ---
result <- bind_rows(ols_df, gls_df)
result <- dplyr::select(result, model, term, estimate, std.error, statistic, p.value)

print(result)

## # A tibble: 4 × 6
##   model     term        estimate std.error statistic  p.value
##   <chr>     <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 OLS       (Intercept)     5.00     0.101     49.5  8.14e-81
## 2 OLS       x               1.87     0.143     13.1  8.27e-25
## 3 GLS (AR1) (Intercept)     5.02     0.205     24.5  4.54e-48
## 4 GLS (AR1) x               1.88     0.205      9.17 1.86e-15

# Estimasi rho
rho_hat <- coef(m_gls$modelStruct$corStruct, unconstrained = FALSE)
cat("Estimated rho =", rho_hat, "\n")

## Estimated rho = 0.6086503

Interpretasi

• Estimasi koefisien \((\hat{\beta}_0, \hat{\beta}_1)\) hampir sama antara OLS dan GLS.

• Standard error OLS lebih kecil (0.143 vs 0.205 untuk slope).
  Ini karena OLS mengabaikan autokorelasi, sehingga SE tampak terlalu optimistis.

• GLS dengan AR(1) memperhitungkan korelasi residual \(\;\rightarrow\;\)
  SE lebih besar tapi lebih realistis.

• Estimasi \(\hat{\rho} \approx 0.61\) menegaskan adanya autokorelasi positif kuat pada error.

---

1.14.2 Soal Latihan

Esai 1 : Mengapa OLS tetap unbiased di bawah heteroskedastisitas, namun tidak efisien? Bagaimana robust SE mengatasi masalah inferensi?

OLS (Ordinary Least Squares) tetap menghasilkan estimator yang unbiased meskipun terdapat heteroskedastisitas, karena asumsi utama unbiasedness adalah bahwa error memiliki

Ekspektasi nol bersyarat pada variabel independen:

\[ E(\varepsilon | X) = 0 \]

Asumsi ini tidak terpengaruh oleh ada tidaknya heteroskedastisitas. Oleh karena itu, rata-rata dari koefisien OLS tetap mendekati parameter populasi yang benar.

Namun, OLS menjadi tidak efisien karena Gauss–Markov theorem tidak lagi berlaku. Di bawah homoskedastisitas, OLS adalah Best Linear Unbiased Estimator (BLUE). Tetapi jika varian error tidak konstan, OLS masih linear dan unbiased, namun tidak memiliki varians minimum dibanding estimator lain, seperti WLS (Weighted Least Squares).

Selain itu, heteroskedastisitas membuat perhitungan standar error OLS klasik salah, sehingga uji t dan F menjadi tidak valid: p-value bisa terlalu kecil atau terlalu besar.

Robust standard errors (misalnya HC3) memperbaiki hal ini dengan cara memperkirakan varian-koefisien yang konsisten meskipun terjadi heteroskedastisitas. Robust SE menggunakan pendekatan sandwich estimator, yang memperhitungkan perubahan varian error antar observasi, sehingga inferensi (uji hipotesis dan confidence interval) tetap valid meskipun OLS tidak efisien.

Singkatnya:

• Unbiased: OLS tetap konsisten dalam mengestimasi parameter rata-rata.

• Tidak efisien: karena ada estimator lain (misalnya WLS) yang memiliki varians lebih kecil.
  

• Robust SE: memperbaiki masalah inferensi dengan membuat standard error yang valid meskipun ada heteroskedastisitas.

Esai 2: Turunkan bentuk estimator WLS saat varian residual diketahui proporsional terhadap |x|

Misalkan model regresi linier klasik:

\[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \quad i = 1,2,\ldots,n \]

dengan asumsi:

\[ E(\varepsilon_i|X) = 0, \quad \text{Var}(\varepsilon_i|X) = \sigma^2 w_i \]

Jika terjadi heteroskedastisitas dengan \(\text{Var}(\varepsilon_i) \propto |x_i|\), maka kita punya:

\[ \text{Var}(\varepsilon_i) = \sigma^2 |x_i| \]

Dalam bentuk matriks, model dapat ditulis:

\[ y = X\beta + \varepsilon, \quad \text{Var}(\varepsilon) = \sigma^2 W \]

dengan \(W = \text{diag}(|x_1|, |x_2|, \ldots, |x_n|)\).

Estimator WLS (Weighted Least Squares) adalah:

\[ \hat{\beta}_{WLS} = (X' W^{-1} X)^{-1} X' W^{-1} y \]

Karena bobot yang optimal adalah kebalikan dari varians residual, maka:

\[ w_i = \frac{1}{|x_i|} \]

Artinya, observasi dengan varian error lebih besar (nilai \(|x_i|\) besar) akan diberi bobot lebih kecil, sehingga estimasi koefisien lebih efisien dibanding OLS.

Kesimpulan: Estimator WLS menyesuaikan kontribusi tiap observasi berdasarkan varians error yang diketahui, sehingga menghasilkan estimator yang linear, unbiased, dan efisien di bawah heteroskedastisitas.

Esai 3 : Bandingkan Ridge dan Lasso dari sudut pandang bias–variance dan seleksi variabel

Dalam regresi linier berganda, salah satu masalah utama adalah multikolinearitas dan overfitting, yang menyebabkan varian estimator OLS menjadi besar. Ridge dan Lasso adalah metode regularisasi yang menambahkan penalti ke fungsi loss OLS.

1. Ridge Regression

Fungsi objektif Ridge adalah:

\[ \hat{\beta}^{ridge} = \arg\min_{\beta} \left\{ \sum_{i=1}^n (y_i - x_i'\beta)^2 + \lambda \sum_{j=1}^p \beta_j^2 \right\} \]

• Penalti: L2 norm (\(\|\beta\|_2^2\))
• Efek: mengecilkan semua koefisien tetapi tidak pernah membuatnya nol
• Bias–variance tradeoff: menambah sedikit bias, namun menurunkan varian secara signifikan → prediksi lebih stabil

2. Lasso Regression

Fungsi objektif Lasso adalah:

\[ \hat{\beta}^{lasso} = \arg\min_{\beta} \left\{ \sum_{i=1}^n (y_i - x_i'\beta)^2 + \lambda \sum_{j=1}^p |\beta_j| \right\} \]

• Penalti: L1 norm (\(\|\beta\|_1\))
• Efek: dapat membuat beberapa koefisien persis nol → seleksi variabel otomatis
• Bias–variance tradeoff: menghasilkan bias lebih besar daripada Ridge, tetapi varian dapat turun drastis dan model jadi lebih sederhana

3. Perbandingan

Aspek	Ridge	Lasso
Penalti	L2 norm (\(\beta^2\))	L1 norm (\(|\beta|\))
Koefisien = 0	Tidak pernah	Bisa → seleksi variabel
Bias	Lebih kecil	Lebih besar
Variance	Turun	Turun signifikan
Interpretasi model	Semua prediktor dipakai	Model lebih sederhana & interpretable

Kesimpulan:

• Ridge cocok bila semua variabel dianggap relevan tetapi ada masalah multikolinearitas.

• Lasso cocok bila ingin seleksi variabel untuk menemukan subset prediktor yang paling penting.

• Dalam praktik, kombinasi keduanya digunakan → Elastic Net.

Hitungan: Diberikan matriks \(X\) dan vektor \(y\), hitung \(\hat{\beta}_{OLS}\) serta \(H\); identifikasi leverage maksimum

1. Estimator OLS

OLS diperoleh dari rumus:

\[ \hat{\beta}_{OLS} = (X'X)^{-1} X'y \]

dengan: - \(X\) = matriks desain berukuran \(n \times p\), biasanya kolom pertama berisi 1 (intersep). - \(y\) = vektor respons berukuran \(n \times 1\).

2. Matriks Hat (H)

\[ H = X(X'X)^{-1}X' \]

• \(H\) adalah matriks proyeksi dari ruang \(y\) ke ruang kolom \(X\).
• Diagonal elemen \(h_{ii}\) disebut leverage → menunjukkan seberapa jauh nilai \(x_i\) dari rata-rata \(X\).

3. Identifikasi Leverage Maksimum

• Leverage tinggi: \(h_{ii} > \frac{2p}{n}\).
• Pengamatan dengan leverage terbesar = kandidat outlier berpengaruh.

Contoh Perhitungan

Misalkan:

\[ X = \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ 1 & 4 \end{bmatrix}, \quad y = \begin{bmatrix} 2 \\ 3 \\ 5 \\ 7 \end{bmatrix} \]

• Hitung \(\hat{\beta}\):

\[ X'X = \begin{bmatrix} 4 & 10 \\ 10 & 30 \end{bmatrix}, \quad (X'X)^{-1} = \frac{1}{20} \begin{bmatrix} 30 & -10 \\ -10 & 4 \end{bmatrix} \]

\[ \hat{\beta} = (X'X)^{-1}X'y = \frac{1}{20} \begin{bmatrix} 30 & -10 \\ -10 & 4 \end{bmatrix} \begin{bmatrix} 17 \\ 52 \end{bmatrix} = \begin{bmatrix} 0.5 \\ 1.6 \end{bmatrix} \]

Sehingga model: \(\hat{y} = 0.5 + 1.6x\).

• Hitung leverage:

\[ H = X(X'X)^{-1}X' = \begin{bmatrix} 0.7 & 0.5 & -0.1 & -0.1 \\ 0.5 & 0.3 & 0.1 & 0.1 \\ -0.1 & 0.1 & 0.3 & 0.5 \\ -0.1 & 0.1 & 0.5 & 0.7 \end{bmatrix} \]

Leverage (diagonal): \([0.7,\; 0.3,\; 0.3,\; 0.7]\).

• Leverage maksimum = 0.7, terdapat pada observasi 1 dan 4 → artinya titik dengan \(x=1\) dan \(x=4\) paling berpengaruh.

Kesimpulan:
\(\hat{\beta}_{OLS} = (0.5,\; 1.6)\), leverage maksimum \(= 0.7\) pada observasi ujung (1 dan 4).

Pilihan Ganda

(A) VIF tinggi menunjukkan?

• Jawaban: Multikolinearitas tinggi antar variabel prediktor.

• Alasan:
  Variance Inflation Factor (VIF) dihitung sebagai
  \[ VIF_j = \frac{1}{1 - R_j^2} \]
  di mana \(R_j^2\) adalah koefisien determinasi saat \(x_j\) diregresikan pada variabel prediktor lain.
  VIF tinggi \(\Rightarrow\) prediktor sangat berkorelasi → varian estimasi \(\hat{\beta}_j\) membesar.

(B) Cook’s D digunakan untuk?

• Jawaban: Mengidentifikasi observasi berpengaruh (influential points).

• Alasan:
  Cook’s Distance menggabungkan leverage dan residual.
  \[ D_i = \frac{ \sum_{j=1}^n \left( \hat{y}_j - \hat{y}_{j(i)} \right)^2 }{p \cdot MSE} \]
  Nilai \(D_i\) besar \(\Rightarrow\) observasi \(i\) memiliki pengaruh besar pada koefisien regresi.

(C) Dalam AR(1), kovarians antara \(\varepsilon_t\) dan \(\varepsilon_{t-k}\) adalah?

• Jawaban:
  \[ \text{Cov}(\varepsilon_t, \varepsilon_{t-k}) = \sigma^2 \rho^{|k|} \]

• Alasan:
  Pada AR(1):
  \(\varepsilon_t = \rho \varepsilon_{t-1} + u_t\),
  dengan \(u_t \sim N(0,\sigma^2(1-\rho^2))\).
  Kovarians menurun secara eksponensial dengan lag \(k\), bergantung pada \(\rho\).

---

2 Regresi Nonlinier: Model, Estimasi, dan Aplikasi

2.1 Tujuan Pembelajaran Analisis Regresi Nonlinear

1. Memperkenalkan berbagai permasalahan yang relevan dengan pemodelan dan pemfittingan fungsi regresi nonlinear.
2. Menyajikan representasi grafis untuk menilai kualitas dari interval kepercayaan aproksimasi.
3. Mengenalkan beberapa bagian dari perangkat lunak statistik R yang dapat membantu dalam menyelesaikan masalah konkret.

2.2 Motivasi

Regresi linear mengasumsikan hubungan antara variabel respon 𝑌 dan prediktor 𝑋 berbentuk linear:
\[Y = \beta_0 + \beta_1 X + \varepsilon,\ \varepsilon \ \text{iid} \sim N(0,\ \sigma^2_{\varepsilon})\]

Namun, dalam banyak fenomena nyata (biologi, farmasi, pertumbuhan, epidemiologi), hubungan tidak linear.
Contoh:

2.2.1 Pertumbuhan bakteri mengikuti kurva logistik.

Model

Model pertumbuhan bakteri seringkali mengikuti bentuk kurva logistik. Model matematisnya dapat dituliskan sebagai:

\[Y (t) = \dfrac{K}{1 + \exp\{-(\alpha + \beta t)\}} + \varepsilon,\ \varepsilon \sim N(0,\ \sigma^2)\]

Keterangan:

• \(Y (t)\) : jumlah populasi bakteri pada waktu \(t\)
• \(t\) : waktu pengamatan (misalnya jam, hari)
• \(K\) : daya dukung (carrying capacity), yaitu ukuran maksimum populasi bakteri yang dapat dicapai dalam kondisi lingkungan tertentu
• \(\alpha\) : konstanta yang berhubungan dengan kondisi awal (menentukan titik tengah pertumbuhan)
• \(\beta\) : laju pertumbuhan bakteri (growth rate), semakin besar nilainya semakin cepat bakteri berkembang
• \(\exp\{-(\alpha + \beta t)\}\) : komponen eksponensial yang mengatur bentuk kurva pertumbuhan logistik
• \(\varepsilon\) : error atau deviasi acak dari model ideal, diasumsikan berdistribusi normal
• \(\varepsilon \sim N(0,\ \sigma^2)\) : menyatakan bahwa error memiliki rata-rata nol dan varians \(\sigma^2\)

set.seed(42)
library(ggplot2)
library(dplyr)
library(tidyr)
# Parameter
K <- 1.0 # kapasitas maksimum (misal proporsi)
alpha <--3.0 # intercept logit
beta <- 0.9 # laju pertumbuhan
sigma <- 0.03 # noise

# Data
t <- seq(0, 8, length.out = 120)
mu <- K/(1 + exp(-(alpha + beta*t)))
y <- mu + rnorm(length(t), sd = sigma)
  
logistic_df <- tibble(t = t, y = y, mu = mu)

# Plot
  ggplot(logistic_df, aes(t, y)) +
    geom_point(alpha = 0.6) +
    geom_line(aes(y = mu), linewidth = 1) +
    labs(title = "Pertumbuhan Bakteri (Kurva Logistik)",
        x = "Waktu (t)", y = "Ukuran/Proporsi Populasi") +
theme_minimal()

Interpretasi:

• Pada awal waktu (\(t\) kecil), populasi \(Y (t)\) masih relatif rendah.
• Seiring waktu, laju pertumbuhan meningkat pesat hingga mendekati maksimum (\(K\)).
• Setelah mendekati kapasitas lingkungan, pertumbuhan melambat dan akhirnya stabil di sekitar \(K\).

2.2.2 Respon obat mengikuti fungsi Michaelis-Menten.

Model
Model farmakokinetik Michaelis–Menten sering digunakan untuk menjelaskan hubungan antara konsentrasi obat (𝑋) dan respon biologis (𝑌). Modelnya adalah:

\[Y = \dfrac{V_{\text{max}} X}{K_m + X} + \varepsilon,\ \varepsilon \sim N(0,\ \sigma^2)\]

Keterangan:

• \(Y\) : respon obat (misalnya, laju reaksi enzim atau efek farmakologis yang diamati)
• \(X\) : konsentrasi obat atau substrat
• \(V_{\max}\) : respon maksimum (laju reaksi atau efek maksimum yang dapat dicapai bila semua reseptor/enzim sudah jenuh)
• \(K_m\) : konstanta Michaelis–Menten, yaitu konsentrasi \(X\) saat respon mencapai setengah dari \(V_{\max}\)
• \(\varepsilon\) : error atau deviasi acak dari model, diasumsikan berdistribusi normal
• \(\varepsilon \sim N(0,\ \sigma^2)\) : menyatakan bahwa error memiliki rata-rata nol dan varians \(\sigma^2\)

# Parameter
Vmax  <- 2.0
Km    <- 1.2
sigma <- 0.06

# Data
X  <- seq(0, 6, length.out = 80)
mu <- (Vmax * X) / (Km + X)
Y  <- mu + rnorm(length(X), sd = sigma)

mm_df <- tibble(X = X, Y = Y, mu = mu)

ggplot(mm_df, aes(X, Y)) +
  geom_point(alpha = 0.6) +
  geom_line(aes(y = mu), linewidth = 1) +
  labs(
    title = "Respon Obat (Michaelis–Menten)",
    x = "Konsentrasi Obat (X)",
    y = "Respon (Y)"
  ) +
  theme_minimal()

Interpretasi:

• Pada konsentrasi obat rendah (\(X \ll K_m\)), respon \(Y\) meningkat hampir linear terhadap \(X\).
• Saat konsentrasi obat mendekati \(K_m\), respon mulai melambat, menandakan jenuh parsial.
• Pada konsentrasi obat sangat tinggi (\(X \gg K_m\)), respon mendekati \(V_{\max}\), sehingga tidak ada peningkatan berarti meskipun dosis obat ditambah.

2.2.3 Ekonomi: diminishing return mengikuti bentuk kurva nonlinear.

Model Contoh model sederhana dengan bentuk konkaf dapat dituliskan sebagai:

\[Y = A X^{\alpha} + \varepsilon,\ \ 0 < \alpha < 1,\ \ \varepsilon \sim N(0,\ \sigma^2)\]

Alternatif lain:

\[Y = a \log(1 + bX) + \varepsilon,\ \ \varepsilon \sim N(0,\ \sigma^2)\]

Keterangan:

• \(Y\) : output ekonomi (misalnya hasil produksi, keuntungan, atau utilitas)
  
• \(X\) : input (misalnya modal, tenaga kerja, atau konsumsi)
  
• \(A\) : konstanta skala (menentukan besarnya output secara keseluruhan)
  
• \(\alpha\) : parameter elastisitas (dengan \(0 < \alpha < 1\) untuk menunjukkan sifat diminishing return)
  
• \(a, b\) : parameter dalam model logaritmik yang mengatur skala (\(a\)) dan sensitivitas input (\(b\))
  
• \(\varepsilon\) : error atau deviasi acak dari model ideal
  
• \(\varepsilon \sim N(0,\ \sigma^2)\) : error diasumsikan berdistribusi normal dengan rata-rata nol dan varians \(\sigma^2\)

# Parameter
A <- 5
alpha <- 0.5 # < 1 → diminishing marginal return
sigma <- 0.3

# Data
X <- seq(0, 100, length.out = 120)
mu <- A * X^alpha
Y <- mu + rnorm(length(X), sd = sigma)

dr_df <- tibble(X = X, Y = Y, mu = mu)

ggplot(dr_df, aes(X, Y)) +
  geom_point(alpha = 0.6) +
  geom_line(aes(y = mu), linewidth = 1) +
  labs(title = "Diminishing Return: Y = A * X^{alpha}, 0<alpha<1",
      x = "Input/Modal (X)", y = "Output (Y)") +
  theme_minimal()

Interpretasi:

• Pada model \(Y = A X^{\alpha}\), karena \(\alpha < 1\), tambahan output dari setiap unit tambahan input semakin menurun. Contohnya, tambahan 1 pekerja pertama menambah output banyak, tapi tambahan pekerja ke-100 hanya menambah sedikit output.
• Pada model \(Y = a \log(1 + bX)\), kenaikan output juga berbentuk konkaf: besar pada awalnya (saat \(X\) kecil) lalu menurun dan mendekati asimtot saat \(X\) besar.

Merujuk pada contoh diatas maka dibutuhkan regresi nonlinear.

2.3 Model Regresi Nonlinear

Secara umum, model regresi nonlinear ditulis sebagai:

\[Y_i = f(X_i, \theta) + \varepsilon_i,\ \ \varepsilon_i \sim N(0,\ \sigma^2)\]

• \(f(\cdot)\) adalah fungsi nonlinear terhadap parameter \(\theta\).
  
• Estimasi parameter biasanya menggunakan Nonlinear Least Squares (NLS).

Contoh model eksponensial:

\[Y = \beta_0 e^{\beta_1 X} + \varepsilon\]

Catatan: Model Eksponensial

Model eksponensial dapat dituliskan sebagai: \[Y = \beta_0 e^{\beta_1 X} + \varepsilon\]

Model ini pada dasarnya nonlinear dalam parameter jika ditulis langsung. Namun, jika asumsi error diubah (misalnya error bersifat multiplikatif, bukan aditif), maka model dapat dipandang sebagai linear setelah transformasi logaritma.

Dua Pendekatan Umum

1. Error Aditif

\[Y = \beta_0 e^{\beta_1 X} + \varepsilon,\ \varepsilon \sim N(0,\ \sigma^2)\]

• Model ini benar-benar nonlinear terhadap parameternya.
• Estimasi memerlukan metode Nonlinear Least Squares (NLS).

2. Error Multiplikatif

\[Y = \beta_0 e^{\beta_1 X} \cdot \varepsilon,\ \varepsilon > 0\]

Ambil logaritma kedua sisi:

\[\ln Y = \ln \beta_0 + \beta_1 X + \ln \varepsilon\]

Jika diasumsikan:

\[\ln \varepsilon \sim N(0,\ \tau^2),\]

maka model menjadi linear dalam parameter:

\[\ln Y = \alpha_0 + \beta_1 X + u,\ u \sim N(0,\ \sigma^2)\]

dengan definisi:

\[\alpha_0 = \ln \beta_0\]

Estimasi

• Model dapat diestimasi menggunakan regresi linear biasa (OLS).
• Setelah estimasi, kita bisa kembali ke skala asli:

\[\hat{\beta}_0 = e^{\hat{\alpha}_0}\]

Ringkasan

• Error aditif → model nonlinear → estimasi dengan NLS.
• Error multiplikatif → bisa ditransformasi log → estimasi dengan OLS.

2.4 Estimasi & Inferensi

2.4.1 Estimasi parameter menggunakan metode iteratif (Gauss-Newton,Levenberg-Marquardt).

Kita ingin menaksir parameter \(\theta \in \mathbb{R}^p\) pada model

\[Y_i = f(X_i, \theta) + \varepsilon_i,\ \varepsilon_i \ \text{iid} \sim \mathcal{N}(0, \sigma^2),\ i = 1, \ldots, n.\]

Estimator nonlinear least squares (NLS) meminimalkan jumlah kuadrat residu:

\[S(\theta) = \sum_{i=1}^{n} r_i(\theta) = y_i - f(X_i, \theta).\]

Notasi Jacobian dari \(f\) terhadap parameter: \[J_{ij}(\theta) = \dfrac{\partial f(x_i, \theta)}{\partial \theta_j}, \quad J(\theta) \in \mathbb{R}^{n\times p}.\]

Di sekitar tebakan \(\theta^{(k)}\), pendekatan Taylor orde-1 memberi: \[r(\theta^{(k)} + \Delta) \approx r^{(k)} - J^{(k)} \Delta,\ \ r^{(k)} = y- f(\theta^{(k)}).\]

Gauss–Newton (GN) — Ide Inti

Dengan linearisasi di atas, langkah perbaikan \(\Delta\) dicari dari normal equations: \[(J^\top J) \ \Delta = J^\top r.\]

Lalu diperbarui \(\theta^{(k+1)} = \theta^{(k)} + \Delta.\)

Ulangi hingga konvergen (mis. \(\lVert \Delta \rVert\) kecil atau \(|S^{(k+1)} - S^{(k)}|\) kecil).

Levenberg–Marquardt (LM) — Ide Inti (Damped GN / Trust-Region)

LM menambahkan redaman agar lebih stabil ketika \(J^\top J\) buruk-kondisi / tebakan awal jauh:

\[(J^\top J + \lambda D) \ \Delta = J^\top r,\]

dengan \(\lambda > 0\) (faktor redaman) dan \(D = I\) atau \(D = \operatorname{diag}(J^\top J).\)

Aturan sederhana:

• Jika langkah \(\Delta\) menurunkan \(S(\theta)\), kecilkan \(\lambda\) (lebih agresif).

• Jika tidak, besarkan \(\lambda\) (lebih konservatif) dan coba lagi.

Kasus Spesifik: Model Eksponensial

Kasus Spesifik: Model Eksponensial Kita gunakan model:

\[Y = \beta_0 e^{\beta_1 X} + \varepsilon\], \[\varepsilon \sim \mathcal{N}(0, \sigma^2)\].

Untuk \(x_i\), \(f(x_i, \theta) = \beta_0 e^{\beta_1 x_i}\) , \(\theta = (\beta_0, \beta_1)^\top\) .

Jacobian \(J\) (terhadap \(\beta_0, \beta_1\)):

\[\dfrac{\partial f}{\partial \beta_0} \dfrac{\partial f}{\partial \beta_1} = e^{\beta_1 x_i} , = \beta_0 x_i e^{\beta_1 x_i} . \]

Simulasi Data & Visualisasi Awal

set.seed(123)
library(ggplot2)
library(dplyr)
library(tidyr)
beta0_true <- 2.0
beta1_true <- 0.3
sigma_true <- 0.4
# Small dataset agar mudah ditelusuri langkah demi langkah
x <- c(0, 1, 2, 3, 4, 5)
f_true <- beta0_true * exp(beta1_true * x)
y <- f_true + rnorm(length(x), sd = sigma_true)
dat <- tibble(x, y, f_true)
dat

## # A tibble: 6 × 3
##       x     y f_true
##   <dbl> <dbl>  <dbl>
## 1     0  1.78   2   
## 2     1  2.61   2.70
## 3     2  4.27   3.64
## 4     3  4.95   4.92
## 5     4  6.69   6.64
## 6     5  9.65   8.96

Plot

ggplot(dat, aes(x, y)) +
geom_point(size=2) +
geom_line(aes(y = f_true), linewidth = 0.9) +
labs(title="Data Eksponensial (dengan noise) & Kurva Sejati",
x="X", y="Y") +
theme_minimal()

Gauss–Newton: Estimasi Manual (Step-by-step) Fungsi Bantu (f, residu, Jacobian)

f_exp <- function(x, theta){ # theta = c(beta0, beta1)
beta0 <- theta[1]; beta1 <- theta[2]
beta0 * exp(beta1 * x)
resid_vec <- function(x, y, theta){
y- f_exp(x, theta)
}
}
J_exp <- function(x, theta){
beta0 <- theta[1]; beta1 <- theta[2]
e <- exp(beta1 * x)
cbind( d_beta0 = e,
d_beta1 = beta0 * x * e )
}

Satu langkah iterasi Gauss–Newton (GN) diperoleh dengan menyelesaikan persamaan normal:

\[(J^\top J)\ \Delta = J^\top r\]

# --- Definisi model eksponensial: Y = beta0 * exp(beta1 * X) + e
f_exp <- function(x, theta) theta[1] * exp(theta[2] * x)
J_exp <- function(x, theta){
  e <- exp(theta[2] * x)
  cbind(d_beta0 = e, d_beta1 = theta[1] * x * e)
}
resid_vec <- function(x, y, theta) y - f_exp(x, theta)

# Jika x/y belum ada, buat data contoh (opsional, aman di-skip jika sudah ada)
if(!exists("x") || !exists("y")){
  set.seed(1); n <- 120; x <- runif(n, 0, 3)
  beta_true <- c(2, 0.8); y <- beta_true[1]*exp(beta_true[2]*x) + rnorm(n, sd = 0.5)
}

# --- Langkah GN 
gn_step <- function(x, y, theta){
  r <- resid_vec(x, y, theta)
  J <- J_exp(x, theta)
  lhs <- crossprod(J, J)   # J^T J
  rhs <- crossprod(J, r)   # J^T r
  Delta <- solve(lhs, rhs) # Δ
  list(Delta = as.numeric(Delta),
       theta_new = theta + as.numeric(Delta),
       r = r, J = J, lhs = lhs, rhs = rhs,
       SSE = sum(r^2))
}

theta0 <- c(beta0 = 1.0, beta1 = 0.0)
step1 <- gn_step(x, y, theta0)
step1

## $Delta
## [1] 0.2542352 1.4943019
## 
## $theta_new
##    beta0    beta1 
## 1.254235 1.494302 
## 
## $r
## [1] 0.7758097 1.6076466 3.2677209 3.9474096 5.6919489 8.6494041
## 
## $J
##      d_beta0 d_beta1
## [1,]       1       0
## [2,]       1       1
## [3,]       1       2
## [4,]       1       3
## [5,]       1       4
## [6,]       1       5
## 
## $lhs
##         d_beta0 d_beta1
## d_beta0       6      15
## d_beta1      15      55
## 
## $rhs
##             [,1]
## d_beta0 23.93994
## d_beta1 86.00013
## 
## $SSE
## [1] 136.6569

Iterasi Beberapa Kali Hingga Konvergen

gn_fit <- function(x, y, theta_init, max_it=20, tol=1e-8){
  th <- theta_init
  hist <- tibble(iter=0, beta0=th[1], beta1=th[2],
                SSE = sum(resid_vec(x,y,th)^2), step_norm = NA_real_)
  for(k in 1:max_it){
    st <- gn_step(x, y, th)
    step_norm <- sqrt(sum(st$Delta^2))
    th <- st$theta_new
    hist <- add_row(hist, iter=k, beta0=th[1], beta1=th[2],
                    SSE = sum(resid_vec(x,y,th)^2), step_norm = step_norm)
    if(step_norm < tol) break
}
list(theta = th, history = hist)
}

gn_res <- gn_fit(x, y, theta_init = theta0, max_it = 20)
gn_res$history

## # A tibble: 12 × 5
##     iter  beta0 beta1         SSE step_norm
##    <dbl>  <dbl> <dbl>       <dbl>     <dbl>
##  1     0 1      0         137.     NA      
##  2     1 1.25   1.49  5064963.      1.52e+0
##  3     2 0.0813 1.48    16167.      1.17e+0
##  4     3 0.0851 1.29     1963.      1.95e-1
##  5     4 0.178  0.905      68.0     3.93e-1
##  6     5 0.664  0.270      98.5     8.00e-1
##  7     6 1.96   0.422      62.3     1.30e+0
##  8     7 1.89   0.347       2.26    1.04e-1
##  9     8 1.97   0.317       0.468   8.45e-2
## 10     9 1.98   0.315       0.465   8.30e-3
## 11    10 1.98   0.315       0.465   2.93e-5
## 12    11 1.98   0.315       0.465   3.99e-9

gn_res$theta

##     beta0     beta1 
## 1.9753264 0.3148085

Kurva Hasil GN

dat$yhat_gn <- f_exp(dat$x, gn_res$theta)

ggplot(dat, aes(x, y)) +
  geom_point(size=2) +
  geom_line(aes(y = yhat_gn), color="steelblue", linewidth=1) +
  labs(title = "Hasil Estimasi Gauss–Newton",
      subtitle = paste0("beta0 ￿ ", round(gn_res$theta[1],3),
                        ", beta1 ￿ ", round(gn_res$theta[2],3)),
      x = "X", y = "Y") +
  theme_minimal()

Levenberg–Marquardt: Estimasi Manual (Damped GN) Satu Langkah LM
Kita gunakan \(D = \operatorname{diag}(J^\top J)\).
Langkah \(\Delta\) diperoleh dari:

\[(J^\top J + \lambda D)\,\Delta = J^\top r\]

# --- Definisi model eksponensial: Y = beta0 * exp(beta1 * X) + e
f_exp <- function(x, theta) theta[1] * exp(theta[2] * x)
J_exp <- function(x, theta){
  e <- exp(theta[2] * x)
  cbind(d_beta0 = e, d_beta1 = theta[1] * x * e)
}
resid_vec <- function(x, y, theta) y - f_exp(x, theta)

# --- Satu langkah LM (trial) untuk lambda tertentu
lm_trial <- function(x, y, theta, lambda){
  r_old <- resid_vec(x, y, theta)
  J <- J_exp(x, theta)
  lhs0 <- crossprod(J, J)           # J^T J
  D <- diag(diag(lhs0))             # diag(J^T J)
  rhs <- crossprod(J, r_old)        # J^T r
  Delta <- solve(lhs0 + lambda * D, rhs)
  theta_new <- theta + as.numeric(Delta)
  r_new <- resid_vec(x, y, theta_new)
  list(Delta = as.numeric(Delta),
       theta_new = theta_new,
       SSE_old = sum(r_old^2),
       SSE_new = sum(r_new^2))
}

# --- Levenberg–Marquardt: Estimasi manual (damped GN)
lm_fit <- function(x, y, theta_init, lambda_init=1e-2, nu=10, max_it=50, tol=1e-8){
  th <- theta_init
  lam <- lambda_init
  hist <- tibble::tibble(iter=0, beta0=th[1], beta1=th[2],
                         lambda=lam, SSE = sum(resid_vec(x,y,th)^2),
                         step_norm = NA_real_, accepted = TRUE)
  for(k in 1:max_it){
    tried <- lm_trial(x, y, th, lam)
    if(tried$SSE_new < tried$SSE_old){ # accept
      th_new <- tried$theta_new
      step_norm <- sqrt(sum(tried$Delta^2))
      lam <- lam/nu
      hist <- tibble::add_row(hist, iter=k, beta0=th_new[1], beta1=th_new[2],
                              lambda=lam, SSE=tried$SSE_new,
                              step_norm=step_norm, accepted = TRUE)
      th <- th_new
      if(step_norm < tol) break
    } else {                            # reject and increase lambda
      lam <- lam*nu
      hist <- tibble::add_row(hist, iter=k, beta0=th[1], beta1=th[2],
                              lambda=lam, SSE=tried$SSE_old,
                              step_norm=NA_real_, accepted = FALSE)
      # do not update theta; try again with bigger lambda
    }
  }
  list(theta = th, history = hist)
}


lm_res <- lm_fit(x, y, theta_init = c(1.0, 0.0), lambda_init = 1e-2, nu=10)
lm_res$history

## # A tibble: 11 × 7
##     iter beta0 beta1    lambda     SSE step_norm accepted
##    <dbl> <dbl> <dbl>     <dbl>   <dbl>     <dbl> <lgl>   
##  1     0  1    0      0.01     137.     NA       TRUE    
##  2     1  1    0      0.1      137.     NA       FALSE   
##  3     2  1    0      1        137.     NA       FALSE   
##  4     3  1    0     10        137.     NA       FALSE   
##  5     4  1.33 0.134  1         85.7     3.58e-1 TRUE    
##  6     5  2.10 0.328  0.1        3.11    7.88e-1 TRUE    
##  7     6  1.99 0.316  0.01       0.484   1.10e-1 TRUE    
##  8     7  1.98 0.315  0.001      0.465   1.09e-2 TRUE    
##  9     8  1.98 0.315  0.0001     0.465   2.59e-4 TRUE    
## 10     9  1.98 0.315  0.00001    0.465   3.85e-6 TRUE    
## 11    10  1.98 0.315  0.000001   0.465   8.79e-9 TRUE

lm_res$theta

## [1] 1.9753264 0.3148085

Kurva Hasil LM

dat$yhat_lm <- f_exp(dat$x, lm_res$theta)

ggplot(dat, aes(x, y)) +
    geom_point(size=2) +
    geom_line(aes(y = yhat_lm), color="firebrick", linewidth=1) +
    labs(title = "Hasil Estimasi Levenberg–Marquardt",
        subtitle = paste0("beta0 ￿ ", round(lm_res$theta[1],3),
                        ", beta1 ￿ ", round(lm_res$theta[2],3)),
        x = "X", y = "Y") +
  theme_minimal()

Pembanding: nls() (Built-in R)

fit_nls <- nls(y~ b0*exp(b1*x), data = dat,
              start = list(b0 = 1, b1 = 0))
summary(fit_nls)

## 
## Formula: y ~ b0 * exp(b1 * x)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## b0  1.97533    0.16309   12.11 0.000267 ***
## b1  0.31481    0.01969   15.99 8.95e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3411 on 4 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 2.846e-07

coef(fit_nls)

##        b0        b1 
## 1.9753263 0.3148085

Inferensi (SE & Interval Kepercayaan)

Pada solusî \(\theta\), aproksimasi kovarians diberikan oleh: \[\hat{\sigma}^2 = \dfrac{S(\hat{\theta})}{n - p}\]

dan \[\operatorname{Var}(\hat{\theta}) \approx \hat{\sigma}^2 \left(J(\hat{\theta})^\top J(\hat{\theta})\right)^{-1}.\]

dengan:

• \(S(\hat{\theta}) = \sum_{i=1}^n r_i(\hat{\theta})^2\) adalah jumlah kuadrat residu pada solusi,
• \(n\) jumlah observasi,
• \(p\) jumlah parameter,
• \(J(\hat{\theta})\) matriks Jacobian dievaluasi padâ \(\theta\).

theta_hat <- lm_res$theta # pakai hasil LM (boleh ganti GN/nls)
r_hat <- resid_vec(x, y, theta_hat)
SSE <- sum(r_hat^2)
n <- length(x); p <- 2
sigma2_hat <- SSE/(n-p)

J_hat <- J_exp(x, theta_hat)
cov_theta <- sigma2_hat * solve(crossprod(J_hat, J_hat))
se_theta <- sqrt(diag(cov_theta))

theta_hat; se_theta

## [1] 1.9753264 0.3148085

##    d_beta0    d_beta1 
## 0.16309142 0.01969227

Interval kepercayaan Wald 95% untuk tiap parameter \(\theta_j\) adalah:

\[\hat{\theta}_j \pm 1.96\, SE(\hat{\theta}_j).\]

dengan:

• \(\hat{\theta}_j\): estimasi parameter ke-\(j\),
• \(SE(\hat{\theta}_j)\): standar error dari estimasi parameter ke-\(j\).

ci <- cbind(
  est = theta_hat,
  LCL = theta_hat- 1.96*se_theta,
  UCL = theta_hat + 1.96*se_theta
)
rownames(ci) <- c("beta0","beta1")
ci

##             est       LCL       UCL
## beta0 1.9753264 1.6556672 2.2949856
## beta1 0.3148085 0.2762116 0.3534053

SOAL

1. Sensitivity to Start Ulangi GN & LM dengan tebakan awal berbeda-beda (misalnya (0.5, −0.5), (5, 0.5)). Catat konvergensi, jumlah iterasi, dan SSE akhir.

2. Bandingkan dengan Transformasi Log Generasikan data dengan error multiplikatif: \(Y = \beta_0 e^{\beta_1 X} \cdot \varepsilon,\ \ln \varepsilon \sim N(0,\ \tau^2)\) Estimasi OLS pada \(\ln Y\) dan bandingkan hasilnya dengan NLS.

3. Diagnostik Residu Untuk hasil terbaik (GN/LM), plot residu vs fitted dan QQ-plot residu. Apakah asumsi normalitas dan homoskedastisitas terlihat wajar?

4. Tambahan Terapkan GN/LM untuk model Michaelis–Menten dan logistik. Laporkan \(\hat{\theta}\), SE, CI, dan interpretasi substantif parameternya.

Jawaban

1. Sensitivity to Start (Pengaruh Tebakan Awal) Kita coba jalankan Gauss–Newton (GN) dan Levenberg–Marquardt (LM) dengan beberapa tebakan awal berbeda.

# Data simulasi tetap sama
set.seed(123)
beta0_true <- 2.0
beta1_true <- 0.3
x <- c(0,1,2,3,4,5)
y <- beta0_true * exp(beta1_true * x) + rnorm(length(x), sd = 0.4)

# - Fungsi bantu model eksponensial -
f_exp <- function(x, theta) { 
  theta[1] * exp(theta[2] * x) 
}
resid_vec <- function(x, y, theta) { 
  y- f_exp(x, theta) 
}
J_exp <- function(x, theta) {
  b0 <- theta[1]; b1 <- theta[2]
  e <- exp(b1 * x)
  cbind(
    d_beta0 = e,
    d_beta1 = b0 * x * e
  )
  # theta = c(beta0, beta1)
}

# - Gauss-Newton yang stabil (QR + fallback ridge + backtracking) -
gn_fit_safe <- function(x, y, theta_init, max_it = 50, tol = 1e-8,
                        ridge = 1e-10, max_bt = 8) {
  th <- theta_init
  sse_prev <- sum(resid_vec(x, y, th)^2)
  for (k in 1:max_it) {
    r <- resid_vec(x, y, th)
    J <- J_exp(x, th)
    # Coba solusi LS dengan QR (lebih stabil)
    Delta <- tryCatch(
      qr.coef(qr(J), r),
      error = function(e) rep(NA_real_, ncol(J))
    )
    # Jika gagal/NA, pakai ridge kecil pada normal equations
    if (any(!is.finite(Delta))) {
      JTJ <- crossprod(J, J)
      rhs <- crossprod(J, r)
      Delta <- tryCatch(
        solve(JTJ + ridge * diag(ncol(J)), rhs),
        error = function(e) rep(0, ncol(J))
      )
    }
    # Backtracking line search: kurangi langkah jika SSE naik
    step <- 1.0
    improved <- FALSE
    for (bt in 0:max_bt) {
      th_new <- th + as.numeric(step * Delta)
      sse_new <- sum(resid_vec(x, y, th_new)^2)
      if (is.finite(sse_new) && (sse_new <= sse_prev)) {
        th <- th_new
        sse_prev <- sse_new
        improved <- TRUE
        break
      }
      step <- step / 2
    }
    # Hentikan jika tidak ada perbaikan atau langkah sangat kecil
    if (!improved || sqrt(sum((step * Delta)^2)) < tol) break
  }
  th
}

# - Coba beberapa starting values (hindari beta0=0) -
starts <- list(c(0.5,-0.5), c(5, 0.5), c(1, 0.0))
results <- lapply(starts, function(st) gn_fit_safe(x, y, st))
names(results) <- c("start (0.5,-0.5)", "start (5,0.5)", "start (1,0)")
results

## $`start (0.5,-0.5)`
## [1] 1.9753264 0.3148085
## 
## $`start (5,0.5)`
## [1] 1.9753264 0.3148085
## 
## $`start (1,0)`
## [1] 1.9753264 0.3148085

Kesimpulan

Hasil akhirnya relatif konsisten (mendekati \(\beta_0 \approx 2, \beta_1 \approx 0.3\)) meskipun starting value berbeda. Namun, jumlah iterasi dan kecepatan konvergensi bisa berbeda.

2. Bandingkan dengan Transformasi Log

Untuk data dengan error multiplikatif, model bisa dilog-kan agar menjadi linear.

set.seed(456)

y_mult <- beta0_true * exp(beta1_true * x) * exp(rnorm(length(x), sd = 0.2))
df <- data.frame(x, y_mult)

# Estimasi OLS pada log(y)
fit_lm <- lm(log(y_mult) ~ x, data = df)
summary(fit_lm)

## 
## Call:
## lm(formula = log(y_mult) ~ x, data = df)
## 
## Residuals:
##        1        2        3        4        5        6 
## -0.20616  0.19319  0.23530 -0.19636 -0.05516  0.02919 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.63060    0.15252   4.135  0.01444 * 
## x            0.29371    0.05038   5.830  0.00431 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2107 on 4 degrees of freedom
## Multiple R-squared:  0.8947, Adjusted R-squared:  0.8684 
## F-statistic: 33.99 on 1 and 4 DF,  p-value: 0.004312

# konversi kembali
beta0_hat <- exp(coef(fit_lm)[1])
beta1_hat <- coef(fit_lm)[2]
c(beta0_hat, beta1_hat)

## (Intercept)           x 
##   1.8787395   0.2937094

Kesimpulan:

Dengan error multiplikatif, model linear log cocok → estimasi lebih sederhana (OLS), hasilnya dekat dengan nilai parameter sebenarnya.

3. Diagnostik Residu

Setelah fitting, cek residu.

fit_nls <- nls(
  y ~ b0 * exp(b1 * x),
  start = list(b0 = 1, b1 = 0),
  trace = FALSE
)

resid_fit  <- resid(fit_nls)
fitted_fit <- fitted(fit_nls)

par(mfrow = c(1, 2))

plot(
  fitted_fit, resid_fit,
  main = "Residu vs Fitted",
  xlab = "Fitted", ylab = "Residuals"
)
abline(h = 0, lty = 2, col = "red")

qqnorm(resid_fit)
qqline(resid_fit, col = "blue")

par(mfrow=c(1,1))

Interpretasi:

Plot residu vs fitted → apakah ada pola sistematis? Jika tidak, asumsi homoskedastisitas cukup terpenuhi. QQ-plot → memeriksa apakah distribusi residu mendekati normal.

4. Tambahan: Model Michaelis–Menten

Gunakan nls() dengan data simulasi Michaelis–Menten.

set.seed(789)

Vmax_true <- 2
Km_true   <- 1.2

X <- seq(0, 6, length.out = 40)
Y <- (Vmax_true * X) / (Km_true + X) + rnorm(length(X), sd = 0.05)

df_mm <- data.frame(X, Y)

fit_mm <- nls(
  Y ~ (Vmax * X) / (Km + X),
  data  = df_mm,
  start = list(Vmax = 1.5, Km = 1)
)

summary(fit_mm)

## 
## Formula: Y ~ (Vmax * X)/(Km + X)
## 
## Parameters:
##      Estimate Std. Error t value Pr(>|t|)    
## Vmax   1.9981     0.0234   85.40   <2e-16 ***
## Km     1.2415     0.0491   25.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03623 on 38 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 1.317e-06

confint(fit_mm)

##          2.5%   97.5%
## Vmax 1.952336 2.04651
## Km   1.146484 1.34395

Interpretasi

• \(\hat{V}_{\max}\) mendekati 2 = kapasitas maksimum.
• \(\hat{K}_m\) mendekati 1.2 = konsentrasi yang menghasilkan setengah \(V_{\max}\).
• CI (interval kepercayaan) memberi gambaran ketidakpastian estimasi.

Ringkasan Jawaban

• GN & LM sensitif terhadap tebakan awal, tetapi dengan tebakan masuk akal tetap konvergen.
• Untuk error multiplikatif, model bisa dilog-kan menjadi linear → estimasi mudah dengan OLS.
• Diagnostik residu penting untuk validasi asumsi normalitas & homoskedastisitas.
• Pada Michaelis–Menten, interpretasi parameter \(V_{\max}\) dan \(K_m\) sangat substantif (kapasitas & efisiensi sistem).

2.4.2 Inferensi: uji signifikansi parameter, interval kepercayaan,goodness-of-fit (AIC, R²).

Setelah parameter model nonlinear̂ 𝜃 diperoleh (misal dengan nonlinear least squares), kita ingin melakukan inferensi:

1. Apakah parameter signifikan (￿ 0)?
2. Bagaimana ketidakpastian estimasi (interval kepercayaan)?
3. Seberapa baik model cocok dengan data (goodness-of-fit)?

2.4.2.1 Uji Signifikansi Parameter

Hipotesis umum untuk tiap parameter \(\theta_j\):

\(H_0 \colon \theta_j = 0\) vs \(H_1 \colon \theta_j \ne 0\).

Statistik uji Wald: \[t_j = \frac{\hat{\theta}_j}{SE(\hat{\theta}_j)},\] dengan \(SE(\hat{\theta}_j)\) diperoleh dari akar diagonal matriks kovarianŝ \(\operatorname{Var}(\hat{\theta})\).

Jika \(|t_j|\) cukup besar, atau \(p\)-value < 0.05, maka \(\theta_j\) signifikan.

2.4.2.2 Interval Kepercayaan

Wald confidence interval 95% untuk tiap parameter: \[\hat{\theta}_j \pm 1.96 \cdot SE(\hat{\theta}_j).\]

Lebih umum, untuk taraf kepercayaan \((1 - \alpha)\): \[\hat{\theta}_j \pm z_{1-\alpha/2} \cdot SE(\hat{\theta}_j).\]

Dengan \(z_{1-\alpha/2}\) = kuantil distribusi normal standar.

2.4.2.3 Goodness-of-Fit

𝑅² (Pseudo-𝑅² pada regresi nonlinear)
\[R^2 = 1 - \dfrac{\text{SSE}}{\text{SST}},\]

dengan:
- SSE = \(\sum (y_i - \hat{y}_i)^2\) (residual sum of squares),
- SST = \(\sum (y_i - \bar{y})^2\) (total sum of squares).

Interpretasi: proporsi variasi 𝑦 yang dijelaskan model. Nilai mendekati 1 → fit baik.

Akaike Information Criterion (AIC)
\[AIC = n \cdot \ln\!\left(\dfrac{\text{SSE}}{n}\right) + 2p,\]

dengan:
- \(n\) = jumlah observasi,
- \(p\) = jumlah parameter.

Interpretasi: lebih kecil → model lebih baik (dengan penalti kompleksitas).

2.4.2.4 Studi Kasus: Model Eksponensial

Kita gunakan model: \[y = \beta_0 e^{\beta_1 x} + \varepsilon,\ \varepsilon \sim N(0,\ \sigma^2).\]

set.seed(123)
library(broom)
library(ggplot2)

# Simulasi data
beta0_true <- 2; beta1_true <- 0.3
x <- 0:10
y <- beta0_true * exp(beta1_true * x) + rnorm(length(x), sd = 0.5)
dat <- data.frame(x, y)

# Estimasi dengan nls()
fit <- nls(
  y ~ b0 * exp(b1 * x),
  data  = dat,
  start = list(b0 = 1, b1 = 0.1)
)
summary(fit)

## 
## Formula: y ~ b0 * exp(b1 * x)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## b0 2.005728   0.096821   20.72 6.66e-09 ***
## b1 0.300091   0.005399   55.59 9.93e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.515 on 9 degrees of freedom
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 4.468e-06

1. Uji Signifikansi Parameter

# Ringkasan dengan uji t dan p-value
tidy(fit, conf.int = TRUE)

## # A tibble: 2 × 7
##   term  estimate std.error statistic  p.value conf.low conf.high
##   <chr>    <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 b0       2.01    0.0968       20.7 6.66e- 9    1.79      2.23 
## 2 b1       0.300   0.00540      55.6 9.93e-13    0.288     0.313

Interpretasi: jika p-value < 0.05 → parameter signifikan berbeda dari 0.

2. Interval Kepercayaan 95%

confint(fit)

##         2.5%     97.5%
## b0 1.7917267 2.2341574
## b1 0.2879737 0.3126612

Interval kepercayaan memberi rentang nilai parameter yang konsisten dengan data.

3. Goodness-of-Fit

n   <- nrow(dat)
p   <- length(coef(fit))
yhat <- fitted(fit)

SSE <- sum(resid(fit)^2)
SST <- sum((y - mean(y))^2)
R2  <- 1 - SSE / SST
AIC <- n * log(SSE / n) + 2 * p

list(SSE = SSE, R2 = R2, AIC = AIC)

## $SSE
## [1] 2.387234
## 
## $R2
## [1] 0.9984577
## 
## $AIC
## [1] -12.80536

2.5 Tugas

2.5.1 Studi Data Nyata: Ambil dataset pertumbuhan tanaman/hewan (atau data ekonomi) yang menunjukkan pola kurva nonlinear.

• Fit-kan model nonlinear dengan nls().

• Laporkan estimasi parameter, interpretasi biologis/ekonomisnya, serta goodness-of-fit.

Jawaban :

1. Pilih dataset nyata

Pakai dataset bawaan R: Puromycin (kinetika enzim, cocok untuk Michaelis–Menten). Ambil subset state == “treated”.

2. Spesifikasi model

\[Y_i = \dfrac{V_{\max}\, X_i}{K_m + X_i} + \varepsilon_i$, dengan parameter $\theta = (V_{\max}, K_m)^\top.\]

3. Residual & Jacobian

• Residual: \(r(\theta)=y-\mu(\theta)\), \(\ \ \mu(\theta)_i = V_{\max} X_i/(K_m + X_i)\)

• Jacobian \(J(\theta)\) (kolom berurutan untuk \(V_{\max}\), \(K_m\)):

\[\displaystyle \frac{\partial \mu_i}{\partial V_{\max}} = \frac{X_i}{K_m + X_i}, \quad \frac{\partial \mu_i}{\partial K_m} = -\frac{V_{\max} X_i}{(K_m + X_i)^2}\]

4. Gauss–Newton (manual)

Ulangi sampai konvergen:

\((J^\top J)\,\Delta = J^\top r\), lalu \(\theta \leftarrow \theta + \Delta\). Hentikan saat \(\lVert \Delta \rVert_\infty\) kecil.

5. Inferensi (manual)

\(\hat{\sigma}^2 = \text{SSE}/(n-p)\). Kovarian \(\operatorname{Var}(\hat{\theta}) \approx \hat{\sigma}^2 (J^\top J)^{-1}\).
SE = akar diagonal, CI Wald 95%: \(\hat{\theta}_j \pm 1.96\, SE(\hat{\theta}_j)\).

6. Goodness-of-fit (manual)

\(R^2 = 1 - \text{SSE}/\text{SST}\), \(\ \ AIC = n\log(\text{SSE}/n) + 2p\).

7. Diagnostik

Plot residu vs fitted dan QQ-plot.

# 1) Data
data(Puromycin)  # dari paket datasets (bawaan R)
dat <- subset(Puromycin, state == "treated")
x <- dat$conc
y <- dat$rate
n <- length(y)

# 2) Model & Jacobian (sesuai langkah manual)
f_mm <- function(x, theta) {
  Vmax <- theta[1]; Km <- theta[2]
  Vmax * x / (Km + x)
}

J_mm <- function(x, theta) {
  Vmax <- theta[1]; Km <- theta[2]
  cbind(
    d_Vmax = x / (Km + x),
    d_Km   = - Vmax * x / (Km + x)^2
  )
}

# 3) Tebakan awal (rule of thumb yang wajar dari data)
V0 <- max(y)                           # kira-kira Vmax
Km0 <- x[which.min(abs(y - V0/2))]     # x saat respon ~ setengah Vmax
theta0 <- c(Vmax = V0, Km = Km0)

# 4) Gauss–Newton (manual, tanpa nls)
gn_fit <- function(x, y, theta_init, tol = 1e-10, maxit = 200) {
  theta <- as.numeric(theta_init)
  for (k in 1:maxit) {
    mu <- f_mm(x, theta)
    r  <- y - mu
    J  <- J_mm(x, theta)

    # Δ = argmin ||J Δ - r||^2  → gunakan QR (lebih stabil)
    Delta <- qr.coef(qr(J), r)

    # Update
    theta_new <- theta + as.numeric(Delta)

    # Cek konvergensi
    if (max(abs(Delta)) < tol) {
      mu <- f_mm(x, theta_new)
      r  <- y - mu
      return(list(
        par = theta_new,
        mu = mu,
        r = r,
        J = J_mm(x, theta_new),
        iter = k,
        converged = TRUE
      ))
    }
    theta <- theta_new
  }
  mu <- f_mm(x, theta)
  r  <- y - mu
  list(par = theta, mu = mu, r = r, J = J_mm(x, theta), iter = maxit, converged = FALSE)
}

fit <- gn_fit(x, y, theta0)

# 5) Inferensi (SE & CI Wald)
p <- 2
SSE <- sum(fit$r^2)
SST <- sum((y - mean(y))^2)
sigma2_hat <- SSE / (n - p)
# (J'J) di evaluasi di θ_hat (gunakan J terbaru)
JTJ <- crossprod(fit$J)
cov_theta <- as.matrix(sigma2_hat * solve(JTJ))
se <- sqrt(diag(cov_theta))
ci95 <- cbind(
  lower = fit$par - 1.96 * se,
  est   = fit$par,
  upper = fit$par + 1.96 * se
)
rownames(ci95) <- c("Vmax", "Km")

# 6) Goodness-of-fit
R2  <- 1 - SSE / SST
AIC <- n * log(SSE / n) + 2 * p

# 7) Ringkasan hasil
cat("\n=== Gauss–Newton (manual) pada Puromycin: treated ===\n")

## 
## === Gauss–Newton (manual) pada Puromycin: treated ===

cat("Konvergen:", fit$converged, " | Iterasi:", fit$iter, "\n")

## Konvergen: TRUE  | Iterasi: 12

cat(sprintf("Vmax_hat = %.4f (SE = %.4f)\n", fit$par[1], se[1]))

## Vmax_hat = 212.6837 (SE = 6.9472)

cat(sprintf("Km_hat   = %.4f (SE = %.4f)\n", fit$par[2], se[2]))

## Km_hat   = 0.0641 (SE = 0.0083)

cat(sprintf("SSE = %.4f,  R^2 = %.4f,  AIC = %.2f\n", SSE, R2, AIC))

## SSE = 1195.4488,  R^2 = 0.9613,  AIC = 59.22

cat("\nCI 95% (Wald):\n"); print(ci95, digits = 4)

## 
## CI 95% (Wald):

##          lower       est     upper
## Vmax 199.06732 212.68374 226.30017
## Km     0.04789   0.06412   0.08035

# 8) Diagnostik sederhana
par(mfrow = c(1, 2))
plot(fit$mu, fit$r,
     xlab = "Fitted", ylab = "Residuals",
     main = "Residu vs Fitted")
abline(h = 0, lty = 2, col = "gray40")

qqnorm(fit$r, main = "QQ-plot Residu")
qqline(fit$r, col = "gray40")

par(mfrow = c(1, 1))

# 9) (Opsional) Visualisasi kurva fit vs data
plot(x, y, pch = 19, main = "Puromycin (treated): Data vs Fit",
     xlab = "conc", ylab = "rate")
ord <- order(x)
lines(x[ord], fit$mu[ord], lwd = 2)

---

2.5.2 Bandingkan dengan Regresi Linear:

• Fit model linear sederhana pada data yang sama.
• Bandingkan SSE, R², dan visualisasi dengan model nonlinear.
• Diskusikan mengapa model nonlinear lebih sesuai.

Jawaban :

Kita akan membandingkan dua pendekatan untuk data pertumbuhan eksponensial:

1. Nonlinear: \(y(t) = b_0 e^{b_1 t} + \varepsilon\) (diestimasi dengan nls()),

2. Linear sederhana: \(y(t) = a_0 + a_1 t + \varepsilon\) (diestimasi dengan OLS).

Karena pola datanya melengkung (eksponensial), model nonlinear biasanya lebih cocok. Kita bandingkan SSE, \(R^2\), dan AIC, serta lihat perbandingan kurva dan diagnostik residu.

Rumus Model nonlinear (eksponensial):

\[ y(t) \;=\; b_0\,e^{\,b_1 t} + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0,\sigma^2). \]

Model linear sederhana:

\[ y(t) \;=\; a_0 + a_1 t + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0,\sigma^2). \]

Metrik perbandingan:

\[ \text{SSE}=\sum_i (y_i-\hat y_i)^2,\quad R^2 = 1 - \frac{\text{SSE}}{\text{SST}},\ \ \text{SST}=\sum_i (y_i-\bar y)^2, \]

\[ \text{AIC} = n \cdot \ln\!\Big(\frac{\text{SSE}}{n}\Big) + 2p, \]

dengan \(n\) jumlah observasi dan \(p\) jumlah parameter dalam model.

set.seed(202)

# --- Simulasi data eksponensial: y = b0 * exp(b1 * t) + e ---
b0_true <- 1.5
b1_true <- 0.4
sigma    <- 0.2

t      <- seq(0, 6, by = 0.2)
y_true <- b0_true * exp(b1_true * t)
y      <- y_true + rnorm(length(t), sd = sigma)
dat    <- data.frame(t = t, y = y)

# --- Fit Nonlinear (nls): y = b0 * exp(b1 * t) ---
fit_nl <- nls(
  y ~ b0 * exp(b1 * t),
  data  = dat,
  start = list(b0 = 1, b1 = 0.1),
  trace = FALSE
)

# --- Fit Linear sederhana: y = a0 + a1 * t ---
fit_lin <- lm(y ~ t, data = dat)

# --- Metrik: SSE, R^2, AIC (keduanya) ---
SST <- sum( (y - mean(y))^2 )

# Nonlinear
yhat_nl <- fitted(fit_nl)
res_nl  <- resid(fit_nl)
SSE_nl  <- sum(res_nl^2)
R2_nl   <- 1 - SSE_nl / SST
AIC_nl  <- length(y) * log(SSE_nl / length(y)) + 2 * length(coef(fit_nl))

# Linear
yhat_lin <- fitted(fit_lin)
res_lin  <- resid(fit_lin)
SSE_lin  <- sum(res_lin^2)
R2_lin   <- 1 - SSE_lin / SST
AIC_lin  <- length(y) * log(SSE_lin / length(y)) + 2 * length(coef(fit_lin))

# --- Tampilkan output ringkas ---
cat("== Ringkasan Koefisien ==\n\n")

## == Ringkasan Koefisien ==

cat("Nonlinear (nls):\n"); print(summary(fit_nl)$coefficients)

## Nonlinear (nls):

##     Estimate  Std. Error   t value     Pr(>|t|)
## b0 1.4827376 0.028640593  51.77049 3.960921e-30
## b1 0.4050742 0.003839203 105.50997 4.790746e-39

cat("\nLinear (lm):\n");    print(summary(fit_lin)$coefficients)

## 
## Linear (lm):

##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) -0.7457131  0.4668751 -1.597243 1.210523e-01
## t            2.3834525  0.1336660 17.831410 3.592662e-17

cat("\n== Perbandingan Metrik ==\n")

## 
## == Perbandingan Metrik ==

cmp <- data.frame(
  Model = c("Nonlinear (nls)", "Linear (OLS)"),
  SSE   = c(SSE_nl, SSE_lin),
  R2    = c(R2_nl,  R2_lin),
  AIC   = c(AIC_nl, AIC_lin)
)
print(cmp, row.names = FALSE)

##            Model       SSE        R2        AIC
##  Nonlinear (nls)  1.018542 0.9983437 -101.88408
##     Linear (OLS) 51.398598 0.9164167   19.67433

# Kurva 'true' (ground truth) dan kurva hasil fit
tt <- seq(min(dat$t), max(dat$t), length.out = 300)

# Kurva true
yy_true <- b0_true * exp(b1_true * tt)

# Kurva fitted Nonlinear
coef_nl <- coef(fit_nl)
yy_nl   <- coef_nl[["b0"]] * exp(coef_nl[["b1"]] * tt)

# Kurva fitted Linear
yy_lin <- predict(fit_lin, newdata = data.frame(t = tt))

# Plot
plot(dat$t, dat$y, pch = 19, xlab = "t", ylab = "y",
     main = "Perbandingan: Nonlinear (nls) vs Linear (OLS)")
lines(tt, yy_true, lwd = 2)             # true
lines(tt, yy_nl,  lwd = 2, lty = 2)     # nonlinear fit
lines(tt, yy_lin, lwd = 2, lty = 3)     # linear fit
legend("topleft",
       legend = c("Data", "True", "Fitted (Nonlinear)", "Fitted (Linear)"),
       pch = c(19, NA, NA, NA), lty = c(NA, 1, 2, 3), lwd = c(NA, 2, 2, 2),
       bty = "n")

par(mfrow = c(2, 2))

# Nonlinear
plot(yhat_nl, res_nl,
     xlab = "Fitted (Nonlinear)", ylab = "Residuals",
     main = "Residu vs Fitted — Nonlinear")
abline(h = 0, lty = 2, col = "gray40")
qqnorm(res_nl, main = "QQ-plot — Nonlinear")
qqline(res_nl, col = "gray40")

# Linear
plot(yhat_lin, res_lin,
     xlab = "Fitted (Linear)", ylab = "Residuals",
     main = "Residu vs Fitted — Linear")
abline(h = 0, lty = 2, col = "gray40")
qqnorm(res_lin, main = "QQ-plot — Linear")
qqline(res_lin, col = "gray40")

par(mfrow = c(1, 1))

• SSE dan AIC biasanya lebih baik pada model nonlinear (nls) karena bentuk eksponensial menangkap kelengkungan data.

• \(R^2\) cenderung lebih tinggi pada model nonlinear untuk pola pertumbuhan seperti ini.

• Visualisasi menunjukkan garis linear sulit mengikuti kenaikan yang semakin cepat, sementara kurva eksponensial dapat menyesuaikan laju pertumbuhan.

• Diagnostik residu: untuk model nonlinear, residu biasanya lebih “acak” di sekitar nol; pada model linear sering terlihat pola (misalnya berbentuk lengkung), tanda misspecification.

Kesimpulan: Untuk data yang mengikuti pola eksponensial (atau umumnya nonlinear), model nonlinear lebih sesuai dibandingkan regresi linear sederhana.

---

2.5.3 Simulasi: Buat simulasi data eksponensial (misal: pertumbuhan bakteri).

• Gunakan nls() untuk estimasi parameter.
• Laporkan hasil dan interpretasi.

Jawaban :

Kita mensimulasikan data dari proses pertumbuhan eksponensial dengan galat aditif, lalu menaksir parameternya menggunakan Nonlinear Least Squares (nls). Setelah itu kita hitung metrik kecocokan (SSE, \(R^2\), AIC), buat interval kepercayaan Wald 95%, dan tampilkan diagnostik residu. Tujuannya: memastikan model menangkap pola kenaikan eksponensial dengan baik.

Rumus model Model yang digunakan:

\[ y(t) \;=\; b_0\,e^{\,b_1 t} \;+\; \varepsilon,\qquad \varepsilon \sim \mathcal{N}(0,\sigma^2). \]

• \(b_0\): level dasar saat \(t=0\) \(\big(y(0)=b_0+\varepsilon\big)\)
  
• \(b_1\): laju pertumbuhan (semakin besar nilainya, pertumbuhan makin cepat)
  
• \(\sigma\): simpangan baku galat

set.seed(202)

# 1) Simulasi data: y = b0 * exp(b1 * t) + e
b0_true <- 1.5
b1_true <- 0.4
sigma    <- 0.2

t      <- seq(0, 6, by = 0.2)
y_true <- b0_true * exp(b1_true * t)
y      <- y_true + rnorm(length(t), sd = sigma)
dat    <- data.frame(t = t, y = y)

# 2) Estimasi NLS
fit <- nls(
  y ~ b0 * exp(b1 * t),
  data  = dat,
  start = list(b0 = 1, b1 = 0.1),
  trace = FALSE
)

# 3) Ringkasan koefisien, SE, dan CI Wald 95%
sm       <- summary(fit)
coef_hat <- coef(fit)
se_hat   <- sm$parameters[, "Std. Error"]
ci95     <- cbind(
  lower = coef_hat - 1.96 * se_hat,
  est   = coef_hat,
  upper = coef_hat + 1.96 * se_hat
)

# 4) Goodness-of-fit: SSE, R^2 (pseudo), AIC
yhat <- fitted(fit)
res  <- resid(fit)
SSE  <- sum(res^2)
SST  <- sum((y - mean(y))^2)
R2   <- 1 - SSE / SST
AIC  <- length(y) * log(SSE / length(y)) + 2 * length(coef_hat)

# 5) Tampilkan output ringkas
cat("== Koefisien (taksiran) dan SE ==\n")

## == Koefisien (taksiran) dan SE ==

print(sm$coefficients)

##     Estimate  Std. Error   t value     Pr(>|t|)
## b0 1.4827376 0.028640593  51.77049 3.960921e-30
## b1 0.4050742 0.003839203 105.50997 4.790746e-39

cat("\n== CI 95% (Wald) ==\n")

## 
## == CI 95% (Wald) ==

print(ci95, digits = 4)

##     lower    est  upper
## b0 1.4266 1.4827 1.5389
## b1 0.3975 0.4051 0.4126

cat(sprintf("\n== Goodness-of-fit ==\nSSE = %.4f  |  R^2 = %.4f  |  AIC = %.2f\n", SSE, R2, AIC))

## 
## == Goodness-of-fit ==
## SSE = 1.0185  |  R^2 = 0.9983  |  AIC = -101.88

plot(dat$t, dat$y, pch = 19, xlab = "t", ylab = "y",
     main = "Data simulasi vs kurva eksponensial")

tt <- seq(min(dat$t), max(dat$t), length.out = 200)

# Kurva 'true' (parameter ground truth)
lines(tt, b0_true * exp(b1_true * tt), lwd = 2)

# Kurva hasil estimasi nls
lines(tt, coef_hat[["b0"]] * exp(coef_hat[["b1"]] * tt), lwd = 2, lty = 2)

legend("topleft",
       legend = c("Data", "True", "Fitted (nls)"),
       pch = c(19, NA, NA), lty = c(NA, 1, 2), lwd = c(NA, 2, 2),
       bty = "n")

par(mfrow = c(1, 2))

plot(yhat, res,
     xlab = "Fitted", ylab = "Residuals",
     main = "Residu vs Fitted")
abline(h = 0, lty = 2, col = "gray40")

qqnorm(res, main = "QQ-plot Residu")
qqline(res, col = "gray40")

par(mfrow = c(1, 1))

• Koefisien \(\hat{b}_0\) dan \(\hat{b}_1\) diharapkan mendekati nilai sebenarnya \((1.5,\, 0.4)\) karena data memang disimulasikan dari model tersebut.

• SSE kecil dan \(R^2\) tinggi menandakan model menangkap pola eksponensial dengan baik. AIC jadi pembanding saat menilai alternatif model.

• CI 95% memberi ketidakpastian taksiran. Jika CI untuk \(\hat{b}_1\) tidak memuat 0, ada bukti kuat bahwa laju pertumbuhan eksponensial signifikan.

• Pada Residu vs Fitted, sebaran acak di sekitar nol memberi indikasi homoskedastisitas yang wajar. Pada QQ-plot, titik yang mengikuti garis mengindikasikan normalitas residu yang cukup baik.

---

2.6 Kesimpulan

• Regresi nonlinear diperlukan ketika hubungan X–Y tidak bisa dijelaskan dengan garis lurus.
• Parameter nonlinear sering punya interpretasi substantif (misalnya: laju pertumbuhan, kapasitas maksimum).
• Estimasi dilakukan dengan algoritma iteratif (nls() di R).
• Visualisasi sangat penting untuk membandingkan data vs kurva model.

2.7 Appendix: A. Gauss–Newton (GN) untuk Model Eksponensial

2.7.1 Model dan Asumsi

Model regresi eksponensial dengan error aditif: \[Y_i = f_i(\theta) + \varepsilon_i = \beta_0 e^{\beta_1 x_i} + \varepsilon_i,\ \varepsilon_i\ \text{iid} \sim \mathcal{N}(0,\ \sigma^2),\ i = 1, \ldots, n.\]

Vektor parameter:

\[\theta = (\beta_0, \beta_1)^\top,\ p = 2.\]

Catatan: Di bawah asumsi normal, NLS = MLE untuk \(\theta\) dan \(\sigma^2\).

2.7.2 Residual, Fungsi Tujuan, dan Jacobian

Residual dan fungsi objektif (jumlah kuadrat residual): \[r_i(\theta) = y_i - \beta_0 e^{\beta_1 x_i}, \quad S(\theta) = \sum_{i=1}^{n} [r_i(\theta)]^2 = \lVert r(\theta) \rVert_2^2.\]

Jacobian \(J(\theta) \in \mathbb{R}^{n \times p}\) dengan baris ke-\(i\):

\[J_{i1} = \dfrac{\partial f_i}{\partial \beta_0} = e^{\beta_1 x_i}, \quad J_{i2} = \dfrac{\partial f_i}{\partial \beta_1} = \beta_0 x_i e^{\beta_1 x_i}.\]

2.7.3 Gradien dan (Aproksimasi) Hessian dari 𝑆(𝜃)

Gradien di \(\theta\): \[\nabla S(\theta) = g(\theta) = -2\, J(\theta)^\top r(\theta).\]

Aproksimasi Hessian Gauss–Newton: \[H_{\text{GN}}(\theta) \approx 2\, J(\theta)^\top J(\theta).\]

(Hessian eksak: \(H = 2 J^\top J - 2 \sum_{i=1}^{n} r_i H_i\), dengan \(H_i\) Hessian \(f_i\).)

2.7.4 Satu Langkah Iterasi GN (Linearized Least Squares)

Pada iterasi ke-\(k\), langkah perbaikan \(\Delta^{(k)}\) diperoleh dari: \[(J^{(k)\top} J^{(k)})\,\Delta^{(k)} = J^{(k)\top} r^{(k)},\] dengan \[r^{(k)} = \bigl(y_i - \beta_0^{(k)} e^{\beta_1^{(k)} x_i}\bigr)_{i=1}^n,\quad J^{(k)} = J(\theta^{(k)}).\]

Pembaruan parameter: \[\theta^{(k+1)} = \theta^{(k)} + \Delta^{(k)}.\]

Kriteria henti (contoh): \[\lVert \Delta^{(k)} \rVert_\infty < \varepsilon_\theta \quad \text{atau} \quad \bigl|S^{(k+1)} - S^{(k)}\bigr| < \varepsilon_S \quad \text{atau} \quad \lVert g^{(k)} \rVert_\infty < \varepsilon_g.\]

Alternatif redaman (Levenberg–Marquardt): \[(J^{(k)\top} J^{(k)} + \lambda^{(k)} D^{(k)})\,\Delta^{(k)} = J^{(k)\top} r^{(k)},\quad D^{(k)} = I \ \text{atau}\ \operatorname{diag}(J^{(k)\top} J^{(k)}).\]

2.7.5 Estimasi Varians Error dan Kovarians Parameter

Setelah konvergen dî \(\theta\), taksiran varians error: \[\hat{\sigma}^2 = \dfrac{S(\hat{\theta})}{n - p}.\]

Kovarians asimtotik parameter (Wald): \[\widehat{\mathrm{Cov}}(\hat{\theta}) \approx \hat{\sigma}^2 \left(J(\hat{\theta})^\top J(\hat{\theta})\right)^{-1}.\]

Galat baku: \[SE(\hat{\theta}_j) = \sqrt{\hat{\sigma}^2 \left[\left(J^\top J\right)^{-1}\right]_{jj}}\ \big|_{\theta=\hat{\theta}}.\]

(Opsional, robust sandwich—untuk heteroskedastisitas): \[\widehat{\text{Covrobust}}(\hat{\theta}) = \left(J^\top J\right)^{-1} J^\top \operatorname{diag}\!\big(r_i(\hat{\theta})^2\big)\, J \left(J^\top J\right)^{-1}.\]

2.7.6 Inferensi: Statistik-𝑧, CI Wald, dan Uji Wald Bersama

Asimtotik normalitas (Wald): \[\hat{\theta} \ \dot{\sim}\ \mathcal{N}\!\big(\theta,\ \sigma^2 (J^\top J)^{-1}\big).\]

Statistik-\(z\) satu parameter untuk hipotesis \(H_0 \colon \theta_j = \theta_{j,0}\): \[z_j = \frac{\hat{\theta}_j - \theta_{j,0}}{SE(\hat{\theta}_j)} \ \dot{\sim}\ \mathcal{N}(0,1).\]

CI Wald \(100(1-\alpha)\%\) untuk \(\theta_j\): \[\hat{\theta}_j \pm z_{1-\alpha/2}\, SE(\hat{\theta}_j),\] dengan \(z_{1-\alpha/2}\) kuantil standar normal (mis. \(z_{0.975} \approx 1.96\)).

(Pada sampel kecil, dapat memakai \(t_{n-p,\ 1-\alpha/2}\).)

Uji Wald bersama (\(q\) kendala linear \(R\theta = r\)):

\[W = (R\hat{\theta} - r)^\top \big[\, R\, \widehat{\mathrm{Cov}}(\hat{\theta})\, R^\top \big]^{-1} (R\hat{\theta} - r)\ \dot{\sim}\ \chi^2_q.\]

2.7.7 Interval untuk Rata-rata Kondisional dan Prediksi

Untuk titik baru \(x^\star\), prediksi mean: \[\hat{y}^\star = f(x^\star, \hat{\theta}) = \hat{\beta}_0 \, e^{\hat{\beta}_1 x^\star}.\]

Linearization delta method: turunan terhadap parameter dî \(\theta\), \[g^\star = \left[\frac{\partial f(x^\star,\theta)}{\partial \beta_0}\ \ \frac{\partial f(x^\star,\theta)}{\partial \beta_1}\right]_{\theta=\hat{\theta}} = \left[ e^{\hat{\beta}_1 x^\star}\ \ \hat{\beta}_0 x^\star e^{\hat{\beta}_1 x^\star} \right].\]

Varian prediksi mean: \[\operatorname{Var}(\hat{y}^\star) \approx {g^\star}^\top\, \widehat{\operatorname{Cov}}(\hat{\theta})\, g^\star.\]

CI untuk mean \(f(x^\star,\theta)\): \[\hat{y}^\star \pm z_{1-\alpha/2}\, \sqrt{\operatorname{Var}(\hat{y}^\star)}.\]

PI (prediction interval) untuk \(Y^\star\): \[\hat{y}^\star \pm z_{1-\alpha/2}\, \sqrt{\operatorname{Var}(\hat{y}^\star) + \hat{\sigma}^2}.\]

2.7.8 Ringkas Algoritma GN (Praktis)

1. Pilih tebakan awal \(\theta^{(0)}\) (mis. via regresi \(\ln y\) pada \(x\) untuk inisialisasi).

2. Ulangi hingga konvergen:

• Hitung \(r^{(k)}\) dan \(J^{(k)}\).
• Pecahkan \((J^{(k)\top} J^{(k)})\,\Delta^{(k)} = J^{(k)\top} r^{(k)}\).
• Set \(\theta^{(k+1)} = \theta^{(k)} + \Delta^{(k)}\).
• Hentikan jika kriteria henti terpenuhi; jika \(S\) tidak turun, gunakan line-search/LM.

3. Padâ \(\theta\), hitung \(\hat{\sigma}^2\), \(\widehat{\text{Cov}}(\hat{\theta})\), SE, \(z\)-stat, dan CI Wald.

Catatan

• Sensitivitas starting value: GN/LM dapat sensitif; gunakan beberapa tebakan dan bandingkan SSE akhir.
• Nonlinearitas kuat: CI Wald bisa kurang akurat; pertimbangkan profil likelihood.
• Skala: Standarisasi \(x\) dapat memperbaiki kondisioning \(J^\top J\).

2.7.9 Perhitungan Manual

Contoh GN Manual (n = 5) untuk Model Eksponensial

Data
Kita gunakan 5 observasi:

\[x = (0, 1, 2, 3, 4)^\top,\]
\[y = (2.100000, 2.649718, 3.724238, 4.819206, 6.790234)^\top.\]

Model
\[y_i = \beta_0 e^{\beta_1 x_i} + \varepsilon_i,\ \ \varepsilon_i \sim \mathcal N(0,\ \sigma^2).\]

Iterasi 0 (tebakan awal dari regresi \(\ln y\) pada \(x\))
\[ \theta^{(0)} = \begin{bmatrix} \beta^{(0)}_0\\[2pt] \beta^{(0)}_1 \end{bmatrix} = \begin{bmatrix} 2.043817\\[2pt] 0.294525 \end{bmatrix}. \]

Nilai fungsi, residual, dan Jacobian di \(\theta^{(0)}\):

\[ f^{(0)} = \beta^{(0)}_0 e^{\beta^{(0)}_1 x} = \begin{bmatrix} 2.043817\\ 2.743801\\ 3.683523\\ 4.945088\\ 6.638725 \end{bmatrix},\qquad r^{(0)} = y - f^{(0)} = \begin{bmatrix} \ \ 0.056183\\ -0.094083\\ \ \ 0.040715\\ -0.125881\\ \ \ 0.151509 \end{bmatrix}. \]

Jacobian (baris ke-\(i\): \([\partial f_i/\partial \beta_0,\ \partial f_i/\partial \beta_1]\)): \[ J^{(0)} = \begin{bmatrix} 1.000000 & \ \ 0.000000\\ 1.342489 & \ \ 2.743801\\ 1.802276 & \ \ 7.367045\\ 2.419536 & 14.835263\\ 3.248200 & 26.554899 \end{bmatrix}. \]

Matriks normal dan vektor kanan: \[ J^{(0)\top}J^{(0)} = \begin{bmatrix} \ 22.455430 & 139.111033\\ 139.111033 & 987.049457 \end{bmatrix},\qquad J^{(0)\top} r^{(0)} = \begin{bmatrix} 0.190815\\ 2.197633 \end{bmatrix}. \]

Langkah GN (persamaan normal): \[ \big(J^{(0)\top}J^{(0)}\big)\,\Delta^{(0)} = J^{(0)\top} r^{(0)} \ \ \Rightarrow\ \ \Delta^{(0)} = \begin{bmatrix} -0.041728\\ \ \ 0.008108 \end{bmatrix}. \]

Pembaruan parameter dan SSE: \[ \theta^{(1)} = \theta^{(0)} + \Delta^{(0)} = \begin{bmatrix} 2.002088\\ 0.302633 \end{bmatrix},\qquad S\big(\theta^{(0)}\big)=0.052467,\ \ S\big(\theta^{(1)}\big)=0.042495. \]

Iterasi 1
Hitung ulang \(J^{(1)},\ r^{(1)}\) pada \(\theta^{(1)}\): \[ J^{(1)} = \begin{bmatrix} 1.000000 & \ \ 0.000000\\ 1.353417 & \ \ 2.709661\\ 1.831738 & \ \ 7.334604\\ 2.479106 & 14.890170\\ 3.355265 & 26.870152 \end{bmatrix},\qquad r^{(1)} = \begin{bmatrix} \ \ 0.097912\\ -0.059943\\ \ \ 0.056936\\ -0.144184\\ \ \ 0.072696 \end{bmatrix}. \]

\[ J^{(1)\top}J^{(1)} = \begin{bmatrix} \ 23.590779 & 144.173187\\ 144.173187 & 1004.860917 \end{bmatrix},\qquad J^{(1)\top} r^{(1)} = \begin{bmatrix} 0.007541\\ 0.061600 \end{bmatrix}. \]

Langkah GN berikutnya: \[ \Delta^{(1)} = \begin{bmatrix} -0.000446\\ \ \ 0.000125 \end{bmatrix},\qquad \theta^{(2)} = \begin{bmatrix} 2.001642\\ 0.302758 \end{bmatrix},\qquad S\big(\theta^{(2)}\big)=0.042491. \]

Terlihat 𝑆(𝜃) menurun dan perubahan parameter makin kecil → konvergen.

2.8 Appendix: Levenberg–Marquardt (LM) untuk Model Eksponensial

2.8.1 Model dan Asumsi

Kita gunakan model regresi eksponensial dengan error aditif: \[Y_i = f_i(\theta) + \varepsilon_i = \beta_0 e^{\beta_1 x_i} + \varepsilon_i,\ \varepsilon_i\ \text{iid} \sim \mathcal{N}(0,\ \sigma^2),\ i = 1, \ldots, n.\]

Vektor parameter: \[\theta = (\beta_0, \beta_1)^\top,\ p = 2.\]

2.8.2 Ide Inti LM (Damped Gauss–Newton)

LM menggabungkan Gauss–Newton dan Gradient Descent.

Iterasi ke-\(k\) menyelesaikan sistem modifikasi:

\[(J^{(k)\top} J^{(k)} + \lambda^{(k)} D^{(k)})\, \Delta^{(k)} = J^{(k)\top} r^{(k)},\]

dengan:

• \(J^{(k)}\) : Jacobian di \(\theta^{(k)}\)
• \(r^{(k)} = (y_i - \beta^{(k)}_0 e^{\beta^{(k)}_1 x_i})_{i=1}^n\)
• \(D^{(k)} = I\) atau \(\operatorname{diag}(J^{(k)\top} J^{(k)})\)
• \(\lambda^{(k)} > 0\) : faktor redaman

Jika \(\lambda^{(k)} \to 0\), metode mendekati GN; jika \(\lambda^{(k)}\) besar, mendekati Gradient Descent.

2.8.3 Residual, Fungsi Tujuan, dan Jacobian

Residual dan fungsi SSE sama dengan GN: \[r_i(\theta) = y_i - \beta_0 e^{\beta_1 x_i}, \quad S(\theta) = \sum_{i=1}^{n} r_i(\theta)^2.\]

Jacobian (baris ke-\(i\)): \[J_{i1} = e^{\beta_1 x_i}, \quad J_{i2} = \beta_0 x_i e^{\beta_1 x_i}.\]

2.8.4 Update Parameter dengan Redaman

Pembaruan parameter: \[\theta^{(k+1)} = \theta^{(k)} + \Delta^{(k)}.\]

Aturan adaptasi \(\lambda\):

• Jika \(S(\theta^{(k+1})) < S(\theta^{(k)})\), kecilkan \(\lambda\) (lebih agresif).
• Jika tidak, besarkan \(\lambda\) (lebih konservatif) dan coba ulangi.

Sehingga: \[ \lambda^{(k+1)} = \begin{cases} \lambda^{(k)}/\nu, & \text{jika } S(\theta^{(k+1})) < S(\theta^{(k)}),\\ \nu\,\lambda^{(k)}, & \text{jika tidak.} \end{cases} \]

(dengan \(\nu > 1\), misalnya \(\nu = 10\)).

2.8.5 Gradien dan (Aproksimasi) Hessian dalam LM

LM memodifikasi matriks normal: \[H_{\text{LM}}^{(k)} = J^{(k)\top} J^{(k)} + \lambda^{(k)} D^{(k)}.\]

Gradien tetap sama: \[g^{(k)} = -2\, J^{(k)\top} r^{(k)}.\]

2.8.6 Estimasi Varians Error dan Kovarians Parameter

Setelah konvergen dî \(\theta\): \[\hat{\sigma}^2 = \dfrac{S(\hat{\theta})}{n - p}.\]

Kovarians asimtotik (Wald): \[\widehat{\mathrm{Cov}}(\hat{\theta}) \approx \hat{\sigma}^2 (J^\top J)^{-1}\big|_{\theta=\hat{\theta}}.\]

SE parameter: \[SE(\hat{\theta}_j) = \sqrt{\hat{\sigma}^2 \,\big[(J^\top J)^{-1}\big]_{jj}}.\]

2.8.7 Inferensi: Statistik-𝑧 dan CI

Distribusi asimtotik: \[\hat{\theta} \ \dot{\sim}\ \mathcal{N}\!\big(\theta,\ \sigma^2 (J^\top J)^{-1}\big).\]

Statistik-\(z\): \[z_j = \frac{\hat{\theta}_j - \theta_{j,0}}{SE(\hat{\theta}_j)} \ \dot{\sim}\ \mathcal{N}(0,1).\]

CI Wald \(100(1 - \alpha)\%\) untuk \(\theta_j\): \[\hat{\theta}_j \pm z_{1-\alpha/2}\, SE(\hat{\theta}_j).\]

2.8.8 Interval Prediksi dan Mean

Prediksi mean di \(x^\star\): \[\hat{y}^\star = f(x^\star, \hat{\theta}) = \hat{\beta}_0 e^{\hat{\beta}_1 x^\star}.\]

Gradien terhadap parameter: \[ g^\star = \begin{bmatrix} e^{\hat{\beta}_1 x^\star} & \hat{\beta}_0 x^\star e^{\hat{\beta}_1 x^\star} \end{bmatrix}. \]

Varian prediksi mean: \[\operatorname{Var}(\hat{y}^\star) \approx {g^\star}^\top\, \widehat{\operatorname{Cov}}(\hat{\theta})\, g^\star.\]

CI mean: \[\hat{y}^\star \pm z_{1-\alpha/2}\, \sqrt{\operatorname{Var}(\hat{y}^\star)}.\]

PI (prediction interval): \[\hat{y}^\star \pm z_{1-\alpha/2}\, \sqrt{\operatorname{Var}(\hat{y}^\star) + \hat{\sigma}^2}.\]

2.8.9 Ringkas Algoritma LM

1. Pilih tebakan awal \(\theta^{(0)}\), pilih \(\lambda^{(0)} > 0\).
2. Iterasi:

• Hitung \(r^{(k)}\) dan \(J^{(k)}\).
• Selesaikan \((J^{(k)\top}J^{(k)} + \lambda^{(k)}D^{(k)})\Delta^{(k)} = J^{(k)\top}r^{(k)}\).
• Update \(\theta^{(k+1)} = \theta^{(k)} + \Delta^{(k)}\).
• Adaptasi \(\lambda\) sesuai aturan sukses/gagal.

3. Hentikan bila kriteria konvergensi tercapai.
4. Setelah konvergen: hitung \(\hat{\sigma}^2\), SE, \(z\)-statistik, dan CI.

Catatan - LM lebih stabil daripada GN ketika \(J^\top J\) hampir singular atau starting value jauh. - \(\lambda\) adaptif menjaga keseimbangan eksplorasi (besar \(\lambda\)) dan efisiensi (kecil \(\lambda\)). - Interval Wald tetap bisa kurang akurat → gunakan profil likelihood bila perlu.