Latihan dan Tugas
Angga Pratama
2025-09-04
Analisis Regresi Tingkat Lanjut
Praktikum Pertama
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:
- Menjelaskan hubungan antara variabel.
- Memprediksi nilai variabel dependen berdasarkan variabel independen.
- 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
• 𝑌𝑖 : variabel dependen (respon)
• 𝑋𝑖 : variabel independen (prediktor)
• 𝛽0, 𝛽1 : parameter regresi
• 𝜀𝑖 : error (residu)
1.1 Konsep Korelasi
Korelasi mengukur derajat keeratan hubungan linier antara dua variabel, tanpa membedakan mana yang dependen dan independen.
Koefisien korelasi Pearson:
\[ 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 𝑟 ∈ [−1, 1].
• 𝑟 > 0: hubungan positif, 𝑟 < 0: hubungan negatif.
• |𝑟| makin mendekati 1 artinya hubungan makin kuat
Berikut tabel perbedaan utama antara korelasi dan regresi:
Tabel Perbandingan
| Aspek | Korelasi | Regresi |
|---|---|---|
| Tujuan | Mengukur keeratan hubungan | Membangun model prediksi/penjelasan |
| Arah hubungan | Simetris (X dan Y setara) | Asimetris (Y sebagai respon, X sebagai prediktor) |
| Output utama | Koefisien korelasi \(r\) | Persamaan regresi \(\hat{Y} = \beta_0 + \beta_1 X\) |
Contoh Aplkasi 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.460Korelasi
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.8676594Regresi
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-10Visualisasi
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 Regesi Linear
Kita memodelkan hubungan antara respons \(y \in \mathbb{R}^n\) dan kovariat \(X \in \mathbb{R}^{n \times p}\):
\[ y = X\beta + \varepsilon,\ \mathbb{E}(\varepsilon)=0,\ \operatorname{Var}(\varepsilon)=\sigma^2 I_n. \]
Dengan \(X = [\mathbf{1},\, x_1,\, \ldots,\, 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:
\[ \widehat{\beta}_{\text{OLS}} = \operatorname*{arg\,min}_{\beta}\,\lVert y - X\beta\rVert^{2} = (X^{\top}X)^{-1} X^{\top} y . \]
Prediksi \(\hat y = X\widehat{\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 = H y, \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), \(\widehat{\beta}_{\text{OLS}}\) adalah BLUE.
Varians koefisien dan estimator ragam:
\[ \operatorname{Var}(\widehat{\beta}_{\text{OLS}})=\sigma^2 (X^{\top}X)^{-1}, \qquad \widehat{\sigma}^2=\frac{\lVert e\rVert^2}{n-p}. \]
Asimtotik: \[ \sqrt{n}\,\big(\widehat{\beta}-\beta\big)\xrightarrow{d}\mathcal{N}(0,\,V), \] sehingga inferensi berbasis normal/\(t\) berlaku untuk \(n\) besar.
1.3 GENERALISASI : WLS & GLS
Jika galat tidak homoskedastik/berkorelasi, maka \(\operatorname{Var}(\varepsilon) = \sigma^2 \Omega \neq \sigma^2 I_n.\)
Weighted least squares (WLS) untuk varians berbeda-beda (\(\Omega\) diagonal): \[ \widehat{\beta}_{\mathrm{WLS}} = (X^{\top} W X)^{-1}\, X^{\top} W y, \qquad W = \Omega^{-1}. \]
Generalized least squares (GLS) untuk \(\Omega\) umum: \[ \widehat{\beta}_{\mathrm{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 (HC0) adalah \[ \widehat{\operatorname{Var}}(\widehat{\beta}) = (X^{\top}X)^{-1} \Bigg(\sum_{i=1}^n x_i x_i^{\top}\,\hat e_i^2\Bigg) (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\): \[ \mathrm{VIF}_j \;=\; \frac{1}{1 - R_j^{2}}, \] dengan \(R_j^{2}\) adalah koefisien determinasi dari regresi \(x_j\) pada semua 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\,\widehat{\sigma}^2}\;\frac{h_{ii}}{(1-h_{ii})^2}, \]
Gunakan plot residu, plot leverage vs residu, QQ-plot, serta inspeksi \(𝐷𝑖\) 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:}\quad \widehat{\beta} = \operatorname*{arg\,min}_{\beta}\;\Big\{\;\|y - X\beta\|_2^2 + \lambda\,\|\beta\|_2^2\;\Big\}. \]
\[ \text{Lasso:}\quad \widehat{\beta} = \operatorname*{arg\,min}_{\beta}\;\Big\{\;\|y - X\beta\|_2^2 + \lambda\,\|\beta\|_1\;\Big\}. \]
1.8 Nonlinearitas: Ekspansi Basis & Spline
Kita dapat menambah term polinomial atau basis spline sehingga tetap linear dalam parameter:
\[ y = \beta_0 + \sum_{k=1}^{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.02031.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")
p2 <- ggplot(diag_df, aes(sample = stdres)) +
stat_qq() + stat_qq_line() + labs(x = "Theoretical", y = "Standardized Residuals")
p1; p2
Figure 1.2: Diagnostik OLS untuk mtcars: residu vs fitted (kiri) dan QQ-plot (kanan).
1.9.4 Multikolinearitas (VIF)
car::vif(fit1)## wt hp disp
## 4.844618 2.736633 7.3245171.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 ' ' 1qplot(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.
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 ' ' 11.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.6430BIC(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
Tanyakan bagian ini belum lengkap
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.281 0.267 0.191
## 2 x1 0.0300 0.0313 0.0267
## 3 x2 0.0435 0.0466 0.0333# Uji Breusch-Pagan untuk heteroskedastisitas
lmtest::bptest(fit_ols)##
## studentized Breusch-Pagan test
##
## data: fit_ols
## BP = 48.294, df = 2, p-value = 3.259e-11ggplot(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")
1.11 Catatan
- Mulai dari OLS + diagnostik; bila heteroskedastik, pilih robust SE atau WLS/GLS sesuai struk- tur ragam.
- Periksa multikolinearitas (VIF), leverage, Cook’s 𝐷, dan normalitas residu.
- Gunakan ekspansi basis (polinomial/spline) saat ada indikasi nonlinearitas.
- Validasi model (CV) dan bandingkan AIC/BIC untuk kompleksitas vs akurasi.
- 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\) dengan \(u_t \sim \mathcal{N}(0,\ \sigma_u^2)\). Kovarians residual memiliki bentuk Toeplitz:\(\Omega_{ts} = \frac{\sigma_u^2}{1-\rho^2}\,\rho^{|t-s|}\) 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.2131.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-16summary(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 residual1.12.3 Diagnostik serial correlation
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 0par(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
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-16Catatan: 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.430232. 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-16broom::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 ' ' 13.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 ' ' 14. Uji Asumsi
4.1 Multikolinearitas
# Variance Inflation Factor (VIF)
car::vif(mod0)## X1 X2 X3
## 1.811923 1.816606 1.003660# Kondisi numerik (opsional): kappa > 30 indikasi masalah
kappa(model.matrix(mod0))## [1] 2.263679Rule 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.06861Jika 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-12qqnorm(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 04.5 Uji spesifikasi (linearitas) – Ramsey RESET
lmtest::resettest(mod0)##
## RESET test
##
## data: mod0
## RESET = 4.2735, df1 = 2, df2 = 194, p-value = 0.015274.6 Plot Diagnostik Standar
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-16broom::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 NaNPilih 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 TRUEcar::outlierTest(mod0) # menguji outlier pada residual studentized dengan koreksi Bonferroni## rstudent unadjusted p-value Bonferroni p
## 162 6.450812 8.5823e-10 1.7165e-07
## 90 4.302679 2.6674e-05 5.3348e-03
## 18 4.188058 4.2583e-05 8.5167e-03Plot 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- 9broom::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-147. 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), multikolin- earitas (VIF), pengaruh (Cook’s D), normalitas (QQ-plot). Laporkan temuan & rekomendasi.
jawaban :
1. Heteroskedastisitas — Uji Breusch–Pagan
Inti Konsep
- Yang diuji: apakah varians residual konstan di
semua level prediksi (homoskedastisitas) atau tidak
(heteroskedastisitas).
- Hipotesis:
- H0: varians residual konstan.
- H1: varians residual bergantung pada satu atau lebih prediktor (atau fitted values).
- H0: varians residual konstan.
Ide Singkat Cara Kerja
- Hitung residual kuadrat \( _i^2 \), lalu regresikan ke prediktor
atau nilai fitted \( \).
- Statistik LM mengikuti distribusi \( ^2 \) dengan df = jumlah
prediktor (tanpa intercept).
- Kaidah keputusan: p-value kecil → tolak H0 (ada heteroskedastisitas).
# Model regresi
data(mtcars)
model <- lm(mpg ~ wt + hp + disp, data = mtcars)
# Uji Breusch–Pagan
library(lmtest)
bptest(model)##
## studentized Breusch-Pagan test
##
## data: model
## BP = 0.9459, df = 3, p-value = 0.81432. Multikolinearitas — VIF (Variance Inflation Factor)
Inti Konsep
Yang diuji: apakah prediktor saling berkorelasi kuat sehingga membuat estimasi koefisien tidak stabil.
Rumus VIF untuk prediktor \(X_j\):
\[ _j = \]
dengan \(R_j^2\) adalah koefisien determinasi saat \(X_j\) diregresikan pada prediktor lain.
Patokan
- VIF ≈ 1 → aman
- VIF 5–10 → waspada
- VIF > 10 → masalah serius
Implementasi di R
# Model regresi
data(mtcars)
model <- lm(mpg ~ wt + hp + disp, data = mtcars)
# Hitung VIF
library(car)
vif(model)## wt hp disp
## 4.844618 2.736633 7.3245173. Observasi Berpengaruh — Cook’s Distance
Inti Konsep
- Tujuan: menemukan observasi yang, bila dihapus,
akan mengubah koefisien/fit model secara nyata.
- Cook’s Distance menggabungkan:
- Leverage: jarak titik dari pusat variabel prediktor
(X).
- Residual: seberapa jauh nilai aktual dari prediksi.
- Leverage: jarak titik dari pusat variabel prediktor
(X).
Patokan Umum
- Ambang awal: \( 4/n \).
- Untuk \( n = 32 \): \( 4/32 = 0.125 \).
- Untuk \( n = 32 \): \( 4/32 = 0.125 \).
- Nilai > 1 → hampir pasti sangat berpengaruh.
Implementasi di R
data(mtcars)
model <- lm(mpg ~ wt + hp + disp, data = mtcars)
# Hitung Cook's Distance
cooksD <- cooks.distance(model)
# Plot
plot(cooksD, type = "h", main="Cook's Distance", ylab="Cook's D")
abline(h = 4/length(cooksD), col="red", lty=2)
# Tampilkan observasi dengan Cook's D terbesar
sort(cooksD, decreasing = TRUE)[1:10]## Maserati Bora Chrysler Imperial Toyota Corolla Fiat 128
## 0.34029114 0.31997067 0.15297714 0.11960192
## Pontiac Firebird Lotus Europa AMC Javelin Dodge Challenger
## 0.06980693 0.05959750 0.04909944 0.04218196
## Merc 280C Toyota Corona
## 0.02606104 0.022158654. Normalitas Residual — QQ-plot
Inti Konsep
Tujuan: mengecek apakah residual mendekati distribusi normal.
Hal ini penting untuk inferensi (uji t, uji F, confidence interval).
QQ-plot:
- Titik residual sebaiknya mengikuti garis 45°.
- Deviasi di ekor → indikasi skewness atau heavy tails.
- Titik residual sebaiknya mengikuti garis 45°.
Implementasi di R
data(mtcars)
model <- lm(mpg ~ wt + hp + disp, data = mtcars)
# QQ-plot residual
qqnorm(rstudent(model), main="QQ-Plot Residuals")
qqline(rstudent(model), col="red")
Ringkasan Temuan
Model yang dianalisis:
\[ mpg wt + hp + disp \]
1.Heteroskedastisitas
- Tidak terdeteksi (Breusch–Pagan p ≈ 0.81–0.84).
- Residual dapat dianggap homoskedastis.
2.Multikolinearitas
- Terdapat multikolinearitas moderate, terutama pada
variabel
disp(VIF ≈ 7.3).
wtjuga relatif tinggi (≈ 4.8), namun masih di bawah ambang serius.
3.Observasi Berpengaruh
- Beberapa observasi memiliki Cook’s Distance di atas ambang
\(4/n\):
- Maserati Bora, Chrysler Imperial,
Toyota Corolla.
- Nilainya tidak ekstrem (>1), tetapi perlu review.
4.Normalitas Residual
- Residual cukup normal.
- Deviasi kecil di ekor masih wajar untuk n = 32.
5. Rekomendasi Praktis
- Model Alternatif
- Coba hilangkan salah satu dari
dispatauwt(pilih variabel yang paling informatif).
- Bandingkan hasil model dengan metrik AIC dan R² adjusted, serta perhatikan stabilitas koefisien.
- Coba hilangkan salah satu dari
- Analisis Sensitivitas
- Refit model tanpa observasi dengan Cook’s Distance besar (misal
Maserati Bora, Chrysler Imperial, Corolla).
- Lihat apakah kesimpulan berubah signifikan.
- Refit model tanpa observasi dengan Cook’s Distance besar (misal
Maserati Bora, Chrysler Imperial, Corolla).
- Inferensi Lebih Kuat
- Jika tujuan utama adalah inferensi (uji t, uji F, CI), gunakan juga robust standard errors (HC3) sebagai pembanding.
Catatan Akhir
Secara keseluruhan, model masih layak
digunakan.
Namun, kolinearitas moderat dan beberapa observasi berpengaruh
menunjukkan perlunya eksplorasi model alternatif agar hasil analisis
lebih stabil dan dapat dipercaya.
2. Robust SE (HC3) vs WLS pada data simulasi heteroskedastik; bandingkan standar error dan p-value.
Jawaban :
1) Setting Masalah (Ringkas)
Kita ingin mengestimasi model linier sederhana:
\[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i \]
dengan asumsi:
- \(\mathbb{E}[\varepsilon_i \mid x_i] =
0\)
- \(\mathrm{Var}(\varepsilon_i \mid x_i) = \sigma^2 \, \omega_i\)
Heteroskedastisitas
- Pada kasus homoskedastisitas, varians error konstan
(\(\omega_i = 1\)).
- Pada kasus heteroskedastisitas, varians error berubah mengikuti \(x_i\), misalnya:
\[ \mathrm{Var}(\varepsilon_i \mid x_i) \propto x_i^2 \]
Implikasi
- Estimator OLS (\(\hat\beta\)) tetap
unbiased.
- Tetapi standard error klasik menjadi salah sehingga
uji t/F dan p-value tidak valid.
- Solusi umum: gunakan robust standard errors (HC3) atau Weighted Least Squares (WLS) bila bobot varians diketahui.
2) OLS (Estimasi) dan Standard Error Klasik
Estimator OLS
Estimator OLS untuk model linier:
\[ \hat{\beta} = (X'X)^{-1}X'y \]
di mana: - \(X\) = matriks desain (n × p), - \(y\) = vektor respon, - \(\hat{\beta}\) = vektor koefisien hasil estimasi.
Varians & Standard Error Klasik
Dengan asumsi homoskedastisitas:
\[ \widehat{\mathrm{Var}}(\hat{\beta}) = \hat{\sigma}^2 (X'X)^{-1} \]
dengan:
\[ \hat{\sigma}^2 = \frac{\sum \hat{e}_i^2}{n - p} \]
dan \(\hat{e}_i = y_i - \hat{y}_i\) adalah residual.
Catatan Penting
- Rumus di atas hanya valid bila varian error
konstan.
- Pada data yang heteroskedastik, standard error ini
bias → uji t, uji F, dan p-value bisa
menyesatkan.
- Oleh karena itu diperlukan robust SE (HC3) atau metode lain seperti WLS.
3) Robust SE (HC3)
Ide Utama
- Koefisien OLS tidak berubah, hanya kovarians/SE
yang dikoreksi.
- HC3 menggunakan sandwich estimator yang memperhitungkan heteroskedastisitas dan leverage.
Rumus HC3
\[ \widehat{\mathrm{Var}}_{\text{HC3}}(\hat{\beta}) = (X'X)^{-1} \Big( X' \, \mathrm{diag}\!\left(\frac{\hat{e}_i^2}{(1-h_{ii})^2}\right) X \Big) (X'X)^{-1} \]
dengan: - \(\hat{e}_i\) = residual OLS, - \(h_{ii}\) = leverage (elemen diagonal matriks hat \(H = X(X'X)^{-1}X'\)).
Intuisi
- HC3 mengoreksi residual kuadrat dengan faktor \((1-h_{ii})^2\).
- Titik dengan leverage tinggi diberi penalti lebih besar.
- Hasil: SE yang lebih konservatif, p-value lebih realistis.
Implikasi Praktis
- Gunakan HC3 bila ada indikasi heteroskedastisitas.
- Cocok ketika bentuk varians error tidak
diketahui.
- Koefisien tetap sama dengan OLS, tapi uji t/F lebih valid.
4) Weighted Least Squares (WLS)
Ide Utama
- Bila kita mengetahui (atau bisa memodelkan) bentuk varians error
\(\mathrm{Var}(\varepsilon_i|x_i)\),
maka bisa digunakan WLS agar estimasi lebih efisien.
- Prinsip: beri bobot kecil pada observasi dengan varians error
besar,
dan bobot besar pada observasi dengan varians kecil.
Rumus WLS
\[ \hat{\beta}_{\text{WLS}} = (X'WX)^{-1} X'Wy \]
dengan: - \(W = \mathrm{diag}(w_i)\) matriks bobot, - biasanya \(w_i = \tfrac{1}{\hat{\sigma}_i^2}\), yaitu kebalikan dari varians error.
Contoh Kasus
- Jika \(\mathrm{Var}(\varepsilon_i) \propto
x_i^2\),
maka bobot yang tepat adalah \(w_i = \tfrac{1}{x_i^2}\).
Implikasi Praktis
- Koefisien bisa berubah dibanding OLS, karena kita
meminimalkan error berbobot.
- Jika bobot yang digunakan benar, WLS menghasilkan estimator yang:
- tak bias,
- efisien,
- SE lebih kecil dibanding OLS maupun OLS+HC3.
- tak bias,
- Tetapi bila bobot salah spesifikasi → hasil bisa menyesatkan.
5) Perbandingan pada Data Simulasi (Heteroskedastik)
Kita simulasikan data dengan varians error meningkat terhadap \(x\) sehingga terjadi heteroskedastisitas. Model kebenaran: \(y = \beta_0 + \beta_1 x + \varepsilon\) dengan \(\beta_0 = 2\), \(\beta_1 = 1.5\). Kita set \(\mathrm{sd}(\varepsilon_i) \propto x_i\) agar varians error membesar saat \(x\) besar.
set.seed(123)
n <- 40
x <- seq(0.5, 4, length.out = n)
sd <- 0.6 * x # sd error proporsional dengan x
eps <- rnorm(n, mean = 0, sd = sd)
y <- 2 + 1.5 * x + eps
dat <- data.frame(x, y, sd)
m_ols <- lm(y ~ x, data = dat)
m_wls <- lm(y ~ x, data = dat, weights = 1 / (x^2))
# SE klasik (default) untuk OLS
se_ols <- sqrt(diag(vcov(m_ols)))
# SE robust HC3
library(sandwich)
library(lmtest)
vc_hc3 <- vcovHC(m_ols, type = "HC3")
se_hc3 <- sqrt(diag(vc_hc3))
# SE WLS (default vcov dari lm dengan weights)
se_wls <- sqrt(diag(vcov(m_wls)))
# Ringkas koefisien, SE, p-value
tab <- data.frame(
Metode = c("OLS (klasik)", "OLS + HC3", "WLS"),
beta0 = c(coef(m_ols)[1], coef(m_ols)[1], coef(m_wls)[1]),
se_b0 = c(se_ols[1], se_hc3[1], se_wls[1]),
p_b0 = c(coeftest(m_ols)[1,4], coeftest(m_ols, vcov.=vc_hc3)[1,4], coeftest(m_wls)[1,4]),
beta1 = c(coef(m_ols)[2], coef(m_ols)[2], coef(m_wls)[2]),
se_b1 = c(se_ols[2], se_hc3[2], se_wls[2]),
p_b1 = c(coeftest(m_ols)[2,4], coeftest(m_ols, vcov.=vc_hc3)[2,4], coeftest(m_wls)[2,4])
)
knitr::kable(tab, digits = 5, caption = "Perbandingan koefisien, standar error, dan p-value")| Metode | beta0 | se_b0 | p_b0 | beta1 | se_b1 | p_b1 |
|---|---|---|---|---|---|---|
| OLS (klasik) | 1.82910 | 0.46313 | 0.00033 | 1.59916 | 0.18697 | 0 |
| OLS + HC3 | 1.82910 | 0.30483 | 0.00000 | 1.59916 | 0.15889 | 0 |
| WLS | 2.07486 | 0.20805 | 0.00000 | 1.48159 | 0.15306 | 0 |
plot(fitted(m_ols), resid(m_ols),
xlab = "Fitted (OLS)", ylab = "Residual (OLS)",
main = "Residual vs Fitted (OLS)")
abline(h = 0, lty = 2)
uji Breusch–Pagan untuk menegaskan heteroskedastisitas.
library(lmtest)
bptest(m_ols)##
## studentized Breusch-Pagan test
##
## data: m_ols
## BP = 4.3684, df = 1, p-value = 0.036613. Spline vs polinomial: modelkan mpg ~ wt dan bandingkan AIC/BIC + 5-fold CV.
jawaban :
- Setting Masalah
- Data:
mtcars(n = 32). - Variabel: respons =
mpg, prediktor =wt. - Tujuan: membandingkan model
polinomial dan spline untuk menangkap
kemungkinan hubungan non-linear antara
wtdanmpg.
Implementasi r
data(mtcars)
mtcars_small <- mtcars[, c("mpg","wt")]
str(mtcars_small)## 'data.frame': 32 obs. of 2 variables:
## $ mpg: num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...summary(mtcars_small)## mpg wt
## Min. :10.40 Min. :1.513
## 1st Qu.:15.43 1st Qu.:2.581
## Median :19.20 Median :3.325
## Mean :20.09 Mean :3.217
## 3rd Qu.:22.80 3rd Qu.:3.610
## Max. :33.90 Max. :5.424- Kandidat Model
Linear
\[
mpg = \beta_0 + \beta_1\, wt + \varepsilon
\]
Polinomial
- Derajat 2: \[
mpg = \beta_0 + \beta_1\,wt + \beta_2\,wt^2 + \varepsilon
\] - Derajat 3: \[
mpg = \beta_0 + \beta_1\,wt + \beta_2\,wt^2 + \beta_3\,wt^3 +
\varepsilon
\]
Natural spline
- df = 3
- df = 4
Catatan: spline memberi fleksibilitas lokal; polinomial memberi fleksibilitas global yang bisa berosilasi di tepi rentang data.
Implementasi r
library(splines)
# Pakai GLM gaussian agar mudah dipakai cv.glm
m_lin <- glm(mpg ~ wt, data = mtcars, family = gaussian())
m_p2 <- glm(mpg ~ poly(wt, 2, raw = TRUE), data = mtcars, family = gaussian())
m_p3 <- glm(mpg ~ poly(wt, 3, raw = TRUE), data = mtcars, family = gaussian())
m_ns3 <- glm(mpg ~ ns(wt, df = 3), data = mtcars, family = gaussian())
m_ns4 <- glm(mpg ~ ns(wt, df = 4), data = mtcars, family = gaussian())
mods <- list(
Linear = m_lin,
Poly2 = m_p2,
Poly3 = m_p3,
NS3 = m_ns3,
NS4 = m_ns4
)
lapply(mods, summary)## $Linear
##
## Call:
## glm(formula = mpg ~ wt, family = gaussian(), data = mtcars)
##
## 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
##
## (Dispersion parameter for gaussian family taken to be 9.277398)
##
## Null deviance: 1126.05 on 31 degrees of freedom
## Residual deviance: 278.32 on 30 degrees of freedom
## AIC: 166.03
##
## Number of Fisher Scoring iterations: 2
##
##
## $Poly2
##
## Call:
## glm(formula = mpg ~ poly(wt, 2, raw = TRUE), family = gaussian(),
## data = mtcars)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.9308 4.2113 11.856 1.21e-12 ***
## poly(wt, 2, raw = TRUE)1 -13.3803 2.5140 -5.322 1.04e-05 ***
## poly(wt, 2, raw = TRUE)2 1.1711 0.3594 3.258 0.00286 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 7.025705)
##
## Null deviance: 1126.05 on 31 degrees of freedom
## Residual deviance: 203.75 on 29 degrees of freedom
## AIC: 158.05
##
## Number of Fisher Scoring iterations: 2
##
##
## $Poly3
##
## Call:
## glm(formula = mpg ~ poly(wt, 3, raw = TRUE), family = gaussian(),
## data = mtcars)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.40370 15.58379 3.106 0.00431 **
## poly(wt, 3, raw = TRUE)1 -11.82598 15.46346 -0.765 0.45081
## poly(wt, 3, raw = TRUE)2 0.68938 4.74034 0.145 0.88541
## poly(wt, 3, raw = TRUE)3 0.04594 0.45070 0.102 0.91954
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 7.273924)
##
## Null deviance: 1126.05 on 31 degrees of freedom
## Residual deviance: 203.67 on 28 degrees of freedom
## AIC: 160.04
##
## Number of Fisher Scoring iterations: 2
##
##
## $NS3
##
## Call:
## glm(formula = mpg ~ ns(wt, df = 3), family = gaussian(), data = mtcars)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.286 1.664 19.406 < 2e-16 ***
## ns(wt, df = 3)1 -13.903 1.938 -7.173 8.32e-08 ***
## ns(wt, df = 3)2 -27.887 3.846 -7.251 6.80e-08 ***
## ns(wt, df = 3)3 -15.627 1.815 -8.612 2.34e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 7.231457)
##
## Null deviance: 1126.05 on 31 degrees of freedom
## Residual deviance: 202.48 on 28 degrees of freedom
## AIC: 159.85
##
## Number of Fisher Scoring iterations: 2
##
##
## $NS4
##
## Call:
## glm(formula = mpg ~ ns(wt, df = 4), family = gaussian(), data = mtcars)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.543 1.834 17.745 < 2e-16 ***
## ns(wt, df = 4)1 -12.942 2.172 -5.960 2.35e-06 ***
## ns(wt, df = 4)2 -15.236 2.565 -5.941 2.47e-06 ***
## ns(wt, df = 4)3 -28.229 5.081 -5.556 6.87e-06 ***
## ns(wt, df = 4)4 -16.272 1.999 -8.139 9.65e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 7.457091)
##
## Null deviance: 1126.05 on 31 degrees of freedom
## Residual deviance: 201.34 on 27 degrees of freedom
## AIC: 161.67
##
## Number of Fisher Scoring iterations: 2- Kriteria AIC/BIC (Konsep)
- AIC: \(\text{AIC} =
-2\ell + 2k\)
- BIC: \(\text{BIC} = -2\ell + k\log n\)
dengan: - \(\ell\) = log-likelihood model, - \(k\) = jumlah parameter efektif, - \(n\) = jumlah observasi.
Semakin kecil AIC/BIC → semakin baik trade-off
fit vs kompleksitas.
Model lebih kompleks biasanya menurunkan \(-2\ell\) (fit membaik) tetapi meningkatkan
penalti (2k atau \(k\log n\)).
Implementasi r
library(dplyr)
metrics_ic <- tibble::tibble(
Model = names(mods),
AIC = sapply(mods, AIC) |> as.numeric(),
BIC = sapply(mods, BIC) |> as.numeric()
) %>% arrange(AIC)
knitr::kable(metrics_ic, digits = 3, caption = "Perbandingan AIC/BIC")| Model | AIC | BIC |
|---|---|---|
| Poly2 | 158.048 | 163.911 |
| NS3 | 159.849 | 167.178 |
| Poly3 | 160.037 | 167.365 |
| NS4 | 161.669 | 170.463 |
| Linear | 166.029 | 170.427 |
- 5-fold Cross-Validation (Konsep)
- Bagi data menjadi 5 lipatan.
- Latih pada 4 lipatan, uji pada 1 lipatan, hitung MSE.
- Ulangi sampai setiap lipatan menjadi data uji, lalu rata-ratakan →
CV-MSE.
- Pilih model dengan CV-MSE terkecil untuk kinerja prediktif out-of-sample terbaik.
cv_mse <- function(data, fit, K = 5) {
# cv.glm: delta[1] = CV estimate (MSE untuk gaussian)
boot::cv.glm(data = data, glmfit = fit, K = K)$delta[1]
}
metrics_cv <- tibble::tibble(
Model = names(mods),
CV_MSE = sapply(mods, function(fit) cv_mse(mtcars, fit, K = 5))
) %>% arrange(CV_MSE)
knitr::kable(metrics_cv, digits = 4, caption = "5-fold CV (MSE)")| Model | CV_MSE |
|---|---|
| Poly3 | 7.8588 |
| Poly2 | 7.9531 |
| NS4 | 8.9777 |
| Linear | 9.7610 |
| NS3 | 9.9831 |
- Analisis Manual (Tanpa Hitung Angka)
Linear: kemungkinan underfit; hubungan
mpg–wtbiasanya tidak murni linear (penurunanmpgmakin tajam pada mobil lebih berat).Polinomial derajat 2: cukup fleksibel untuk lengkungan halus → AIC biasanya turun nyata dibanding linear.
Polinomial derajat 3: lebih fleksibel, tapi berisiko overfit; BIC sering lebih ketat dan bisa tidak mendukung kenaikan derajat.
Natural spline df=3/4: fleksibel dan stabil di tepi rentang; sering kompetitif atau lebih baik dari polinomial orde tinggi.
Prediksi manual urutan performa:
AIC: Linear > Poly-2 ≈ NS(df=3) < Poly-3 ≈ NS(df=4) \ (lebih kecil lebih baik; model lebih fleksibel cenderung menang, tetapi tergantung penalti)
BIC: cenderung memilih model lebih sederhana → Poly-2 atau NS(df=3) sering unggul.
5-fold CV: Poly-2 dan NS(df=3) biasanya memberi CV-MSE paling kecil; Linear underfit, sementara Poly-3/NS(df=4) tidak selalu memberi perbaikan stabil pada n kecil (32).
library(ggplot2)
wt_grid <- data.frame(wt = seq(min(mtcars$wt), max(mtcars$wt), length.out = 200))
pred_df <- dplyr::bind_rows(
data.frame(wt = wt_grid$wt,
mpg_hat = predict(m_lin, newdata = wt_grid),
Model = "Linear"),
data.frame(wt = wt_grid$wt,
mpg_hat = predict(m_p2, newdata = wt_grid),
Model = "Poly-2"),
data.frame(wt = wt_grid$wt,
mpg_hat = predict(m_p3, newdata = wt_grid),
Model = "Poly-3"),
data.frame(wt = wt_grid$wt,
mpg_hat = predict(m_ns3, newdata = wt_grid),
Model = "NS df=3"),
data.frame(wt = wt_grid$wt,
mpg_hat = predict(m_ns4, newdata = wt_grid),
Model = "NS df=4")
)
ggplot(mtcars, aes(wt, mpg)) +
geom_point() +
geom_line(data = pred_df, aes(wt, mpg_hat, color = Model), linewidth = 1) +
labs(title = "mpg ~ wt: Spline vs Polinomial",
x = "wt", y = "mpg")
6. Rekomendasi
- Linear: underfit (AIC/BIC besar, CV error
besar).
- Poly-2: keseimbangan bagus antara fit dan
kompleksitas.
- Poly-3: bisa menurunkan AIC, tapi BIC/CV belum
tentu membaik karena penalti & risiko overfit.
- NS(df=3): alternatif kuat setara Poly-2, sering
stabil dan halus.
- NS(df=4): sedikit lebih fleksibel, tapi penalti BIC lebih besar; CV sering tidak jauh lebih baik dari df=3.
Pilihan praktis:
- Untuk interpretasi sederhana: Polinomial derajat
2.
- Untuk fleksibilitas prediksi yang halus: 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 :
1) Pemodelan & Tujuan
Model mean (data bulanan):
\[ y_t = \beta_0 + \beta_1\,x_t + \varepsilon_t,\quad x_t=\sin\!\Big(\tfrac{2\pi t}{12}\Big),\quad t=1,\dots,n. \]
Struktur error AR(1):
\[ \varepsilon_t = \rho\,\varepsilon_{t-1} + u_t,\quad u_t \sim \text{i.i.d. }(0,\sigma_u^2),\ |\rho|<1. \]
Tujuan
1) Estimasi efek musiman \(\beta_1\)
dengan GLS yang memperhitungkan autokorelasi.
2) Menafsirkan \(\hat\rho\) sebagai
ukuran persistensi gangguan antarbulan.
3) Membandingkan standar error OLS vs GLS untuk \(\beta_1\).
2) Baseline OLS (mengabaikan AR)
Estimator OLS:
\[ \hat\beta^{\text{OLS}} = (X'X)^{-1}X'y,\quad X=[\mathbf{1},\,x]. \]
Jika \(x_t\) eksogen, \(\hat\beta^{\text{OLS}}\) tak bias.
Masalahnya, autokorelasi AR(1) membuat kovarian OLS klasik tidak valid
untuk uji t/F.
Kovarian error “benar” di bawah AR(1):
\[ \Omega=\frac{\sigma_u^2}{1-\rho^2}\big[\rho^{|t-s|}\big]_{t,s=1}^n. \]
3) GLS untuk AR(1): Ide & Transformasi
Secara teoretis:
\[ \hat\beta^{\text{GLS}}=(X'\Omega^{-1}X)^{-1}X'\Omega^{-1}y. \]
Praktik umum menggunakan FGLS via transformasi Prais–Winsten:
Untuk \(t\ge 2\): \[ y_t^*=y_t-\rho\,y_{t-1},\qquad x_t^*=x_t-\rho\,x_{t-1}. \]
Untuk \(t=1\): \[ y_1^*=\sqrt{1-\rho^2}\,y_1,\qquad x_1^*=\sqrt{1-\rho^2}\,x_1. \]
Kemudian OLS pada data bertanda bintang:
\[ \hat\beta^{\text{GLS}}=\big((X^*)'X^*\big)^{-1}(X^*)'y^*,\qquad \widehat{\mathrm{Var}}(\hat\beta^{\text{GLS}})=\hat\sigma_*^2\big((X^*)'X^*\big)^{-1}. \]
4) Estimasi \(\rho\) dan Interpretasi \(\hat\rho\)
Langkah ringkas:
- Estimasi OLS \(y_t=\beta_0+\beta_1
x_t+\varepsilon_t\), simpan residual \(\hat e_t\).
- Regress \(\hat e_t\) pada \(\hat e_{t-1}\) untuk \(\hat\rho\): \[ \hat\rho=\frac{\sum_{t=2}^n \hat e_t \hat e_{t-1}}{\sum_{t=2}^n \hat e_{t-1}^2}. \]
- Lakukan transformasi Prais–Winsten memakai \(\hat\rho\), lalu OLS pada \((y^*,X^*)\).
- (Opsional) Ulangi hingga konvergen.**
Makna \(\hat\rho\):
- \(\hat\rho\approx 0\): hampir tidak
ada autokorelasi → OLS dan GLS mirip.
- \(\hat\rho>0\): gangguan
persisten. Contoh \(\hat\rho=0.7\)
berarti ~70% shock bertahan ke bulan berikutnya.
- \(\hat\rho<0\): pola residual berselang-seling.
Catatan musiman sinus:
Untuk \(x_t=\sin(2\pi t/12)\), korelasi
lag-1 dari sinyalnya \(\cos(2\pi/12)=\cos(\pi/6)\approx 0.866\).
Jika error juga persisten (\(\hat\rho>0\)), koreksi AR(1) penting
agar inferensi tidak terlalu optimistis.
5) SE OLS vs SE GLS: Apa yang Diharapkan
- OLS (SE klasik) mengasumsikan i.i.d. error. Saat
\(\hat\rho>0\), SE klasik cenderung
terlalu kecil sehingga p-value terlihat terlalu
kecil.
- GLS memperhitungkan autokorelasi → SE untuk \(\beta_1\) biasanya lebih
besar dan lebih realistis.
- Dibanding SE robust HAC/Newey–West, GLS
“tepat-spesifikasi” untuk AR(1) sehingga efisien bila model AR(1)
benar.
- Intuisi: makin besar \(\hat\rho>0\), efektif sampel menurun, SE “benar” meningkat relatif OLS klasik.
6) Apa yang Perlu Dilaporkan
- Spesifikasi model: \(y_t=\beta_0+\beta_1\sin(2\pi
t/12)+\varepsilon_t\), \(\varepsilon_t\) AR(1).
- Nilai \(\hat\rho\) dari residual
OLS dan tafsirannya.
- Estimasi \(\hat\beta_1\) OLS dan
GLS.
- Perbandingan SE dan p-value: OLS
(klasik) vs GLS (Prais–Winsten).
- Opsional: sertakan SE robust (HAC) sebagai pembanding.
- Opsional: sertakan SE robust (HAC) sebagai pembanding.
- Catatan: jika \(\hat\rho>0\) dan SE GLS > SE OLS, jelaskan bahwa OLS mengabaikan autokorelasi sehingga meng-understate ketidakpastian.
7) Checklist Langkah
- OLS → dapatkan residual \(\hat
e_t\).
- Estimasi \(\hat\rho\) via regresi
\(\hat e_t\) pada \(\hat e_{t-1}\).
- Transformasi Prais–Winsten dengan \(\hat\rho\) untuk \(y_t\) dan \(x_t\).
- OLS pada \((y^*,X^*)\) → \(\hat\beta^{\text{GLS}}\), SE(GLS).
1.14.2 Tugas
• Esai: Mengapa OLS tetap unbiased di bawah heteroskedastisitas, namun tidak efisien? Bagaimana robust SE mengatasi masalah inferensi?
Jawaban
1. Mengapa OLS tetap unbiased?
Kita mulai dari model linier:
\[ y = X\beta + \varepsilon,\qquad \mathbb{E}[\varepsilon \mid X] = 0,\quad \mathrm{Var}(\varepsilon \mid X)=\Omega \]
dengan \(\Omega\) tidak harus \(\sigma^2 I\) (boleh heteroskedastik).
Estimator OLS adalah:
\[ \hat{\beta}_{\text{OLS}} = (X'X)^{-1}X'y. \]
Ambil harapan bersyarat pada \(X\):
\[ \begin{aligned} \mathbb{E}[\hat{\beta}_{\text{OLS}} \mid X] &= \mathbb{E}\!\left[(X'X)^{-1}X'(X\beta + \varepsilon)\,\middle|\,X\right] \\ &= (X'X)^{-1}X'X\,\beta \;+\; (X'X)^{-1}X'\,\mathbb{E}[\varepsilon \mid X] \\ &= \beta \;+\; (X'X)^{-1}X'\,0 \\ &= \beta. \end{aligned} \]
Kesimpulan: selama asumsi eksogenitas berlaku (\(\mathbb{E}[\varepsilon \mid X]=0\)), OLS tetap tak bias meski ada heteroskedastisitas. Heteroskedastisitas memengaruhi varians \(\hat{\beta}\), bukan rata-ratanya.
Catatan penting - Heteroskedastisitas tidak menggeser nilai tengah \(\hat{\beta}\), tetapi membuat \(\mathrm{Var}(\hat{\beta}_{\text{OLS}} \mid X) = (X'X)^{-1}X'\Omega X(X'X)^{-1}\) bukan \(\sigma^2(X'X)^{-1}\). - Akibatnya OLS tidak efisien (bukan lagi BLUE) dan SE klasik jadi tidak valid untuk inferensi.
Implementasi R
1.1) Simulasi: OLS tetap tak-bias (eksogenitas terpenuhi)
Model kebenaran:
\(y = \beta_0 + \beta_1 x +
\varepsilon\), dengan \(E[\varepsilon
\mid x] = 0\) dan \(sd(\varepsilon) =
0.3 + 0.7x\).
beta0 <- 1
beta1 <- 2
nsim <- 1000 # jumlah replikasi simulasi
n <- 300 # ukuran sampel per simulasi
beta1_hat <- numeric(nsim)
for (s in 1:nsim) {
x <- runif(n, 0, 3)
sd <- 0.3 + 0.7*x # heteroskedastisitas: sd naik saat x naik
eps <- rnorm(n, 0, sd)
y <- beta0 + beta1*x + eps
fit <- lm(y ~ x)
beta1_hat[s] <- coef(fit)["x"]
}
c(mean_beta1_hat = mean(beta1_hat),
true_beta1 = beta1,
bias = mean(beta1_hat) - beta1)## mean_beta1_hat true_beta1 bias
## 2.001085494 2.000000000 0.0010854941.2) Mengapa inferensi OLS klasik bisa salah? (coverage SE klasik vs robust)
Kita bandingkan cakupan 95% untuk \(\beta_1\) menggunakan:
- SE klasik (asumsi homoskedastisitas)
- SE robust HC3 (konsisten di bawah heteroskedastisitas)
Implementasi R
library(sandwich)
library(lmtest)
nsim <- 1000
n <- 300
cover_classic <- logical(nsim)
cover_hc3 <- logical(nsim)
sebar_classic <- numeric(nsim)
sebar_hc3 <- numeric(nsim)
for (s in 1:nsim) {
x <- runif(n, 0, 3)
sd <- 0.3 + 0.7*x
eps <- rnorm(n, 0, sd)
y <- beta0 + beta1*x + eps
fit <- lm(y ~ x)
# SE klasik
se_c <- sqrt(diag(vcov(fit)))["x"]
ci_c <- coef(fit)["x"] + c(-1,1) * 1.96 * se_c
cover_classic[s] <- (beta1 >= ci_c[1] && beta1 <= ci_c[2])
sebar_classic[s] <- se_c
# SE robust HC3
vc_hc3 <- vcovHC(fit, type = "HC3")
se_r <- sqrt(diag(vc_hc3))["x"]
ci_r <- coef(fit)["x"] + c(-1,1) * 1.96 * se_r
cover_hc3[s] <- (beta1 >= ci_r[1] && beta1 <= ci_r[2])
sebar_hc3[s] <- se_r
}
results <- data.frame(
Metrik = c("Rata-rata SE", "Coverage 95%"),
Klasik = c(mean(sebar_classic), mean(cover_classic)),
HC3 = c(mean(sebar_hc3), mean(cover_hc3))
)
knitr::kable(results, digits = 4, caption = "SE klasik vs SE robust (HC3): rata-rata SE dan coverage 95%")| Metrik | Klasik | HC3 |
|---|---|---|
| Rata-rata SE | 0.0986 | 0.1051 |
| Coverage 95% | 0.9240 | 0.9440 |
Interpretasi yang diharapkan
Coverage klasik sering < 0.95 (terlalu optimistis) saat heteroskedastisitas ada.
Coverage HC3 mendekati 0.95 → inferensi valid.
Rata-rata SE HC3 biasanya lebih besar dari klasik (lebih “jujur”).
1.3) Visual satu realisasi (opsional)
# Satu contoh data untuk visual
x <- runif(n, 0, 3)
sd <- 0.3 + 0.7*x
eps <- rnorm(n, 0, sd)
y <- beta0 + beta1*x + eps
fit <- lm(y ~ x)
par(mfrow = c(1,2))
plot(x, y, pch = 19, cex = 0.6,
main = "Scatter y vs x (heteroskedastik)")
abline(fit, col = "red", lwd = 2)
plot(fitted(fit), resid(fit), pch = 19, cex = 0.6,
main = "Residual vs Fitted (OLS)",
xlab = "Fitted", ylab = "Residual")
abline(h = 0, lty = 2)
par(mfrow = c(1,1))Ringkas:
Simulasi menunjukkan OLS tak-bias (rata-rata \(\hat{\beta}_1 \approx \beta_1\)).
SE klasik gagal menjaga coverage ketika ada heteroskedastisitas.
Robust SE (HC3) mengoreksi kovarians sehingga uji t/CI kembali valid tanpa mengubah koefisien OLS.
2) Mengapa tidak efisien?
Varians sebenarnya dari \(\hat{\beta}_{OLS}\) (bersyarat pada \(X\)) di bawah heteroskedastisitas:
\(Var(\hat{\beta}_{OLS} \mid X) = (X'X)^{-1} X' \, \Omega \, X \, (X'X)^{-1}.\)
Di bawah homoskedastisitas (\(\Omega =
\sigma^2 I\)) ini menjadi \(\sigma^2
(X'X)^{-1}\) dan OLS adalah BLUE (Gauss–Markov).
Di bawah heteroskedastisitas, bentuk di atas bukan \(\sigma^2 (X'X)^{-1}\) dan OLS bukan
lagi estimator linear tak-bias dengan varians minimum.
Intuisi: jika beberapa titik punya varians error besar, menyamakan bobot semua titik “mencampur” banyak noise. GLS/WLS memakai bobot \(w_i \propto 1/\sigma_i^2\) untuk menurunkan kontribusi observasi yang sangat berisik, sehingga varians koefisien lebih kecil (lebih efisien).
Implementasi R
Mendemonstrasikan bahwa:
OLS tetap tak-bias (eksogenitas terpenuhi), tetapi tidak efisien di bawah heteroskedastisitas.
WLS dengan bobot yang sesuai var(error) memiliki varians estimator koefisien lebih kecil (lebih efisien).
Desain simulasi
Model kebenaran:
\(y = \beta_0 + \beta_1 x + \varepsilon, \; E[\varepsilon \mid x] = 0, \; sd(\varepsilon) = 0.3 + 0.7x,\)
sehingga \(Var(\varepsilon) = (0.3 + 0.7x)^2\) (heteroskedastik meningkat saat \(x\) besar).
Estimator dibandingkan:
- OLS:
lm(y ~ x)(bobot seragam).
- WLS:
lm(y ~ x, weights = 1 / (0.3 + 0.7x)^2)(bobot “benar”).
Metrik:
- Rata-rata \(\hat{\beta}_1\)
(bias),
- Varians/SD empirik \(\hat{\beta}_1\),
- Rata-rata SE terestimasi (klasik untuk OLS dan WLS),
- Rasio efisiensi: \(Var(OLS)/Var(WLS) > 1 \;\;\Rightarrow\;\;\) WLS lebih efisien.
beta0 <- 1
beta1 <- 2
nsim <- 2000 # jumlah replikasi
n <- 300 # ukuran sampel per replikasi
b1_ols <- numeric(nsim)
b1_wls <- numeric(nsim)
se_ols <- numeric(nsim) # SE klasik OLS
se_wls <- numeric(nsim) # SE WLS (default vcov pada lm berbobot)
for (s in 1:nsim) {
x <- runif(n, 0, 3)
sd <- 0.3 + 0.7*x
eps <- rnorm(n, 0, sd)
y <- beta0 + beta1*x + eps
fit_ols <- lm(y ~ x)
fit_wls <- lm(y ~ x, weights = 1/(sd^2))
b1_ols[s] <- coef(fit_ols)["x"]
b1_wls[s] <- coef(fit_wls)["x"]
se_ols[s] <- coef(summary(fit_ols))["x","Std. Error"]
se_wls[s] <- coef(summary(fit_wls))["x","Std. Error"]
}
# Ringkasan kinerja
mean_ols <- mean(b1_ols); mean_wls <- mean(b1_wls)
bias_ols <- mean_ols - beta1; bias_wls <- mean_wls - beta1
sd_ols <- sd(b1_ols); sd_wls <- sd(b1_wls)
var_ratio <- (sd_ols^2)/(sd_wls^2)
tab <- data.frame(
Metrik = c("Rata-rata beta1_hat", "Bias", "SD empirik beta1_hat",
"Rata-rata SE (terlapor)", "Rasio Var(OLS)/Var(WLS)"),
OLS = c(mean_ols, bias_ols, sd_ols, mean(se_ols), var_ratio),
WLS = c(mean_wls, bias_wls, sd_wls, mean(se_wls), NA)
)
knitr::kable(tab, digits = 5, caption = "Perbandingan OLS vs WLS di bawah heteroskedastisitas")| Metrik | OLS | WLS |
|---|---|---|
| Rata-rata beta1_hat | 1.99904 | 1.99675 |
| Bias | -0.00096 | -0.00325 |
| SD empirik beta1_hat | 0.10923 | 0.07607 |
| Rata-rata SE (terlapor) | 0.09870 | 0.07470 |
| Rasio Var(OLS)/Var(WLS) | 2.06191 | NA |
Interpretasi yang diharapkan:
- Bias keduanya ≈ 0 (eksogenitas → tak-bias).
- SD/Varians empirik WLS lebih kecil dari OLS → WLS lebih
efisien.
- Rasio Var(OLS)/Var(WLS) > 1 mengonfirmasi efisiensi WLS.
- Rata-rata SE yang dilaporkan WLS biasanya < OLS (konsisten dengan efisiensi).
Visual distribusi estimator (opsional)
library(ggplot2)
df_plot <- rbind(
data.frame(beta1_hat = b1_ols, Metode = "OLS"),
data.frame(beta1_hat = b1_wls, Metode = "WLS")
)
ggplot(df_plot, aes(beta1_hat, after_stat(density), color = Metode)) +
geom_density(linewidth = 1) +
geom_vline(xintercept = beta1, linetype = 2) +
labs(title = "Distribusi Sampling \\u03B2\u2081 Hat: OLS vs WLS",
x = expression(hat(beta)[1]), y = "Kerapatan") +
theme_minimal()
Catatan penting
Hasil ini bergantung pada bobot yang benar. Jika bobot salah spesifikasi, WLS bisa tidak efisien atau bahkan bias (jika ada endogenitas bobot).
OLS + robust SE (HC3) memperbaiki inferensi (SE/p-value valid) tanpa membuat estimator lebih efisien. Jika tujuanmu efisiensi (varians kecil) dan kamu tahu bentuk var(error), gunakan WLS/GLS.
3) Mengapa inferensi OLS klasik jadi salah?
Rumus SE “klasik” mengasumsikan \(\Omega = \sigma^2 I\):
\(\widehat{Var}_{klasik}(\hat{\beta}) = \hat{\sigma}^2 (X'X)^{-1}, \quad \hat{\sigma}^2 = \tfrac{1}{n-k} \sum \hat{e}_i^2.\)
Jika \(\Omega\) heteroskedastik,
estimator ini tidak konsisten: sering meng-underestimate
ketidakpastian.
Akibatnya, uji t/F dan p-value bisa terlalu optimistis (sering “terlihat
signifikan” padahal tidak).
Implementasi R
Tujuan
- Menunjukkan bahwa SE/p-value klasik (asumsi homoskedastisitas) bisa
menyesatkan saat error heteroskedastik.
- Menunjukkan bahwa robust SE (HC3) memberi inferensi yang lebih valid (ukuran uji mendekati level nominal).
Desain simulasi (uji ukuran/Type I error)
Kita uji hipotesis \(H_0 : \beta_1 = 0\) saat benar \(\beta_1 = 0\) dengan error heteroskedastik:
\(y = \beta_0 + \beta_1 x + \varepsilon, \; \beta_0 = 1, \; \beta_1 = 0, \; sd(\varepsilon) = 0.3 + 0.7x.\)
Jika inferensi benar, tingkat penolakan pada \(\alpha = 0.05\) harus \(\approx 0.05\).
Kita bandingkan p-value klasik vs p-value robust HC3.
library(sandwich)
library(lmtest)
nsim <- 2000 # replikasi
n <- 300 # sampel per replikasi
alpha <- 0.05
beta0 <- 1
beta1 <- 0 # Null benar
rej_classic <- logical(nsim)
rej_hc3 <- logical(nsim)
avg_se_classic <- numeric(nsim)
avg_se_hc3 <- numeric(nsim)
for (s in 1:nsim) {
x <- runif(n, 0, 3)
sd <- 0.3 + 0.7*x # heteroskedastik (var naik saat x naik)
eps <- rnorm(n, 0, sd)
y <- beta0 + beta1*x + eps
fit <- lm(y ~ x)
# p-value OLS klasik
p_classic <- coef(summary(fit))["x", "Pr(>|t|)"]
rej_classic[s] <- (p_classic < alpha)
avg_se_classic[s] <- coef(summary(fit))["x", "Std. Error"]
# p-value robust (HC3)
vc_hc3 <- vcovHC(fit, type = "HC3")
test_hc3 <- coeftest(fit, vcov. = vc_hc3)
p_hc3 <- test_hc3["x", "Pr(>|t|)"]
rej_hc3[s] <- (p_hc3 < alpha)
avg_se_hc3[s] <- sqrt(diag(vc_hc3))["x"]
}
# Ringkasan ukuran uji (Type I error) dan rata-rata SE
hasil <- data.frame(
Metrik = c("Empirical size Pr(reject|H0)", "Rata-rata SE slope"),
OLS_klasik = c(mean(rej_classic), mean(avg_se_classic)),
Robust_HC3 = c(mean(rej_hc3), mean(avg_se_hc3))
)
knitr::kable(hasil, digits = 4, caption = "Uji ukuran & rata-rata SE: OLS klasik vs Robust HC3")| Metrik | OLS_klasik | Robust_HC3 |
|---|---|---|
| Empirical size Pr(reject|H0) | 0.0695 | 0.0570 |
| Rata-rata SE slope | 0.0987 | 0.1054 |
Interpretasi yang diharapkan:
- Empirical size OLS klasik sering > 0.05 (uji “terlalu sering
menolak” → optimistis).
- Empirical size HC3 mendekati 0.05 (uji valid).
- Rata-rata SE HC3 umumnya lebih besar dari klasik → ketidakpastian diukur lebih “jujur”.
Visual distribusi p-value (opsional)
# Sampel p-value representatif dari loop terakhir (atau jalankan ulang khusus untuk plotting)
# Di sini kita simulasikan ulang sekali untuk mem-plot p-value-nya
B <- 2000
p_class <- p_hc3_vec <- numeric(B)
for (s in 1:B) {
x <- runif(n, 0, 3)
sd <- 0.3 + 0.7*x
eps <- rnorm(n, 0, sd)
y <- beta0 + beta1*x + eps
fit <- lm(y ~ x)
p_class[s] <- coef(summary(fit))["x","Pr(>|t|)"]
p_hc3_vec[s] <- coeftest(fit, vcov.=vcovHC(fit, type="HC3"))["x","Pr(>|t|)"]
}
par(mfrow=c(1,2))
hist(p_class, breaks=20, main="P-value OLS klasik", xlab="p", col="grey")
abline(v=alpha, col=2, lty=2)
hist(p_hc3_vec, breaks=20, main="P-value Robust HC3", xlab="p", col="grey")
abline(v=alpha, col=2, lty=2)
par(mfrow=c(1,1))Contoh satu realisasi (perbandingan SE & p-value)
x <- runif(n, 0, 3)
sd <- 0.3 + 0.7*x
eps <- rnorm(n, 0, sd)
y <- beta0 + beta1*x + eps
fit <- lm(y ~ x)
# Klasik
se_c <- coef(summary(fit))["x","Std. Error"]
p_c <- coef(summary(fit))["x","Pr(>|t|)"]
# Robust HC3
vc_hc3 <- vcovHC(fit, type="HC3")
se_r <- sqrt(diag(vc_hc3))["x"]
p_r <- coeftest(fit, vcov.=vc_hc3)["x","Pr(>|t|)"]
data.frame(
Metode = c("OLS klasik", "Robust HC3"),
SE = c(se_c, se_r),
p_val = c(p_c, p_r)
)## Metode SE p_val
## OLS klasik 0.1014053 0.6257025
## x Robust HC3 0.1164521 0.6709984Kesimpulan singkat
Di bawah heteroskedastisitas, SE klasik cenderung underestimate, meningkatkan false positive (Type I error > 5%).
Robust SE (HC3) memperbaiki kovarians (tanpa mengubah koefisien OLS), sehingga uji t/CI menjadi valid (empirical size ≈ 5%).
4) Bagaimana robust SE memperbaiki inferensi?
Solusi praktis: pakai heteroskedasticity-consistent (HC) covariance (“sandwich” White/Huber). Versi dasar (HC0):
\(\widehat{Var}_{HC0}(\hat{\beta}) = (X'X)^{-1} (X' \, diag(\hat{e}_i^2) \, X)(X'X)^{-1}.\)
Variasi populer:
- HC1: skala \(\tfrac{n}{n-k}\)
(koreksi derajat bebas).
- HC2/HC3: koreksi leverage, mis. HC3 memakai \(\tfrac{\hat{e}_i^2}{(1-h_{ii})^2}\).
Yang berubah hanya SE/kovariansnya, bukan \(\hat{\beta}\). Dengan robust SE, statistik t, p-value, dan CI jadi konsisten terhadap heteroskedastisitas (untuk \(n\) cukup besar), sehingga keputusan inferensi kembali bisa dipercaya.
Ringkas:
- OLS + robust SE → inferensi valid, tapi efisiensi tidak otomatis
optimal.
- WLS/GLS → jika bentuk \(\Omega\) “diketahui/termodelkan” dengan baik, bisa lebih efisien daripada OLS-robust.
Implementasi R
Contoh satu realisasi: SE/p-value klasik vs HC0–HC3
library(sandwich)
library(lmtest)
library(dplyr)
# Satu dataset heteroskedastik (true beta1 = 2)
n <- 300
beta0 <- 1
beta1 <- 2
x <- runif(n, 0, 3)
sd <- 0.3 + 0.7*x
eps <- rnorm(n, 0, sd)
y <- beta0 + beta1*x + eps
fit <- lm(y ~ x)
# Kovarians & SE
vc_classic <- vcov(fit) # klasik (homoskedastis)
vc_hc0 <- vcovHC(fit, type = "HC0")
vc_hc1 <- vcovHC(fit, type = "HC1")
vc_hc2 <- vcovHC(fit, type = "HC2")
vc_hc3 <- vcovHC(fit, type = "HC3")
get_row <- function(vcovM, label) {
se <- sqrt(diag(vcovM))["x"]
bhat <- coef(fit)["x"]
tval <- bhat / se
pval <- 2*pt(-abs(tval), df = n - 2)
ci <- bhat + c(-1,1) * 1.96 * se
data.frame(Metode = label, beta1_hat = bhat, SE = se, t = tval, p = pval,
CI_low = ci[1], CI_high = ci[2])
}
tab1 <- bind_rows(
get_row(vc_classic, "Klasik"),
get_row(vc_hc0, "HC0"),
get_row(vc_hc1, "HC1"),
get_row(vc_hc2, "HC2"),
get_row(vc_hc3, "HC3")
)
knitr::kable(tab1, digits = 5, caption = "Satu realisasi: koefisien sama, SE/p-value berbeda (Klasik vs HC0–HC3)")| Metode | beta1_hat | SE | t | p | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| x…1 | Klasik | 1.96431 | 0.10080 | 19.48705 | 0 | 1.76674 | 2.16188 |
| x…2 | HC0 | 1.96431 | 0.11010 | 17.84159 | 0 | 1.74852 | 2.18010 |
| x…3 | HC1 | 1.96431 | 0.11047 | 17.78202 | 0 | 1.74780 | 2.18082 |
| x…4 | HC2 | 1.96431 | 0.11065 | 17.75224 | 0 | 1.74743 | 2.18119 |
| x…5 | HC3 | 1.96431 | 0.11121 | 17.66326 | 0 | 1.74634 | 2.18228 |
Catatan cepat
- Koefisien OLS sama untuk semua baris (kita tidak mengubah
estimatornya).
- SE klasik biasanya lebih kecil → p-value lebih optimistis.
- HC (khususnya HC3) meningkatkan SE ke nilai yang lebih “jujur”.
Simulasi cakupan 95% CI: Klasik vs HC0–HC3
Kita cek coverage 95% CI untuk \(\beta_1 = 2\) di bawah heteroskedastisitas. CI yang valid seharusnya mengandung nilai benar sekitar 95% dari replikasi.
nsim <- 1500
n <- 300
beta0 <- 1
beta1 <- 2
cover <- matrix(FALSE, nrow = nsim, ncol = 5)
colnames(cover) <- c("Klasik","HC0","HC1","HC2","HC3")
avg_se <- matrix(NA_real_, nrow = nsim, ncol = 5)
colnames(avg_se) <- c("Klasik","HC0","HC1","HC2","HC3")
for (s in 1:nsim) {
x <- runif(n, 0, 3)
sd <- 0.3 + 0.7*x
eps <- rnorm(n, 0, sd)
y <- beta0 + beta1*x + eps
fit <- lm(y ~ x)
# Semua kovarians
vcs <- list(
Klasik = vcov(fit),
HC0 = vcovHC(fit, type="HC0"),
HC1 = vcovHC(fit, type="HC1"),
HC2 = vcovHC(fit, type="HC2"),
HC3 = vcovHC(fit, type="HC3")
)
b1 <- coef(fit)["x"]
j <- 1
for (nm in names(vcs)) {
se <- sqrt(diag(vcs[[nm]]))["x"]
ci <- b1 + c(-1,1)*1.96*se
cover[s, j] <- (beta1 >= ci[1] && beta1 <= ci[2])
avg_se[s, j] <- se
j <- j + 1
}
}
res <- data.frame(
Metrik = c("Coverage 95%", "Rata-rata SE"),
Klasik = c(mean(cover[, "Klasik"]), mean(avg_se[, "Klasik"])),
HC0 = c(mean(cover[, "HC0"]), mean(avg_se[, "HC0"])),
HC1 = c(mean(cover[, "HC1"]), mean(avg_se[, "HC1"])),
HC2 = c(mean(cover[, "HC2"]), mean(avg_se[, "HC2"])),
HC3 = c(mean(cover[, "HC3"]), mean(avg_se[, "HC3"]))
)
knitr::kable(res, digits = 4, caption = "Empirical coverage 95% dan rata-rata SE: Klasik vs HC0–HC3")| Metrik | Klasik | HC0 | HC1 | HC2 | HC3 |
|---|---|---|---|---|---|
| Coverage 95% | 0.9233 | 0.9393 | 0.9407 | 0.9407 | 0.9420 |
| Rata-rata SE | 0.0988 | 0.1044 | 0.1048 | 0.1049 | 0.1055 |
Interpretasi yang diharapkan
- Coverage klasik biasanya < 0.95 (undercoverage) → uji/CI terlalu
optimistis.
- HC (terutama HC3) mendekati 0.95 → inferensi valid di bawah
heteroskedastisitas.
- Rata-rata SE HC > SE klasik (menggambarkan ketidakpastian yang sebenarnya).
Ringkas
- Robust SE tidak mengubah koefisien, hanya kovariansnya.
- Di bawah heteroskedastisitas, robust SE (HC0–HC3) membuat uji t/CI
valid.
- HC3 sering direkomendasikan karena paling “ketat” (mengoreksi leverage), terutama untuk sampel kecil–menengah.
5) Kapan pakai yang mana?
- Tidak yakin bentuk heteroskedastisitas → laporkan robust SE (HC3)
secara default.
- Ada model varians yang kredibel (mis. ragam ∝ \(x_i^2\), atau dari diagnostik) →
pertimbangkan WLS/GLS untuk efisiensi.
- Jika data berkelompok/panel → pertimbangkan cluster-robust SE (korelasi intra-cluster).
- Esai Turunkan bentuk estimator WLS saat varian residual diketahui proporsional terhadap |𝑥|.
Jawaban :
estimator WLS untuk regresi sederhana \[ y_i=\beta_0+\beta_1 x_i+\varepsilon_i,\qquad \mathbb{E}[\varepsilon_i|x_i]=0,\quad \mathrm{Var}(\varepsilon_i|x_i)=\sigma^2\,|x_i|. \]
1) Bobot WLS dari struktur varians
Jika \(\mathrm{Var}(\varepsilon_i|x_i)=\sigma^2\,|x_i|\), maka bobot optimum \[ w_i=\frac{1}{\mathrm{Var}(\varepsilon_i|x_i)}=\frac{1}{\sigma^2|x_i|}\ \ \propto\ \ \frac{1}{|x_i|}. \] Faktor \(\sigma^{-2}\) tidak berpengaruh pada solusi, jadi cukup ambil \(w_i=1/|x_i|\).
2) Estimator WLS (bentuk matriks)
Tuliskan \(X=[\mathbf{1},\,x]\) (kolom intercept dan \(x\)), \(y\) vektor respons, dan \(W=\mathrm{diag}(w_i)\). Estimator WLS meminimalkan \[ Q(\beta)=(y-X\beta)'W(y-X\beta), \] turunannya nol memberi normal equation berbobot: \[ X'WX\,\hat\beta=X'Wy\quad\Rightarrow\quad \boxed{\ \hat\beta_{\text{WLS}}=(X'WX)^{-1}X'Wy\ }. \]
Dengan bobot \(w_i=1/|x_i|\), ini langsung menjadi bentuk yang diinginkan.
3) Estimator eksplisit untuk \(\beta_0,\beta_1\) (kasus sederhana)
Definisikan jumlah berbobot: \[ S_w=\sum_i w_i,\quad S_x=\sum_i w_i x_i,\quad S_y=\sum_i w_i y_i,\quad S_{xx}=\sum_i w_i x_i^2,\quad S_{xy}=\sum_i w_i x_i y_i. \] Maka \[ \boxed{\ \hat\beta_1=\frac{S_w S_{xy}-S_x S_y}{S_w S_{xx}-S_x^2}\ },\qquad \boxed{\ \hat\beta_0=\frac{S_y-\hat\beta_1 S_x}{S_w}\ }. \] Ini identik dengan bentuk “mean berbobot”. Jika \(\bar x_w=S_x/S_w\) dan \(\bar y_w=S_y/S_w\), maka \[ \hat\beta_1=\frac{\sum_i w_i(x_i-\bar x_w)(y_i-\bar y_w)}{\sum_i w_i(x_i-\bar x_w)^2},\qquad \hat\beta_0=\bar y_w-\hat\beta_1\,\bar x_w. \]
Dengan \(w_i=1/|x_i|\) (proporsional dari varian), keduanya memberi solusi yang sama.
4) Cara transformasi setara (homoskedastisasi)
Karena \(w_i=1/|x_i|\), ambil \(z_i=\sqrt{w_i}=1/\sqrt{|x_i|}\). Bentuk variabel bertanda bintang: \[ y_i^*=z_i y_i,\quad x_i^*=z_i x_i,\quad \text{dan kolom intercept}~1^*=z_i. \] Lalu jalankan OLS pada model \[ y_i^*=\beta_0\,1_i^*+\beta_1\,x_i^*+u_i^*. \] Hasil \((\hat\beta_0,\hat\beta_1)\) sama dengan WLS di atas.
5) Catatan penting
- Jika ada \(x_i=0\), maka \(\mathrm{Var}(\varepsilon_i|x_i)=0\) dan \(w_i=1/|x_i|\) tak terdefinisi. Praktiknya, berikan toleransi kecil \(w_i=1/(|x_i|+\epsilon)\) atau tangani titik tersebut secara khusus.
- Semua turunan di atas mengasumsikan model mean benar dan eksogenitas \(\mathbb{E}[\varepsilon_i|x_i]=0\).
Esai: Bandingkan Ridge dan Lasso dari sudut pandang bias–variance dan seleksi variabel.
Jawaban :
Gambaran besar
Ridge dan Lasso sama-sama mengecilkan koefisien untuk mengurangi overfitting. Pengecilan menambah bias tetapi menurunkan varians, sehingga galat uji bisa turun. Bedanya ada pada koefisiennya:
- Ridge (L2): mengecilkan koefisien mendekati nol.
- Lasso (L1): bisa membuat sebagian koefisien persis nol.
Bias dan varians
- Ridge
- Bias: naik mulus seiring λ membesar.
- Varians: turun stabil. Sering lebih kuat menurunkan varians dibanding Lasso untuk tingkat penyusutan yang sama.
- Efek keseluruhan: bagus saat banyak efek kecil/menengah atau ada multikolinearitas tinggi.
- Lasso
- Bias: naik, kadang tajam saat melewati titik di mana koefisien jadi nol.
- Varians: turun, tetapi langkah seleksi dapat membuat jalurnya lebih berubah-ubah daripada Ridge.
- Efek keseluruhan: unggul saat sinyal benar-benar jarang (hanya sedikit prediktor yang relevan).
Seleksi variabel
- Ridge: tidak melakukan seleksi dalam arti ketat. Koefisien jarang persis nol, jadi semua prediktor tetap ada (hanya bobotnya mengecil). Dengan prediktor yang berkorelasi, Ridge cenderung membagi bobot di antaranya (grouping effect).
- Lasso: melakukan seleksi. Sebagian koefisien jadi nol sehingga model lebih sederhana dan mudah ditafsirkan. Pada prediktor yang sangat berkorelasi, Lasso sering memilih satu dan mengabaikan lainnya, yang bisa tidak stabil antar-sampel.
Intuisi geometri: mengapa Lasso bisa nol
- Ridge memakai kendala bola L2 (permukaan halus). Solusi biasanya jatuh di permukaan yang bukan sumbu, jadi koefisien menyusut tetapi tidak nol.
- Lasso memakai belah ketupat L1 (sudut tajam di sumbu). Solusi sering mendarat tepat di sudut sumbu, membuat sebagian koefisien nol.
Panduan praktis
- Perkirakan banyak sinyal kecil atau multikolinearitas berat: pilih Ridge untuk prediksi yang stabil.
- Menginginkan seleksi fitur dan menduga model sparse: pilih Lasso.
- Jika prediktor sangat berkorelasi dan Anda tetap ingin seleksi, pertimbangkan Elastic Net (L1 + L2) agar mendapat sparsity ala Lasso dan stabilitas ala Ridge.
Perbandingan cepat
| Aspek | Ridge (L2) | Lasso (L1) |
|---|---|---|
| Bias (naik saat λ naik) | Naik mulus | Bisa melonjak di titik seleksi |
| Varians | Turun kuat, sangat stabil | Turun, tetapi seleksi bisa menambah ketidakstabilan |
| Seleksi variabel | Tidak (jarang nol) | Ya (nol persis) |
| Fitur berkorelasi | Berbagi bobot (grouping) | Cenderung pilih satu, buang lainnya |
| Cocok saat… | Banyak efek kecil; multikolinearitas | Model sparse; butuh interpretasi |
Pilih λ dengan cross-validation. Jika stabilitas seleksi penting, kombinasikan Lasso dengan stability selection atau refit OLS/Ridge pada fitur terpilih untuk mengurangi bias.
Hitungan: Diberikan matriks 𝑋 dan vektor 𝑦, hitunĝ 𝛽 OLS serta 𝐻; identifikasi leverage mak- simum.
Jawaban :
Estimator OLS
Asumsikan \(X\) berukuran \(n \times p\) dengan rank penuh.
Estimator OLS:
\[ \hat{\beta} = (X'X)^{-1}X'y \]
Matriks Hat
Matriks hat (proyeksi ke ruang kolom \(X\)):
\[ H = X (X'X)^{-1} X' \]
Leverage
Leverage observasi ke-\(i\) adalah diagonal \(h_{ii}\) dari \(H\).
Leverage maksimum:
\[ \max_i h_{ii}, \quad \text{dan indeksnya } \arg\max_i h_{ii} \]
Implementasi R
## INPUT: X (n x p), y (n x 1)
## Jika butuh intercept, tambahkan kolom 1 di X: X <- cbind(1, X)
ols_with_hat <- function(X, y) {
stopifnot(nrow(X) == length(y))
XtX <- crossprod(X) # X'X
Xty <- crossprod(X, y) # X'y
beta <- solve(XtX, Xty) # (X'X)^(-1) X'y -- stabil: solve(XtX, Xty)
# Matriks hat H = X (X'X)^(-1) X'
# H full (n x n) bisa besar; untuk leverage cukup diagonalnya:
XtX_inv <- solve(XtX)
# leverage h = diag( X (X'X)^(-1) X' ) = rowSums( (X %*% XtX_inv) * X )
h <- rowSums( (X %*% XtX_inv) * X )
# Jika Anda tetap ingin H lengkap:
H <- X %*% XtX_inv %*% t(X)
# Leverage maksimum
h_max_val <- max(h)
h_max_idx <- which.max(h)
list(
beta = as.vector(beta),
H = H,
leverage = h,
h_max_value = h_max_val,
h_max_index = h_max_idx
)
}
## Contoh kecil (hapus atau ganti dengan data Anda):
# X <- cbind(1, c(1,2,3,4)) # intercept + satu prediktor
# y <- c(2, 3, 5, 4)
# out <- ols_with_hat(X, y)
# out$beta
# out$h_max_value; out$h_max_index2 Regresi Nonlinier: Model, Estimasi, dan Aplikasi
2.1 Tujuan Pembelajaran Analisis Regresi Nonlinear
- Memperkenalkan berbagai permasalahan yang relevan dengan pemodelan dan pemfittingan fungsi regresi nonlinear.
- Menyajikan representasi grafis untuk menilai kualitas dari interval kepercayaan aproksimasi.
- 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.96Plot
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.6569Iterasi 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-9gn_res$theta## beta0 beta1
## 1.9753264 0.3148085Kurva 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 TRUElm_res$theta## [1] 1.9753264 0.3148085Kurva 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-07coef(fit_nls)## b0 b1
## 1.9753263 0.3148085Inferensi (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.01969227Interval 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.3534053SOAL
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.
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.
Diagnostik Residu Untuk hasil terbaik (GN/LM), plot residu vs fitted dan QQ-plot residu. Apakah asumsi normalitas dan homoskedastisitas terlihat wajar?
Tambahan Terapkan GN/LM untuk model Michaelis–Menten dan logistik. Laporkan \(\hat{\theta}\), SE, CI, dan interpretasi substantif parameternya.
Jawaban
- 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.3148085Kesimpulan
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.
- 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.2937094Kesimpulan:
Dengan error multiplikatif, model linear log cocok → estimasi lebih sederhana (OLS), hasilnya dekat dengan nilai parameter sebenarnya.
- 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.
- 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-06confint(fit_mm)## 2.5% 97.5%
## Vmax 1.952336 2.04651
## Km 1.146484 1.34395Interpretasi
- \(\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:
- Apakah parameter signifikan ( 0)?
- Bagaimana ketidakpastian estimasi (interval kepercayaan)?
- 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- 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.313Interpretasi: jika p-value < 0.05 → parameter signifikan berbeda dari 0.
- Interval Kepercayaan 95%
confint(fit)## 2.5% 97.5%
## b0 1.7917267 2.2341574
## b1 0.2879737 0.3126612Interval kepercayaan memberi rentang nilai parameter yang konsisten dengan data.
- 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.805362.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.
Jawab :
Kita memodelkan tinggi pohon (height)
sebagai fungsi umur (age) pada dataset
Loblolly. Pola pertumbuhan pohon umumnya
nonlinear: cepat di awal, melambat, lalu mendekati
batas (kapasitas). Bentuk ini cocok dimodelkan dengan kurva
logistik.
Rumus Model
Kita gunakan logistik 3-parameter:
\[ h(a) \;=\; \frac{\text{Asym}}{1 + \exp\!\left(\frac{\text{xmid} - a}{\text{scal}}\right)} \;+\; \varepsilon, \quad \varepsilon \sim \mathcal{N}(0,\sigma^2), \] dengan:
Asym: tinggi maksimum (asimtot),
xmid: umur saat tinggi ≈ Asym/2 (titik belok),
scal: skala yang mengatur curamnya kurva (sekitar xmid, kemiringan ≈ Asym/(4·scal)).
Kode:
Fit nls(), ringkasan parameter, goodness-of-fit, dan
plot
suppressPackageStartupMessages({
library(dplyr); library(ggplot2); library(broom); library(knitr)
})
# 1) Data ------------------------------------------------------
data(Loblolly)
# Pilih satu pohon (Seed) agar kurva tidak tercampur antar individu
# Cek daftar Seed yang tersedia:
unique(Loblolly$Seed)## [1] 301 303 305 307 309 311 315 319 321 323 325 327 329 331
## 14 Levels: 329 < 327 < 325 < 307 < 331 < 311 < 315 < 321 < 319 < ... < 305seed_pilih <- "329" # ganti sesuai kebutuhan, mis. "301", "307", "319", dst.
dat <- Loblolly %>%
filter(Seed == seed_pilih) %>%
arrange(age)
if (nrow(dat) < 6) stop("Observasi terlalu sedikit. Pilih Seed lain.")
# 2) Fit NLS (logistik self-start) -----------------------------
# SSlogis memberi start value otomatis yang stabil
fit <- nls(height ~ SSlogis(age, Asym, xmid, scal), data = dat)
# Ringkasan koefisien dan SE (Wald)
sm <- summary(fit)
coef_hat <- coef(fit)
se_hat <- sm$parameters[, "Std. Error"]
# CI 95% (Wald ~ hampiran normal)
ci95 <- cbind(
lower = coef_hat - 1.96 * se_hat,
est = coef_hat,
upper = coef_hat + 1.96 * se_hat
)
# 3) Goodness-of-fit -------------------------------------------
yhat <- fitted(fit)
res <- resid(fit)
SSE <- sum(res^2)
SST <- sum( (dat$height - mean(dat$height))^2 )
R2 <- 1 - SSE/SST
AIC <- nrow(dat) * log(SSE / nrow(dat)) + 2 * length(coef_hat)
cat("== RINGKASAN PARAMETER (Seed =", seed_pilih, ") ==\n")## == RINGKASAN PARAMETER (Seed = 329 ) ==print(sm$coefficients)## Estimate Std. Error t value Pr(>|t|)
## Asym 57.337166 3.9261432 14.603942 0.0006962704
## xmid 11.738683 1.0244671 11.458331 0.0014266745
## scal 4.292549 0.7299197 5.880851 0.0098105679cat(sprintf("\nSSE = %.3f | R^2 = %.4f | AIC = %.2f\n", SSE, R2, AIC))##
## SSE = 24.517 | R^2 = 0.9889 | AIC = 14.45cat("\nCI 95% (Wald):\n"); print(ci95, digits = 4)##
## CI 95% (Wald):## lower est upper
## Asym 49.642 57.337 65.032
## xmid 9.731 11.739 13.747
## scal 2.862 4.293 5.723# 4) Visualisasi: data vs kurva logistik -----------------------
age_grid <- seq(min(dat$age), max(dat$age), length.out = 200)
pred_grid <- data.frame(
age = age_grid,
height_hat = with(as.list(coef_hat),
Asym / (1 + exp((xmid - age_grid)/scal)))
)
ggplot(dat, aes(age, height)) +
geom_point(size = 2) +
geom_line(data = pred_grid, aes(age, height_hat), linewidth = 1) +
labs(
title = paste("Loblolly (Seed", seed_pilih, "): Data vs Kurva Logistik (nls)"),
x = "Umur (tahun)", y = "Tinggi (unit bawaan data)"
) +
theme_minimal()
# 5) Diagnostik residu -----------------------------------------
par(mfrow = c(1, 2))
plot(yhat, res,
xlab = "Fitted", ylab = "Residuals",
main = "Residu vs Fitted (Logistik)")
abline(h = 0, lty = 2, col = "gray40")
qqnorm(res, main = "QQ-plot Residu (Logistik)")
qqline(res, col = "gray40")
par(mfrow = c(1, 1))Studi Kasus: Loblolly — Fit Model Nonlinear (Logistik) dengan
nls()
Interpretasi biologis
- Asym: perkiraan tinggi maksimum pohon (satuannya
sama dengan
height). Makin besar nilainya, makin tinggi potensi akhir pohon. - xmid: umur ketika tinggi pohon kira-kira mencapai setengah Asym. Titik ini menandai peralihan dari fase tumbuh cepat ke melambat.
- scal: seberapa tajam perubahan di sekitar xmid. Nilai scal yang lebih kecil berarti transisi dari lambat → cepat → lambat terjadi dalam rentang usia yang lebih sempit.
Goodness-of-fit
Jika SSE kecil, \(R^2\) tinggi, dan AIC lebih rendah daripada model linear, kurva logistik biasanya lebih pas untuk pola pertumbuhan yang cenderung menjenuh.
Cek residu: sebaran yang acak di sekitar nol (tidak membentuk pola U atau ∩) dan QQ-plot yang mengikuti garis lurus menunjukkan asumsi galat normal yang cukup wajar.
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.
Jawab
Loblolly: Linear vs Nonlinear (Logistik) dengan
nls()
Tujuan
Fit model linear:
height ~ ageFit model nonlinear logistik:
height ~ SSlogis(age, Asym, xmid, scal)Bandingkan SSE, R², AIC dan visualisasi
suppressPackageStartupMessages({
library(dplyr); library(ggplot2); library(broom); library(knitr)
})
# 1) Ambil satu pohon (Seed) agar kurva tidak tercampur antar individu
data(Loblolly)
seed_pilih <- "329" # ganti bila perlu: "301", "307", "319", dst.
dat <- Loblolly %>% filter(Seed == seed_pilih) %>% arrange(age)
if (nrow(dat) < 6) stop("Observasi terlalu sedikit. Pilih Seed lain.")
# 2) Fit model
fit_lin <- lm(height ~ age, data = dat)
fit_nls <- nls(height ~ SSlogis(age, Asym, xmid, scal), data = dat)
# 3) Metrik: SSE, R^2 (pseudo untuk nls), AIC
SST <- sum( (dat$height - mean(dat$height))^2 )
# Linear
yhat_lin <- fitted(fit_lin)
res_lin <- resid(fit_lin)
SSE_lin <- sum(res_lin^2)
R2_lin <- 1 - SSE_lin / SST # sama dengan summary(lm)$r.squared
AIC_lin <- nrow(dat) * log(SSE_lin / nrow(dat)) + 2 * length(coef(fit_lin))
# Nonlinear (logistik)
yhat_nls <- fitted(fit_nls)
res_nls <- resid(fit_nls)
SSE_nls <- sum(res_nls^2)
R2_nls <- 1 - SSE_nls / SST
AIC_nls <- nrow(dat) * log(SSE_nls / nrow(dat)) + 2 * length(coef(fit_nls))
# 4) Ringkasan parameter + metrik
cat("== Koefisien Linear (lm) ==\n"); print(summary(fit_lin)$coefficients)## == Koefisien Linear (lm) ==## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.260856 2.2718677 -0.5549865 6.084748e-01
## age 2.428784 0.1495859 16.2367179 8.418892e-05cat("\n== Koefisien Logistik (nls) ==\n"); print(summary(fit_nls)$coefficients)##
## == Koefisien Logistik (nls) ==## Estimate Std. Error t value Pr(>|t|)
## Asym 57.337166 3.9261432 14.603942 0.0006962704
## xmid 11.738683 1.0244671 11.458331 0.0014266745
## scal 4.292549 0.7299197 5.880851 0.0098105679cmp <- tibble(
Model = c("Linear (lm)", "Logistik (nls)"),
SSE = c(SSE_lin, SSE_nls),
R2 = c(R2_lin, R2_nls),
AIC = c(AIC_lin, AIC_nls)
)
cat("\n== Perbandingan Metrik (Seed =", seed_pilih, ") ==\n")##
## == Perbandingan Metrik (Seed = 329 ) ==print(knitr::kable(cmp, digits = 4))##
##
## |Model | SSE| R2| AIC|
## |:--------------|-------:|------:|-------:|
## |Linear (lm) | 33.1164| 0.9851| 14.2496|
## |Logistik (nls) | 24.5171| 0.9889| 14.4457|# 5) Visualisasi: data + dua kurva fit
# Grid umur untuk kurva halus
age_grid <- seq(min(dat$age), max(dat$age), length.out = 200)
# Prediksi dari kedua model di grid
pred_nls <- predict(fit_nls, newdata = data.frame(age = age_grid))
pred_lin <- predict(fit_lin, newdata = data.frame(age = age_grid))
# Gabungkan ke data frame
df_grid <- data.frame(
age = age_grid,
y_nls = pred_nls,
y_lin = pred_lin
)
# Plot: titik data + dua kurva
library(ggplot2)
ggplot(dat, aes(age, height)) +
geom_point(size = 2) +
geom_line(data = df_grid, aes(x = age, y = y_nls),
linewidth = 1, inherit.aes = FALSE) +
geom_line(data = df_grid, aes(x = age, y = y_lin),
linewidth = 1, linetype = 2, color = "gray40",
inherit.aes = FALSE) +
labs(
title = paste("Loblolly (Seed", seed_pilih, "): Linear vs Logistik"),
subtitle = "Garis penuh: logistik (nls) | Garis putus: linear (lm)",
x = "Umur (tahun)", y = "Tinggi"
) +
theme_minimal()
# 6) Diagnostik residu: kedua model
par(mfrow = c(2, 2))
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")
plot(yhat_nls, res_nls,
xlab = "Fitted (Logistik)", ylab = "Residuals",
main = "Residu vs Fitted — Logistik")
abline(h = 0, lty = 2, col = "gray40")
qqnorm(res_nls, main = "QQ-plot — Logistik"); qqline(res_nls, col = "gray40")
par(mfrow = c(1, 1))Kesimpulan
- SSE & AIC: biasanya lebih baik
(SSE lebih kecil, AIC lebih rendah) pada logistik,
karena pertumbuhan tinggi pohon cenderung melengkung dan
menjenuh.
- R²: umumnya lebih tinggi pada
logistik untuk pola pertumbuhan seperti ini.
- Visual: garis linear sulit mengikuti pelandaian
pada umur besar, sedangkan logistik menyesuaikan bentuk saturasi.
- Residu: pada model linear sering muncul pola (U/∩), tanda spesifikasi kurang pas. Residu logistik cenderung lebih acak di sekitar nol.
2.5.3 Simulasi: Buat simulasi data eksponensial
- Gunakan nls() untuk estimasi parameter.
- Laporkan hasil dan interpretasi.
Jawab :
Simulasi data eksponensial + estimasi nls()
Kita mensimulasikan data dari model eksponensial dengan galat aditif,
lalu menaksir parameternya menggunakan Nonlinear Least
Squares (nls). Setelah itu kita laporkan
koefisien, SE/CI, serta metrik kecocokan (SSE,
R², AIC) dan melihat plot data vs kurva hasil estimasi serta
diagnostik residu.
Model
\[
y(t) \;=\; b_0\,e^{\,b_1 t} \;+\; \varepsilon,\quad \varepsilon \sim
\mathcal{N}(0,\sigma^2).
\] - \(b_0\): level dasar saat
\(t=0\)
- \(b_1\): laju pertumbuhan
- \(\sigma\): simpangan baku galat
set.seed(202)
# 1) Simulasi data
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)
# 2) Estimasi NLS
fit <- nls(
y ~ b0 * exp(b1 * t),
start = list(b0 = 1, b1 = 0.1),
trace = FALSE,
control = nls.control(maxiter = 200, warnOnly = TRUE)
)
# 3) Ringkasan, SE, CI 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
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 hasil
cat("== KOEFISIEN (taksiran) & SE ==\n")## == KOEFISIEN (taksiran) & 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-39cat("\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.4126cat(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.88Visualisasi kurva (data vs hasil nls)
# Kurva halus dari hasil nls dan kurva 'true' sebagai pembanding
tt <- seq(min(t), max(t), length.out = 200)
yy_hat <- coef_hat[["b0"]] * exp(coef_hat[["b1"]] * tt)
plot(t, y, pch = 19, xlab = "t", ylab = "y",
main = "Data simulasi vs kurva eksponensial (NLS)")
lines(tt, yy_hat, lwd = 2) # hasil nls
lines(tt, b0_true * exp(b1_true * tt), lty = 2) # kurva true (opsional)
legend("topleft",
legend = c("Data", "Fitted (nls)", "True"),
pch = c(19, NA, NA), lty = c(NA, 1, 2), lwd = c(NA, 2, 1),
bty = "n")
Diagnostik residu
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))Interpretasi
- Taksiran \(\hat{b}_0\) dan \(\hat{b}_1\) diharapkan mendekati nilai simulasi \((1.5,\,0.4)\).
- SSE kecil dan \(R^2\) tinggi menandakan model menangkap pola eksponensial dengan baik; AIC dipakai untuk membandingkan dengan model alternatif.
- Residu yang menyebar acak di sekitar nol dan QQ-plot yang mendekati garis lurus mendukung asumsi galat normal yang wajar.
- Jika \(\hat{b}_1\) signifikan (CI tidak meliputi nol), ada bukti kuat pertumbuhan eksponensial pada data simulasi ini.
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)
Pilih tebakan awal \(\theta^{(0)}\) (mis. via regresi \(\ln y\) pada \(x\) untuk inisialisasi).
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.
- 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
- Pilih tebakan awal \(\theta^{(0)}\), pilih \(\lambda^{(0)} > 0\).
- 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.
- Hentikan bila kriteria konvergensi tercapai.
- 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.
3 Robust Regression
Bab ini menyajikan materi komprehensif tentang outlier dalam analisis regresi, mencakup:
- (i) definisi dan klasifikasi,
- (ii) teknik identifikasi grafis dan statistik,
- (iii) strategi penanganan, serta
- (iv) ide dasar dan praktik robust regression berikut contoh kode R yang siap dijalankan.
3.1 Konsep Dasar Outlier
Outlier adalah observasi yang menyimpang jauh dari pola umum data. Dalam regresi, pembedaan penting:
- Outlier di ruang respons (Y-space): residual sangat besar.
- Outlier di ruang prediktor (X-space): nilai prediktor jauh dari pusat sebaran.
- Leverage point: observasi dengan nilai X ekstrem (potensi memengaruhi slope).
- Influential point: observasi yang jika dihapus, koefisien berubah nyata (diukur mis. Cook’s distance).
Secara umum, regresi OLS meminimalkan jumlah kuadrat residual
\[\hat{\beta} = \arg\min_{\beta} \sum_i (y_i - x_i^\top \beta)^2\]
sehingga residual besar akan dikuadratkan dan diberi bobot tidak proporsional — ini yang membuat OLS rentan terhadap outlier.
3.2 Identifikasi Outlier
3.2.1 Pemeriksaan Grafis
Scatterplot (untuk regresi sederhana) atau added-variable plot (untuk multivariat).
Residual plot untuk melihat pola non-acak dan titik residual ekstrem.
QQ-plot residual untuk melihat deviasi dari normalitas.
set.seed(123)
x <- 1:60
y <- 2 + 0.6*x + rnorm(length(x), 0, 3)
y[c(7, 35, 50)] <- y[c(7, 35, 50)] + c(18,-22, 25) # tiga outlier
ols <- lm(y~ x)
par(mfrow=c(1,3))
plot(x, y, pch=19, main="Scatterplot y vs x"); abline(ols, lwd=2)
plot(ols, which=1) # Residuals vs Fitted
qqnorm(rstandard(ols), main="QQ-plot Standardized Residual"); qqline(rstandard(ols))
par(mfrow=c(1,1))3.2.2 Diagnostik Statistik
Notasi: \(\hat{e}_i\) residual OLS, \(\hat{\sigma}^2\) ragam residual, \(h_{ii}\) nilai diagonal hat matrix \(H = X(X^\top X)^{-1}X^\top\).
- Standardized/Studentized residual
\[r_i = \dfrac{\hat{e}_i}{\hat{\sigma}\,\sqrt{1 - h_{ii}}}, \quad t_i = \dfrac{\hat{e}_i}{\hat{\sigma}_{(i)}\,\sqrt{1 - h_{ii}}}.\]
Kriteria praktis: \(|r_i| > 2\) atau \(|t_i| > 3\) untuk kandidat outlier.
- Leverage
\[h_{ii} = x_i^\top (X^\top X)^{-1} x_i$, ambang populer: $h_{ii} > 2p\].
dengan 𝑝 jumlah parameter (termasuk intercept), 𝑛 ukuran sampel.
- Cook’s Distance
\[D_i = \dfrac{r_i^2}{p}\cdot \dfrac{h_{ii}}{1 - h_{ii}}.\] Aturan praktis: 𝐷𝑖 > 1 perlu ditinjau.
- DFFITS dan DFBETAS
\[\mathrm{DFFITS}_i \;=\; \dfrac{\hat{y}_i - \hat{y}_i^{(i)}}{\hat{\sigma}_{(i)}\sqrt{h_{ii}}}, \qquad \mathrm{DFBETAS}_{ij} \;=\; \dfrac{\hat{\beta}_j - \hat{\beta}_{j}^{(i)}}{SE_{(i)}(\hat{\beta}_j)}.\]
library(broom)
library(dplyr)
aug <- augment(ols) %>%
mutate(
h_ii = hatvalues(ols),
rstd = rstandard(ols),
cooks = cooks.distance(ols)
)
head(aug, 10)## # A tibble: 10 × 11
## y x .fitted .resid .hat .sigma .cooksd .std.resid h_ii rstd
## <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.919 1 3.05 -2.13 0.0650 5.59 5.49e-3 -0.397 0.0650 -0.397
## 2 2.51 2 3.65 -1.14 0.0618 5.60 1.49e-3 -0.213 0.0618 -0.213
## 3 8.48 3 4.26 4.22 0.0587 5.57 1.91e-2 0.784 0.0587 0.784
## 4 4.61 4 4.86 -0.249 0.0557 5.60 6.27e-5 -0.0461 0.0557 -0.0461
## 5 5.39 5 5.46 -0.0756 0.0528 5.60 5.46e-6 -0.0140 0.0528 -0.0140
## 6 10.7 6 6.07 4.68 0.0500 5.56 1.97e-2 0.865 0.0500 0.865
## 7 25.6 7 6.67 18.9 0.0474 4.97 3.03e-1 3.49 0.0474 3.49
## 8 3.00 8 7.27 -4.27 0.0448 5.57 1.45e-2 -0.787 0.0448 -0.787
## 9 5.34 9 7.88 -2.54 0.0424 5.59 4.83e-3 -0.467 0.0424 -0.467
## 10 6.66 10 8.48 -1.82 0.0400 5.59 2.33e-3 -0.334 0.0400 -0.334
## # ℹ 1 more variable: cooks <dbl>par(mfrow=c(1,2))
plot(aug$rstd, pch=19, main="Standardized Residual", ylab="r_i"); abline(h=c(-2,2), lty=2)
plot(aug$cooks, pch=19, main="Cook's Distance", ylab="D_i"); abline(h=1, lty=2)
par(mfrow=c(1,1))Catatan: Gunakan beberapa metrik sekaligus untuk keputusan yang lebih andal (mis. kandidat pengaruh jika rstd besar dan \(h_{ii}\) tinggi dan \(D_i\) besar).
3.2.3 Strategi Penanganan Outlier
3.2.3.1 Verifikasi Kualitas Data
Koreksi data entry error atau coding error.
Konfirmasi konteks substantif (outlier bermakna versus noise).
3.2.3.2 Transformasi dan Rekayasa Fitur
Transformasi log, sqrt, Box–Cox/Yeo–Johnson untuk menstabilkan ragam dan mengurangi pengaruh ekstrem.
Binning atau winsorizing (mengganti nilai di luar persentil 1–99).
y_log <- log(pmax(y - min(y) + 1, 1)) # contoh log; pastikan argumen positif
par(mfrow = c(1, 2))
plot(x, y, pch = 19, main = "Asli")
plot(x, y_log, pch = 19, main = "Transformasi log(y)")
par(mfrow=c(1,1))3.2.3.3 Model & Estimasi yang Tahan
- Weighted least squares (WLS): beri bobot kecil pada observasi dengan varians besar.
- Model robust (M-estimator, LAD, RLM-IRLS).
- Quantile regression: memusat pada median/kuantil lain, lebih stabil terhadap outlier pada mean.
- Model alternatif: GLM dengan link dan family sesuai sifat error (bila masalahnya bukan sekadar outlier).
3.2.3.4 Dokumentasi Keputusan
- Laporkan kriteria, titik yang terpengaruh, perubahan koefisien, serta implikasi substantif.
- Jangan sekadar menghapus data tanpa alasan yang dapat diaudit.
3.2.4 Ide & Metode Robust Regression
3.2.4.1 Mengapa Robust?
Fungsi kuadrat memberi pengaruh tak terbatas (unbounded influence) pada residual besar. Robust regression mengganti fungsi kerugian kuadrat dengan kerugian less than quadratic atau bounded, sehingga pengaruh outlier dibatasi.
3.2.4.2 Keluarga M‑Estimator
M‑estimasi meminimalkan
\[ \hat{\beta} \;=\; \arg\min_{\beta}\; \sum_{i=1}^{n} \rho\!\left(\frac{y_i - x_i^{\top}\beta}{\sigma}\right). \]
dengan fungsi loss 𝜌(⋅) tertentu dan fungsi influence 𝜓 = 𝜌′. Dua fungsi populer:
- Huber loss dengan ambang 𝑐:
$$ - Tukey’s bisquare (redescending):
\[ \rho(e) = \begin{cases} \dfrac{c^{2}}{6}\left[1 - \left(1 - \left(\dfrac{e}{c}\right)^{2}\right)^{3}\right], & |e|\le c,\\[8pt] \dfrac{c^{2}}{6}, & |e|>c, \end{cases} \qquad \psi(e) = \begin{cases} e\left(1 - \left(\dfrac{e}{c}\right)^{2}\right)^{2}, & |e|\le c,\\[8pt] 0, & |e|>c. \end{cases} \] Estimasi praktis sering menggunakan IRLS (Iteratively Reweighted Least Squares) dengan bobot
\[ w_i = \frac{\psi(e_i)}{e_i}. \] sehingga solusi pada iterasi \(t\) adalah pemecahan WLS dengan bobot \(w_i^{(t)}\).
3.2.4.3 Menentukan nilai \(c\)**
Nilai \(c\) adalah tuning constant yang menentukan kapan sebuah residual dianggap outlier.
Jika \(|e| \le c\), maka residual dianggap normal → diperlakukan seperti OLS.
Jika \(|e| > c\), maka residual dianggap outlier → diberi pengaruh terbatas.
Cara memilih \(c\)
- Huber Loss
Nilai \(c \approx 1.345\,\hat{\sigma}\)
sering digunakan karena menghasilkan estimator dengan efisiensi ~95%
dibanding OLS saat error normal.
- Tukey’s Bisquare
Nilai \(c \approx 4.685\,\hat{\sigma}\)
banyak dipakai agar estimator efisien mendekati OLS di distribusi
normal.
Artinya: \(c\) dikalibrasi berdasarkan skala error (\(\hat{\sigma}\)), supaya fungsi robust tetap efisien kalau data normal, tapi membatasi outlier.
Fungsi influence
Fungsi influence (\(\psi(e)\)) adalah
turunan dari fungsi loss \(\rho(e)\):
\[
\psi(e) = \frac{d}{de}\,\rho(e).
\] Interpretasi
- \(\psi(e)\) menunjukkan seberapa besar pengaruh sebuah residual terhadap estimasi koefisien.
- Pada OLS: \(\psi(e)=e\). Residual makin besar → pengaruh makin besar (tanpa batas).
- Pada Huber dan Tukey: \(\psi(e)\) dibatasi. Residual besar tidak bisa “mengendalikan” hasil.
Contoh - Huber: \[ \psi(e) = \begin{cases} e, & |e|\le c,\\[6pt] c\,\operatorname{sign}(e), & |e|>c \end{cases} \]
Efek outlier dibatasi pada \(\pm c\).
Tukey: \[ \psi(e) = \begin{cases} e\left(1 - \left(\dfrac{e}{c}\right)^2\right)^2, & |e|\le c,\\[8pt] 0, & |e|>c \end{cases} \]
outlier ekstrem bahkan diberi bobot nol.
Bobot dalam Iteratively Reweighted Least Squares (IRLS)
Dalam IRLS, setiap iterasi menghitung bobot \(w_i\) berdasarkan residual saat itu:
\[ w_i \;=\; \dfrac{\psi(e_i)}{e_i} \] Sehingga regresi robust menjadi Weighted Least Squares (WLS) dengan bobot berbeda di tiap observasi.
Contoh Perhitungan Bobot
Misal residual dari model OLS awal adalah: \(e=\{-0.5,\; 1.2,\; 5.0\}\), \(c=1.345\).
Huber loss:
- Untuk \(e = -0.5\): \(|e| < c \Rightarrow \psi(e) = -0.5 \Rightarrow w = 1.\)
- Untuk \(e = 1.2\): \(|e| < c \Rightarrow \psi(e) = 1.2 \Rightarrow w = 1.\)
- Untuk \(e = 5.0\): \(|e| > c \Rightarrow \psi(e) = 1.345 \Rightarrow w = 1.345/5 = 0.269.\)
Observasi ke-3 (outlier) hanya dapat bobot \(0.269\), jauh lebih kecil dibanding observasi lain.
Iterasi IRLS
- Fit OLS \(\rightarrow\) dapat residual awal.
- Hitung \(w_i\) untuk semua observasi.
- Fit ulang dengan WLS (berat \(w_i\)).
- Update residual \(\rightarrow\) hitung bobot baru \(\rightarrow\) ulangi.
- Berhenti saat perubahan koefisien kecil (konvergen).
Contoh R (Huber dengan IRLS)
library(MASS)
# Data dengan outlier
set.seed(123)
x <- 1:20
y <- 2*x + rnorm(20, 0, 2)
y[c(5, 15)] <- y[c(5, 15)] + 20 # tambahkan outlier
# Robust regression (Huber M-estimator)
fit_rob <- rlm(y~ x, psi = psi.huber)
# Lihat bobot yang diberikan
weights(fit_rob)## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1Output weights(fit_rob) akan menunjukkan:
- Hampir semua observasi bobot ~1
- Observasi outlier → bobot jauh lebih kecil (<1)
Estimasi Skala \(\hat{\sigma}\)
\(\hat{\sigma}\) pada robust regression
tidak sama dengan standar deviasi OLS.
Karena OLS sensitif terhadap outlier, dipakai estimasi skala robust.
Median Absolute Deviation (MAD)
\[
\hat{\sigma} = 1.4826 \cdot \operatorname{median}\bigl(\,|e_i -
\operatorname{median}(e)|\,\bigr)
\]
Konstanta 1.4826 membuatnya konsisten dengan simpangan baku normal.
Update Iteratif dalam IRLS
Setelah setiap iterasi, skala \(\sigma\) diupdate dengan formula berbasis
\(\psi(e)\).
Implementasi di R (fungsi rlm() dari paket
MASS):
- Menggunakan MAD residual OLS sebagai starting value.
- Update \(\hat{\sigma}\) dengan metode robust (Huber’s Proposal 2).
fit_rob$s # nilai sigma robust## [1] 2.487055summary(fit_rob)##
## Call: rlm(formula = y ~ x, psi = psi.huber)
## Residuals:
## Min 1Q Median 3Q Max
## -4.6241 -1.6307 -0.0158 1.7847 19.3656
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 0.9707 1.2954 0.7493
## x 1.9845 0.1081 18.3511
##
## Residual standard error: 2.487 on 18 degrees of freedomCatatan:
- Nilai \(c\): ambang residual dianggap outlier (biasanya proporsional dengan \(\hat{\sigma}\)).
- Influence function: ukuran pengaruh residual, dibatasi dalam robust regression.
- Bobot IRLS: membatasi kontribusi outlier, membuat model stabil.
- \(\hat{\sigma}\): estimasi skala robust (MAD atau update IRLS), bukan simpangan baku OLS.
3.2.4.4 LAD (Least Absolute Deviation)
Meminimalkan: \[ \min_{\beta}\ \sum_{i=1}^{n} \big|\, y_i - x_i^{\top}\beta \,\big| \]
- Lebih tahan terhadap outlier pada \(Y\), namun kurang efisien jika error
benar-benar normal.
- Dapat diselesaikan dengan pemrograman linear atau algoritme khusus.
3.2.4.5 Quantile Regression
Menggeneralisasi LAD ke kuantil \(\tau \in (0, 1)\) dengan meminimalkan check loss \[ \sum_{i=1}^{n} \rho_\tau(\varepsilon_i), \qquad \rho_\tau(u) = u \cdot \big(\tau - \mathbf{1}\{u < 0\}\big). \] Memungkinkan analisis efek heterogen di berbagai titik distribusi respons (median, P90, dst.).
3.2.5 Praktik di R: OLS vs Robust
library(MASS) # rlm (Huber & bisquare)
library(quantreg) # quantile regression
library(dplyr); library(broom)3.2.5.1 Simulasi & OLS
set.seed(42)
n <- 120
x1 <- runif(n, 0, 10)
x2 <- rnorm(n)
y <- 3 + 2*x1- 1.5*x2 + rnorm(n, 0, 2)
# Tambah outlier respons & leverage
y[c(10, 55)] <- y[c(10, 55)] + c(25,-30)
x1[c(80, 95)] <- x1[c(80, 95)] + c(15, 18)
dat <- data.frame(y, x1, x2)
fit_ols <- lm(y~ x1 + x2, dat)
glance(fit_ols)## # 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.476 0.467 5.53 53.1 3.96e-17 2 -374. 756. 767.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>Diagnosis pengaruh
h <- hatvalues(fit_ols)
cooks <- cooks.distance(fit_ols)
rstd <- rstandard(fit_ols)
cand_idx <- which( (abs(rstd) > 2) | (h > 2*length(coef(fit_ols))/n) | (cooks > 1) )
length(cand_idx); cand_idx## [1] 7## 10 44 55 58 68 80 95
## 10 44 55 58 68 80 953.2.5.2 Robust (Huber & Bisquare)
fit_hub <- rlm(y ~ x1 + x2, data = dat, psi = psi.huber) # Huber
fit_bis <- rlm(y ~ x1 + x2, data = dat, psi = psi.bisquare) # Tukey bisquare
bind_rows(
broom::tidy(fit_ols) %>% mutate(model = "OLS"),
broom::tidy(fit_hub) %>% mutate(model = "Huber"),
broom::tidy(fit_bis) %>% mutate(model = "Bisquare")
) %>%
arrange(term, model)## # A tibble: 9 × 6
## term estimate std.error statistic p.value model
## <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 (Intercept) 2.82 0.353 8.00 NA Bisquare
## 2 (Intercept) 3.22 0.360 8.92 NA Huber
## 3 (Intercept) 5.53 0.926 5.97 2.62e- 8 OLS
## 4 x1 2.04 0.0532 38.2 NA Bisquare
## 5 x1 1.95 0.0543 35.9 NA Huber
## 6 x1 1.44 0.140 10.3 4.24e-18 OLS
## 7 x2 -1.64 0.215 -7.63 NA Bisquare
## 8 x2 -1.51 0.220 -6.88 NA Huber
## 9 x2 -0.605 0.564 -1.07 2.86e- 1 OLSPerbandingan garis regresi (kasus sederhana y~x1):
fit_ols1 <- lm(y~ x1, dat)
fit_hub1 <- rlm(y~ x1, dat, psi=psi.huber)
fit_bis1 <- rlm(y~ x1, dat, psi=psi.bisquare)
plot(dat$x1, dat$y, pch=19, main="OLS vs Robust (y ~ x1)", xlab="x1", ylab="y")
abline(fit_ols1, lwd=2, lty=2)
abline(fit_hub1, lwd=2)
abline(fit_bis1, lwd=2, lty=3)
legend("topleft", c("OLS","Huber","Bisquare"), lty=c(2,1,3), lwd=2, bty="n")
3.2.5.3 Quantile Regression (Median)
fit_q50 <- rq(y~ x1 + x2, tau=0.5, data=dat)
broom::tidy(fit_q50)## # A tibble: 3 × 5
## term estimate conf.low conf.high tau
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.22 2.61 4.16 0.5
## 2 x1 1.97 1.74 2.08 0.5
## 3 x2 -1.50 -1.68 -0.924 0.53.2.5.4 Ringkasan Diagnostik & Rekomendasi
dx <- tibble(
idx = 1:n,
rstd = rstd,
leverage = h,
cooks = cooks
)
flagged <- (abs(rstd) > 2) | (h > 2 * length(coef(fit_ols)) / n) | (cooks > 1)
dx %>%
filter(flagged) %>%
arrange(desc(abs(rstd)))## # A tibble: 7 × 4
## idx rstd leverage cooks
## <int> <dbl> <dbl> <dbl>
## 1 55 -5.82 0.0496 0.589
## 2 95 -5.15 0.315 4.07
## 3 80 -4.64 0.0656 0.504
## 4 10 4.12 0.0151 0.0869
## 5 58 -0.768 0.0999 0.0218
## 6 68 0.725 0.0571 0.0106
## 7 44 -0.157 0.0528 0.000456Tahapan Analisis
- Eksplorasi grafis: scatter/QQ/residual plots.
- Hitung metrik: \(r_i\), \(h_{ii}\), \(D_i\), DFFITS/DFBETAS.
- Verifikasi data: benahi kesalahan input.
- Modelisasi: bandingkan OLS vs robust
(Huber/Bisquare/LAD/Quantile).
- Sensitivitas: laporkan perubahan koefisien/prediksi
setelah menandai/menangani outlier.
- Dokumentasi: catat kriteria, ambang, dan alasan keputusan.
Template Pelaporan - Data &
konteks: sumber, pembersihan, asumsi.
- Identifikasi outlier: metrik yang digunakan dan titik
yang ditandai.
- Penanganan: transformasi, winsorizing, atau model
robust; alasan pemilihan.
- Hasil: koefisien, interval kepercayaan/kuantil,
akurasi/predictive performance.
- Uji sensitivitas: perbandingan OLS vs robust.
- Catatan keterbatasan: potensi bias, ukuran sampel,
generalisasi.
Catatan Praktis - Gunakan beberapa indikator
bersamaan; jangan andalkan satu metrik.
- Ambang rule-of-thumb (mis. \(|r_i|>2\), \(D_i>1\)) bersifat panduan, bukan aturan
keras.
- Outlier dapat menyimpan informasi penting; pahami konteks sebelum
menghapus.
- Robust \(\neq\) kebal terhadap semua
masalah; tetap cek spesifikasi model, heteroskedastisitas, dan
multikolinearitas.
Referensi - Rousseeuw, P. J., & Leroy, A. M.
(2003). Robust Regression and Outlier Detection. Wiley.
- Huber, P. J., & Ronchetti, E. M. (2009). Robust
Statistics. Wiley.
- Fox, J. (2016). Applied Regression Analysis and Generalized Linear
Models. Sage.
- Koenker, R. (2005). Quantile Regression. Cambridge University
Press.
3.3 Tugas
Petunjuk
- Bagian A dan B dapat dijawab di kertas/jurnal Anda.
- Bagian C memerlukan R. Semua kode disediakan dan dapat langsung dijalankan.
- Kunci jawaban ada di bagian akhir dokumen.
3.3.1 Bagian A — Konsep (Jawaban Singkat)
- Jelaskan apa itu outlier pada regresi linier. Bedakan response outlier, predictor outlier, leverage point, dan influential point.
- Mengapa OLS sensitif terhadap outlier? Jelaskan memakai fungsi kerugian OLS.
- Jelaskan apa yang dimaksud influence function 𝜓(𝑒) dalam robust regression dan bedanya pada OLS, Huber, dan Tukey’s bisquare.
- Bagaimana memilih nilai ambang 𝑐 pada Huber dan Tukey? Jelaskan peran ̂𝜎 (skala robust).
- Outlier pada regresi linier
- Outlier respons (Y-space): residu sangat
besar.
- Outlier prediktor (X-space): nilai prediktor jauh
dari pusat sebaran.
- Leverage point: titik dengan X ekstrem, berpotensi
memengaruhi kemiringan; diukur dengan \(h_{ii}\).
- Influential point: jika dihapus, koefisien berubah nyata; ukur dengan Cook’s distance, DFFITS, DFBETAS.
- Mengapa OLS sensitif terhadap outlier
OLS meminimalkan \(\sum_{i=1}^n e_i^2\). Karena dikuadratkan, residu besar mendapat bobot berlebih sehingga satu titik ekstrem bisa mendominasi hasil. Secara influence function, OLS punya \(\psi(e)=e\) yang tidak dibatasi (unbounded).
- Influence function \(\psi(e)\) pada robust
regression
\(\psi(e)=\dfrac{d}{de}\rho(e)\) mengukur pengaruh marjinal residu terhadap estimasi.
OLS: \(\psi(e)=e\) (tidak dibatasi).
Huber: (e)=
\[\begin{cases} e,& |e|\le c\\ c\,\mathrm{sign}(e),& |e|>c \end{cases}\]Pengaruh outlier dibatasi pada \(\pm c\).
Tukey’s bisquare:
(e)=
\[\begin{cases} e\left(1-(e/c)^2\right)^2,& |e|\le c\\ 0,& |e|>c \end{cases}\]Outlier jauh hampir tidak berpengaruh.
- Memilih ambang \(c\) dan
peran \(\hat\sigma\)
- Huber: \(c \approx
1.345\,\hat\sigma\) (efisiensi $$95% saat error normal).
- Tukey: \(c \approx
4.685\,\hat\sigma\) (efisiensi mendekati OLS saat normal).
\(\hat\sigma\) adalah skala robust (mis. MAD atau pembaruan IRLS), sehingga \(c\) terkalibrasi dengan tingkat noise data: tetap efisien saat data normal, namun membatasi pengaruh outlier.
3.3.2 Bagian B — Analisis (Teori)
5) Bobot Huber dan identifikasi outlier
Diberikan residual OLS: \(e = \{-0.5, 1.2, 5.0\}\) dan \(c = 1.345\).
- Hitung bobot \(w_i\) untuk Huber loss.
- Observasi mana yang dianggap outlier oleh pendekatan ini?
6) Cook’s Distance Diberikan Cook’s Distance untuk tiga observasi: \(\{0.2, 0.8, 1.5\}\).
- Mana yang paling berpengaruh?
- Mengapa Cook’s Distance berguna dalam deteksi observasi berpengaruh?
- Bobot Huber dan identifikasi outlier — analisis teori
Diberikan: residual OLS \(e=\{-0.5,\;1.2,\;5.0\}\), konstanta Huber \(c=1.345\).
Definisi Huber (skor/influence) dan bobot IRLS: \[ \psi(e)= \begin{cases} e, & |e|\le c\\ c\,\mathrm{sign}(e), & |e|>c \end{cases} \qquad\text{dan}\qquad w(e)=\dfrac{\psi(e)}{e} \;=\; \begin{cases} 1, & |e|\le c\\ \dfrac{c}{|e|}, & |e|>c \end{cases} \]
Perhitungan bobot: - \(e=-0.5\): \(|e|=0.5 < c \Rightarrow w=1\). - \(e=1.2\): \(|e|=1.2 < c \Rightarrow w=1\). - \(e=5.0\): \(|e|=5.0 > c \Rightarrow w=c/|e|=1.345/5=0.269\).
Kesimpulan outlier (menurut ambang \(|e|>c\)): observasi dengan \(e=5.0\) diklasifikasikan outlier dan otomatis di-downweight (bobot 0.269), sehingga pengaruhnya pada estimasi menjadi terbatas.
- Cook’s Distance — analisis teori
Diberikan: nilai Cook’s Distance \(\{0.2,\;0.8,\;1.5\}\).
- Paling berpengaruh: \(D=1.5\) (terbesar). Aturan praktis: \(D_i>1\) sering dianggap sangat berpengaruh; alternatifnya \(D_i>4/n\) (bergantung ukuran sampel).
- Mengapa berguna: Cook’s Distance menggabungkan dua
komponen kunci:
- Residual terstandar besar (deviasi di \(Y\)-space), dan
- Leverage tinggi \(h_{ii}\) (ekstrem di \(X\)-space).
Secara teori, \(D_i\) mengaproksimasi perubahan keseluruhan fit/koefisien bila observasi \(i\) dihapus. Karena langsung mengukur dampak pada model (bukan hanya “keunikan” di \(X\) atau besar residual di \(Y\) saja), \(D_i\) efektif untuk menandai influential points yang benar-benar dapat mengubah hasil regresi.
- Residual terstandar besar (deviasi di \(Y\)-space), dan
3.3.3 Bagian C — Praktik dengan R
C.1 Simulasi data dengan outlier
set.seed(123)
x <- 1:20
y <- 2*x + rnorm(20, 0, 2)
y[c(5, 15)] <- y[c(5, 15)] + 20 # tambahkan outlier
dat <- data.frame(y, x)Estimasi OLS vs Robust (Huber)
library(MASS) # rlm
fit_ols <- lm(y~ x, dat)
fit_rob <- rlm(y~ x, dat, psi = psi.huber)
coef(fit_ols)## (Intercept) x
## 2.935974 1.937836coef(fit_rob)## (Intercept) x
## 0.9707001 1.9844546C.3 Visualisasi garis OLS vs Robust
plot(dat$x, dat$y, pch=19, xlab="x", ylab="y", main="OLS vs Robust (Huber)")
abline(fit_ols, lwd=2, lty=2, col="gray40")
abline(fit_rob, lwd=2, col="black")
legend("topleft", c("OLS","Robust (Huber)"),
lty=c(2,1), lwd=2, bty="n")
C.4 Diagnostik standar (leverage, residual, Cook’s D)
h_ii <- hatvalues(fit_ols)
rstd <- rstandard(fit_ols)
cooks <- cooks.distance(fit_ols)
idx_flag <- which( (abs(rstd)>2) | (h_ii > 2*length(coef(fit_ols))/nrow(dat)) | (cooks > 1) )
data.frame(idx=1:nrow(dat), rstd=round(rstd,3), h_ii=round(h_ii,3), cooks=round(cooks,3),
flagged = 1:nrow(dat) %in% idx_flag)[idx_flag, ]## idx rstd h_ii cooks flagged
## 5 5 2.906 0.095 0.446 TRUE
## 15 15 2.760 0.080 0.333 TRUEC.5 Bobot IRLS dari model robust
w <- weights(fit_rob)
data.frame(idx=1:length(w), weight=round(w,3))[order(w), ][1:5, ] # tampilkan 5 bobot terkecil## idx weight
## 1 1 1
## 2 2 1
## 3 3 1
## 4 4 1
## 5 5 13.4 Quis
- Waktu pengerjaan: 15–20 menit.
- Tuliskan jawaban singkat dan jelas.
- Bagian Kunci Jawaban & Rubrik dapat disembunyikan saat dibagikan ke mahasiswa.
3.4.1 Bagian A — Pilihan Ganda (PG)
(Masing-masing 2 poin, total 10 poin)
A1. Mengapa OLS sangat sensitif terhadap outlier?
A. Menggunakan fungsi absolute residual.
B. Menggunakan fungsi kuadrat residual.
C. Menggunakan median residual.
D. Tidak menggunakan residual.
A2. Dalam regresi, leverage point adalah …
A. Observasi dengan nilai Y ekstrem.
B. Observasi dengan nilai X jauh dari pusat
distribusi.
C. Observasi dengan \(|r_i| >
2\).
D. Observasi dengan Cook’s \(D >
1\).
A3. Jika residual \(e_i = 5\) dan ambang Huber \(c = 1.345\), maka bobot \(w_i\) kira-kira …
A. 1
B. 0.269
C. 5
D. 1.345
A4. Pada Tukey’s bisquare, apa yang terjadi jika \(|e| > c\)?
A. Residual dihitung penuh.
B. Residual diabaikan (bobot nol).
C. Residual diganti dengan rata-rata.
D. Residual dibagi dua.
A5. Estimasi skala \(\hat{\sigma}\) yang tahan outlier biasanya menggunakan …
A. Rata-rata residual.
B. SD residual OLS.
C. Median Absolute Deviation (MAD).
D. Variansi OLS.
3.4.2 Bagian B — Isian Singkat
(Masing-masing 5 poin, total 15 poin)
B6. Tuliskan formula Cook’s Distance dan jelaskan arti jika \(D_i > 1\).
B7. Sebutkan dua metode robust regression selain Huber loss.
B8. Diberikan residual \(\{-0.5, 1.2, 5.0\}\) dengan \(c = 1.345\). Tentukan residual mana yang dianggap outlier oleh metode Huber dan berikan alasannya.
Jawab
B6. Cook’s Distance
Rumus: \[
D_i \;=\; \frac{r_i^2}{p}\,\frac{h_{ii}}{1-h_{ii}}
\] dengan \(r_i\) residual
terstandar, \(h_{ii}\) leverage, dan
\(p\) jumlah parameter (termasuk
intersep).
Arti \(D_i > 1\):
secara praktik sering dianggap sangat berpengaruh;
observasi \(i\) berpotensi mengubah
koefisien/fit jika dihapus.
B7. Dua metode robust regression selain Huber
- Tukey’s bisquare (bisquare M-estimator).
- Least Absolute Deviations (LAD) / regresi L1
(alternatif: Quantile Regression, MM-estimator, RANSAC).
B8. Huber outlier check
Diberikan residual \(\{-0.5,\, 1.2,\,
5.0\}\) dan \(c=1.345\).
Kriteria Huber: titik dianggap outlier bila \(|e| > c\).
- \(|-0.5|=0.5 < 1.345 \Rightarrow\)
bukan outlier (bobot \(=1\)).
- \(|1.2|=1.2 < 1.345 \Rightarrow\)
bukan outlier (bobot \(=1\)).
- \(|5.0|=5.0 > 1.345 \Rightarrow\)
outlier (bobot \(\approx
c/|e|=1.345/5=0.269\), pengaruh dibatasi).