knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, dpi = 150,
fig.align = "center", fig.width = 7, fig.height = 4.5)
library(tidyverse)
library(broom)
library(sandwich) # robust/sandwich SE
library(lmtest) # coeftest
library(car) # VIF
library(modelr) # CV utilities
library(splines) # basis spline
set.seed(123)
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:
Model regresi linear sederhana dituliskan sebagai:
\[ Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i, \quad i=1,\ldots,n \]
dengan
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}} \]
Perbedaan Utama
| 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 Aplikasi di R
Data
Gunakan dataset mtcars untuk melihat hubungan
mpg (miles per gallon) dengan wt (berat
mobil).
data(mtcars)
head(mtcars[, c("mpg", "wt")])
## mpg wt
## Mazda RX4 21.0 2.620
## Mazda RX4 Wag 21.0 2.875
## Datsun 710 22.8 2.320
## Hornet 4 Drive 21.4 3.215
## Hornet Sportabout 18.7 3.440
## Valiant 18.1 3.460
Korelasi
cor.test(mtcars$wt, mtcars$mpg)
##
## Pearson's product-moment correlation
##
## data: mtcars$wt and mtcars$mpg
## t = -9.559, df = 30, p-value = 1.294e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9338264 -0.7440872
## sample estimates:
## cor
## -0.8676594
Regresi
model <- lm(mpg ~ wt, data = mtcars)
summary(model)
##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10
Visualiasi
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)")
Kita memodelkan hubungan antara respons \(y \in \mathbb{R}^n\) dan kovariat \(X \in \mathbb{R}^{n\times p}\):
\[ y = X\beta + \varepsilon, \qquad \mathbb{E}(\varepsilon)=0, \qquad \operatorname{Var}(\varepsilon)=\sigma^2 I_n.\ (\#eq:model) \]
Dengan \(X = [\mathbf{1}, x_1,\ldots,x_{p-1}]\) (termasuk intersep), \(\beta\) adalah vektor koefisien, dan \(\varepsilon\) adalah galat.
Estimator ordinary least squares (OLS) meminimalkan jumlah kuadrat residu:
\[ \widehat{\beta}_{\text{OLS}}=\arg\min_{\beta} \, \lVert y-X\beta \rVert^2 = (X^{\top}X)^{-1}X^{\top}y.\ (\#eq:ols) \]
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 = Hy, \qquad e = (I-H)y. \]
Leverage pengamatan ke-\(i\) adalah \(h_{ii}\) (elemen diagonal \(H\)).
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}.\ (\#eq:varbeta) \]
Asimtotik: \(\sqrt{n}\,(\widehat{\beta}-\beta) \overset{d}{\to} \mathcal{N}(0, V)\), sehingga inferensi berbasis normal/t berlaku untuk \(n\) besar.
Jika galat tidak homoskedastik/berkorelasi, maka \(\operatorname{Var}(\varepsilon)=\sigma^2\,\Omega \neq \sigma^2 I\).
Weighted least squares (WLS) untuk varians berbeda-beda (\(\Omega\) diagonal):
\[ \widehat{\beta}_{\text{WLS}}=(X^{\top}WX)^{-1}X^{\top}Wy,\qquad W=\Omega^{-1}.\ (\#eq:wls) \]
Generalized least squares (GLS) (\(\Omega\) umum):
\[ \widehat{\beta}_{\text{GLS}}=(X^{\top}\Omega^{-1}X)^{-1}X^{\top}\Omega^{-1}y.\ (\#eq:gls) \]
Tanpa menspesifikasi bentuk \(\Omega\), kovarians heteroskedasticity-consistent (HC) ala White adalah:
\[ \widehat{\operatorname{Var}}(\widehat{\beta})=(X^{\top}X)^{-1} \, \Big(\sum_{i=1}^{n} x_i x_i^{\top} \, \hat e_i^{2}\Big) \, (X^{\top}X)^{-1}.\ (\#eq:white) \]
Ini memberi robust SE yang valid terhadap heteroskedastisitas (HC0–HC3).
Variance Inflation Factor (VIF) untuk koefisien ke-\(j\):
\[ \mathrm{VIF}_j = \frac{1}{1-R_j^2}.\ (\#eq:vif) \]
dengan \(R_j^2\) adalah koefisien determinasi dari regresi \(x_j\) pada kovariat lain. Nilai besar (\(\gtrsim 10\)) mengindikasikan multikolinearitas kuat.
Ukuran pengaruh Cook’s distance:
\[ D_i = \frac{e_i^2}{p\,\widehat{\sigma}^{2}}\cdot\frac{h_{ii}}{(1-h_{ii})^{2}}.\ (\#eq:cook) \]
Gunakan plot residu, plot leverage vs residu, QQ-plot, serta inspeksi \(D_i\) untuk memeriksa asumsi dan outlier.
Model makin kompleks (banyak parameter/basis) cenderung bias rendah namun varians tinggi. Ridge dan Lasso menambahkan penalti:
\[ \text{Ridge:}\quad \widehat{\beta}=\arg\min_{\beta} \lVert y-X\beta \rVert^{2}+\lambda\lVert\beta\rVert_2^{2},\\ \text{Lasso:}\quad \widehat{\beta}=\arg\min_{\beta} \lVert y-X\beta \rVert^{2}+\lambda\lVert\beta\rVert_1.\ (\#eq:ridge-lasso) \]
Kita dapat menambah term polinomial atau basis spline sehingga tetap linear dalam parameter:
\[ y = \beta_0 + \sum_{k} \beta_k\, b_k(x) + \varepsilon.\ (\#eq:spline) \]
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
fit1 <- lm(mpg ~ wt + hp + disp, data = dat1)
summary(fit1)
##
## Call:
## lm(formula = mpg ~ wt + hp + disp, data = dat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.891 -1.640 -0.172 1.061 5.861
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.105505 2.110815 17.579 < 2e-16 ***
## wt -3.800891 1.066191 -3.565 0.00133 **
## hp -0.031157 0.011436 -2.724 0.01097 *
## disp -0.000937 0.010350 -0.091 0.92851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.639 on 28 degrees of freedom
## Multiple R-squared: 0.8268, Adjusted R-squared: 0.8083
## F-statistic: 44.57 on 3 and 28 DF, p-value: 8.65e-11
# Koefisien dengan CI 95%
tidy(fit1, conf.int = TRUE)
## # A tibble: 4 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 37.1 2.11 17.6 1.16e-16 32.8 41.4
## 2 wt -3.80 1.07 -3.56 1.33e- 3 -5.98 -1.62
## 3 hp -0.0312 0.0114 -2.72 1.10e- 2 -0.0546 -0.00773
## 4 disp -0.000937 0.0103 -0.0905 9.29e- 1 -0.0221 0.0203
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
Diagnostik OLS untuk mtcars: residu vs fitted (kiri) dan QQ-plot (kanan).
Diagnostik OLS untuk mtcars: residu vs fitted (kiri) dan QQ-plot (kanan).
qplot(seq_along(diag_df$cook), diag_df$cook, geom = c("line","point"),
xlab = "Index", ylab = "Cook's D")
Cook’s distance per observasi.
car::vif(fit1)
## wt hp disp
## 4.844618 2.736633 7.324517
coeftest(fit1) # SE OLS
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.10550527 2.11081525 17.5788 < 2.2e-16 ***
## wt -3.80089058 1.06619064 -3.5649 0.001331 **
## hp -0.03115655 0.01143579 -2.7245 0.010971 *
## disp -0.00093701 0.01034974 -0.0905 0.928507
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit1, vcov = vcovHC(fit1, type = "HC3")) # SE robust (HC3)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.10550527 2.71026981 13.6907 6.257e-14 ***
## wt -3.80089058 1.08086207 -3.5165 0.00151 **
## hp -0.03115655 0.01437457 -2.1675 0.03886 *
## disp -0.00093701 0.01012138 -0.0926 0.92690
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
wtfit_poly <- lm(mpg ~ poly(wt, 2, raw = TRUE) + hp + disp, data = dat1)
fit_spl <- lm(mpg ~ bs(wt, df = 5) + hp + disp, data = dat1)
AIC(fit1, fit_poly, fit_spl) %>% arrange(AIC)
## df AIC
## fit_poly 6 150.7305
## fit_spl 9 154.3323
## fit1 5 158.6430
BIC(fit1, fit_poly, fit_spl) %>% arrange(BIC)
## df BIC
## fit_poly 6 159.5250
## fit1 5 165.9717
## fit_spl 9 167.5239
# Kurva prediksi
grid_wt <- tibble(wt = seq(min(dat1$wt), max(dat1$wt), length.out = 200),
hp = mean(dat1$hp), disp = mean(dat1$disp))
pred <- grid_wt %>%
mutate(OLS = predict(fit1, newdata = grid_wt),
Poly = predict(fit_poly, newdata = grid_wt),
Spline = predict(fit_spl, newdata = grid_wt)) %>%
pivot_longer(OLS:Spline, names_to = "model", values_to = "yhat")
GG2 <- ggplot() +
geom_point(data = dat1, aes(wt, mpg), alpha = 0.3) +
geom_line(data = pred, aes(wt, yhat, linetype = model), linewidth = 0.9) +
labs(x = "wt", y = "mpg", title = "Perbandingan OLS vs Polinomial vs Spline")
GG2
k <- 5
cv_res <- crossv_kfold(dat1, k = k)
rmse_fun <- function(train, test, formula) {
m <- lm(formula, data = as.data.frame(train))
sqrt(mean((as.data.frame(test)$mpg - predict(m, newdata = as.data.frame(test)))^2))
}
forms <- list(
OLS = mpg ~ wt + hp + disp,
Poly = mpg ~ poly(wt, 2, raw = TRUE) + hp + disp,
Spline = mpg ~ bs(wt, df = 5) + hp + disp
)
cv_tbl <- map_dfr(names(forms), function(nm){
fml <- forms[[nm]]
tibble(model = nm, fold = seq_len(k),
rmse = map2_dbl(cv_res$train, cv_res$test, ~ rmse_fun(.x, .y, fml)))
})
cv_tbl %>% group_by(model) %>% summarise(mean_rmse = mean(rmse), sd_rmse = sd(rmse)) %>% arrange(mean_rmse)
## # A tibble: 3 × 3
## model mean_rmse sd_rmse
## <chr> <dbl> <dbl>
## 1 Poly 2.35 0.980
## 2 Spline 2.74 0.807
## 3 OLS 2.92 1.01
n <- 400
x1 <- runif(n, 0, 10)
x2 <- rnorm(n, 5, 2)
sigma_i <- 0.6 + 0.2 * x1 # var naik dengan x1
err <- rnorm(n, 0, sigma_i)
y <- 2 + 1.5*x1 - 1.2*x2 + err
sim <- tibble(y, x1, x2, w = 1/sigma_i^2)
fit_ols <- lm(y ~ x1 + x2, data = sim)
fit_wls <- lm(y ~ x1 + x2, data = sim, weights = w)
# Tabel perbandingan SE
comp_se <- tibble(
term = names(coef(fit_ols)),
OLS_SE = coef(summary(fit_ols))[, "Std. Error"],
Robust_SE = sqrt(diag(vcovHC(fit_ols, type = "HC3"))),
WLS_SE = coef(summary(fit_wls))[, "Std. Error"]
)
comp_se
## # A tibble: 3 × 4
## term OLS_SE Robust_SE WLS_SE
## <chr> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.250 0.212 0.175
## 2 x1 0.0293 0.0282 0.0261
## 3 x2 0.0415 0.0399 0.0313
# Uji Breusch-Pagan untuk heteroskedastisitas
lmtest::bptest(fit_ols)
##
## studentized Breusch-Pagan test
##
## data: fit_ols
## BP = 29.415, df = 2, p-value = 4.098e-07
ggplot(data.frame(fitted=fitted(fit_ols), resid=resid(fit_ols)), aes(fitted, resid)) +
geom_point(alpha=.6) + geom_smooth(se=FALSE, method="loess") +
geom_hline(yintercept = 0, linetype=2) + labs(x="Fitted", y="Residuals")
Heteroskedastisitas: Residual vs Fitted (OLS).
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} = \sigma_u^2\,\rho^{|t-s|}/(1-\rho^2)\). GLS menimbang data dengan \(\Omega^{-1}\) (atau transformasi ekuivalen) agar estimasi efisien.
set.seed(42)
T <- 400
x <- runif(T, -2, 2)
rho <- 0.6
u <- rnorm(T, 0, 1)
err <- numeric(T)
for (t in 2:T) err[t] <- rho*err[t-1] + u[t]
y <- 1 + 2*x + err
dar1 <- tibble::tibble(t = 1:T, x, y)
head(dar1)
## # A tibble: 6 × 3
## t x y
## <int> <dbl> <dbl>
## 1 1 1.66 4.32
## 2 2 1.75 4.83
## 3 3 -0.855 0.661
## 4 4 1.32 6.53
## 5 5 0.567 2.49
## 6 6 0.0764 0.213
library(nlme)
# OLS biasa
fit_ols_ar1 <- lm(y ~ x, data = dar1)
# GLS dengan korelasi AR(1) di residual
fit_gls_ar1 <- gls(y ~ x, data = dar1,
correlation = corAR1(form = ~ t))
summary(fit_ols_ar1)
##
## Call:
## lm(formula = y ~ x, data = dar1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3360 -0.8231 0.0026 0.8211 3.9013
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.95441 0.06043 15.79 <2e-16 ***
## x 2.06045 0.05162 39.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.209 on 398 degrees of freedom
## Multiple R-squared: 0.8002, Adjusted R-squared: 0.7997
## F-statistic: 1594 on 1 and 398 DF, p-value: < 2.2e-16
summary(fit_gls_ar1)
## Generalized least squares fit by REML
## Model: y ~ x
## Data: dar1
## AIC BIC logLik
## 1139.541 1155.487 -565.7706
##
## Correlation Structure: AR(1)
## Formula: ~t
## Parameter estimate(s):
## Phi
## 0.5772332
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.9552097 0.11651869 8.19791 0
## x 2.0449090 0.03743358 54.62766 0
##
## Correlation:
## (Intr)
## x 0.004
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.752477543 -0.674750151 -0.007333444 0.679546648 3.225627947
##
## Residual standard error: 1.210605
## Degrees of freedom: 400 total; 398 residual
par(mfrow = c(1,2))
acf(resid(fit_ols_ar1), main = "ACF Residuals (OLS)")
acf(resid(fit_gls_ar1, type = "normalized"), main = "ACF Residuals (GLS)")
ACF residu dari OLS vs GLS
par(mfrow = c(1,1))
lmtest::dwtest(fit_ols_ar1) # indikasi korelasi serial pada OLS
##
## Durbin-Watson test
##
## data: fit_ols_ar1
## DW = 0.84993, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0
# Estimasi rho dari regresi resid_t ~ resid_{t-1}
res <- resid(fit_ols_ar1)
rho_hat <- coef(lm(res[-1] ~ 0 + res[-length(res)]))[1]
# Quasi-differencing (transformasi data)
y_star <- dar1$y[-1] - rho_hat*dar1$y[-T]
x_star <- dar1$x[-1] - rho_hat*dar1$x[-T]
fit_qd <- lm(y_star ~ x_star)
summary(fit_qd)
##
## Call:
## lm(formula = y_star ~ x_star)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6691 -0.6723 -0.0232 0.6970 3.2363
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.40603 0.04955 8.194 3.52e-15 ***
## x_star 2.04499 0.03757 54.439 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9898 on 397 degrees of freedom
## Multiple R-squared: 0.8819, Adjusted R-squared: 0.8816
## F-statistic: 2964 on 1 and 397 DF, p-value: < 2.2e-16
Catatan:
gls()memperkirakan \(\rho\) dan ragam secara simultan dalam kerangka likelihood sehingga sering lebih stabil dibanding pendekatan dua tahap.
set.seed(123)
# Paket yang diperlukan
pkgs <- c("dplyr","ggplot2","car","lmtest","sandwich","MASS","boot","broom")
# Memuat semua paket tanpa menampilkan TRUE
invisible(lapply(pkgs, require, character.only = TRUE))
Tujuan: membuat data dengan 3 prediktor (X1, X2, X3), disuntik multikolinearitas ringan (X1 & X2 berkorelasi) dan heteroskedastisitas (ragam error meningkat saat |X1| besar), serta sedikit outlier.
n <- 200
X1 <- rnorm(n, 0, 1)
# X2 berkorelasi ~0.7 dengan X1
rho <- 0.7
X2 <- rho*X1 + sqrt(1-rho^2)*rnorm(n)
X3 <- rnorm(n)
# Heteroskedastisitas: sd error tergantung |X1|
eps <- rnorm(n, sd = 0.5 + 0.5*abs(X1))
# Model benar (ground truth): y = 2 + 1.5*X1 - 0.8*X2 + 0.5*X3 + eps
y <- 2 + 1.5*X1 - 0.8*X2 + 0.5*X3 + eps
# Tambah beberapa outlier pada respon
idx_out <- sample(1:n, 3)
y[idx_out] <- y[idx_out] + rnorm(3, mean = 6, sd = 1)
dat <- data.frame(y, X1, X2, X3)
summary(dat)
## y X1 X2 X3
## Min. :-3.194 Min. :-2.30917 Min. :-2.36565 Min. :-2.80977
## 1st Qu.: 1.157 1st Qu.:-0.62576 1st Qu.:-0.68926 1st Qu.:-0.55753
## Median : 2.040 Median :-0.05874 Median : 0.03268 Median : 0.07583
## Mean : 2.066 Mean :-0.00857 Mean : 0.02408 Mean : 0.03178
## 3rd Qu.: 2.895 3rd Qu.: 0.56840 3rd Qu.: 0.69036 3rd Qu.: 0.68098
## Max. : 9.427 Max. : 3.24104 Max. : 3.00049 Max. : 2.43023
mod0 <- lm(y ~ X1 + X2 + X3, data = dat)
summary(mod0)
##
## Call:
## lm(formula = y ~ X1 + X2 + X3, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3070 -0.6446 -0.0414 0.5284 7.1821
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.07742 0.08726 23.806 < 2e-16 ***
## X1 1.31536 0.12466 10.551 < 2e-16 ***
## X2 -0.75833 0.12303 -6.164 3.97e-09 ***
## X3 0.57737 0.09070 6.366 1.35e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.232 on 196 degrees of freedom
## Multiple R-squared: 0.4404, Adjusted R-squared: 0.4318
## F-statistic: 51.41 on 3 and 196 DF, p-value: < 2.2e-16
broom::glance(mod0) # ringkas: R^2, Adj R^2, AIC, dll.
## # A tibble: 1 × 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.440 0.432 1.23 51.4 1.47e-24 3 -324. 657. 674.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>
Catatan interpretasi cepat: Koefisien menunjukkan perubahan rata-rata y untuk kenaikan 1 unit pada prediktor (dengan prediktor lain konstan). Lihat Pr(>|t|) untuk signifikansi.
3.1 Uji serentak (overall F-test)
Sudah tersedia pada output summary(mod0) (Signifikansi model secara keseluruhan).
3.2 Uji parsial (t-test per koefisien)
Juga ada di summary(mod0).
3.3 Uji hipotesis gabungan (contoh:
car::linearHypothesis(mod0, c("X2 = 0", "X3 = 0"))
##
## Linear hypothesis test:
## X2 = 0
## X3 = 0
##
## Model 1: restricted model
## Model 2: y ~ X1 + X2 + X3
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 198 423.38
## 2 196 297.59 2 125.79 41.423 9.897e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
3.4 Uji dengan SE robust (menghadapi heteroskedastisitas)
# Koefisien + SE robust (HC3)
coeftest(mod0, vcov. = sandwich::vcovHC(mod0, type = "HC3"))
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.077419 0.087811 23.6578 < 2.2e-16 ***
## X1 1.315357 0.163042 8.0676 7.027e-14 ***
## X2 -0.758334 0.123626 -6.1341 4.644e-09 ***
## X3 0.577370 0.107728 5.3595 2.330e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 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.263679
Rule of thumb: VIF > 10 (atau > 5) → indikasi kolinearitas meningkat.
4.2 Heteroskedastisitas
# Breusch-Pagan
lmtest::bptest(mod0)
##
## studentized Breusch-Pagan test
##
## data: mod0
## BP = 1.0283, df = 3, p-value = 0.7944
# White test (via bptest dengan kuadrat fitted)
lmtest::bptest(mod0, ~ fitted(mod0) + I(fitted(mod0)^2))
##
## studentized Breusch-Pagan test
##
## data: mod0
## BP = 5.3586, df = 2, p-value = 0.06861
Jika signifikan, gunakan SE robust (lihat 3.4) untuk inferensi yang lebih andal.
4.3 Normalitas residual
res <- rstandard(mod0)
shapiro.test(res) # sensitif untuk n besar; gunakan juga QQ-plot
##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.87436, p-value = 7.729e-12
qqnorm(res); qqline(res)
4.4 Autokorelasi residual (opsional)
(Biasanya untuk data runtun waktu)
lmtest::dwtest(mod0)
##
## Durbin-Watson test
##
## data: mod0
## DW = 2.0342, p-value = 0.5998
## alternative hypothesis: true autocorrelation is greater than 0
4.5 Uji spesifikasi (linearitas) – Ramsey RESET
lmtest::resettest(mod0)
##
## RESET test
##
## data: mod0
## RESET = 4.2735, df1 = 2, df2 = 194, p-value = 0.01527
4.6 Plot Diagnostik Standar
par(mfrow=c(2,2))
plot(mod0)
par(mfrow=c(1,1))
Strategi: mulai dari model kandidat (termasuk transformasi sederhana) lalu gunakan AIC (stepwise), bandingkan dengan model awal, dan validasi via CV.
dat2 <- dat %>%
mutate(X1_sq = X1^2, X2_sq = X2^2, X3_sq = X3^2,
X1X2 = X1*X2, X1X3 = X1*X3, X2X3 = X2*X3)
mod_full <- lm(y ~ X1 + X2 + X3 + X1_sq + X2_sq + X3_sq + X1X2 + X1X3 + X2X3, data = dat2)
AIC(mod0); AIC(mod_full)
## [1] 657.0574
## [1] 660.7653
# Stepwise berbasis AIC (dua arah)
mod_step <- MASS::stepAIC(mod_full, direction = "both", trace = FALSE)
summary(mod_step)
##
## Call:
## lm(formula = y ~ X1 + X2 + X3 + X2_sq + X1X2, data = dat2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2286 -0.6853 -0.0298 0.4825 7.0242
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.19439 0.10487 20.925 < 2e-16 ***
## X1 1.31032 0.12513 10.472 < 2e-16 ***
## X2 -0.72799 0.12278 -5.929 1.37e-08 ***
## X3 0.57331 0.08996 6.373 1.32e-09 ***
## X2_sq -0.28501 0.11895 -2.396 0.0175 *
## X1X2 0.23671 0.13546 1.748 0.0821 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.22 on 194 degrees of freedom
## Multiple R-squared: 0.4571, Adjusted R-squared: 0.4431
## F-statistic: 32.67 on 5 and 194 DF, p-value: < 2.2e-16
broom::glance(mod_step)[,c("r.squared","adj.r.squared","AIC","BIC")]
## # A tibble: 1 × 4
## r.squared adj.r.squared AIC BIC
## <dbl> <dbl> <dbl> <dbl>
## 1 0.457 0.443 655. 678.
Validasi silang (10-fold CV) untuk membandingkan RMSE:
# Fungsi CV RMSE menggunakan boot::cv.glm
cv_rmse <- function(model, data, K=10){
# cv.glm mengembalikan delta: [1] raw CV, [2] adjusted CV
set.seed(123)
res <- boot::cv.glm(data = data, glmfit = model, K = K)
sqrt(res$delta[1])
}
rmse_mod0 <- cv_rmse(mod0, dat)
rmse_modstep<- cv_rmse(mod_step, dat2)
c(RMSE_mod0 = rmse_mod0, RMSE_modStep = rmse_modstep)
## RMSE_mod0 RMSE_modStep
## NaN NaN
Pilih model dengan trade-off terbaik antara kesederhanaan dan kinerja (RMSE, AIC/BIC, Adj, R-squared.
Gunakan studentized residuals, leverage (hat), dan Cook’s distance. Tambahkan Bonferroni outlier test.
infl <- influence.measures(mod0)
# Ringkasan pengaruh
summary(infl)
## Potentially influential observations of
## lm(formula = y ~ X1 + X2 + X3, data = dat) :
##
## dfb.1_ dfb.X1 dfb.X2 dfb.X3 dffit cov.r cook.d hat
## 16 -0.21 -0.06 -0.32 0.54 -0.78_* 0.94 0.15 0.07_*
## 18 0.30 -0.30 -0.24 0.09 0.74_* 0.75_* 0.13 0.03
## 27 0.00 0.00 0.00 0.01 0.01 1.07_* 0.00 0.04
## 56 -0.05 -0.05 0.01 0.11 -0.14 1.07_* 0.00 0.06
## 65 -0.14 0.36 -0.34 0.17 -0.45_* 0.97 0.05 0.04
## 90 0.31 0.25 0.04 -0.12 0.50_* 0.72_* 0.06 0.01
## 97 0.00 0.00 -0.01 0.00 -0.01 1.08_* 0.00 0.05
## 135 -0.03 0.07 -0.01 -0.08 -0.12 1.07_* 0.00 0.05
## 147 0.19 -0.21 -0.01 -0.25 0.42 0.92_* 0.04 0.03
## 149 0.23 0.54 -0.20 0.60 0.86_* 0.89_* 0.18 0.07_*
## 160 -0.02 0.04 -0.04 0.01 -0.05 1.06_* 0.00 0.04
## 162 0.43 -0.56 0.28 0.26 0.78_* 0.48_* 0.13 0.01
## 164 -0.20 -0.38 -0.25 -0.35 -0.87_* 0.95 0.18 0.08_*
## 191 0.03 0.02 -0.02 -0.08 0.08 1.07_* 0.00 0.05
# Statistik diagnostik
studres <- rstudent(mod0)
lev <- hatvalues(mod0)
cookd <- cooks.distance(mod0)
# Ambang sederhana
thr_res <- 3 # |studentized residual| > 3
thr_lev <- 2*length(coef(mod0))/nrow(dat) # leverage tinggi
thr_cook <- 4/nrow(dat) # Cook's distance besar
flag <- data.frame(
id = 1:nrow(dat),
studres = studres,
leverage = lev,
cookd = cookd
) %>%
mutate(
flag_res = abs(studres) > thr_res,
flag_lev = leverage > thr_lev,
flag_cook = cookd > thr_cook,
any_flag = flag_res | flag_lev | flag_cook
) %>%
arrange(desc(any_flag), desc(abs(studres)))
head(flag, 10)
## id studres leverage cookd flag_res flag_lev flag_cook any_flag
## 162 162 6.450812 0.01441863 0.12607179 TRUE FALSE TRUE TRUE
## 90 90 4.302679 0.01343365 0.05785187 TRUE FALSE TRUE TRUE
## 18 18 4.188058 0.03056365 0.12748735 TRUE FALSE TRUE TRUE
## 149 149 3.243249 0.06634585 0.17821039 TRUE TRUE TRUE TRUE
## 164 164 -2.855691 0.08454000 0.18164105 FALSE TRUE TRUE TRUE
## 16 16 -2.803489 0.07241620 0.14821120 FALSE TRUE TRUE TRUE
## 147 147 2.584024 0.02573283 0.04284913 FALSE FALSE TRUE TRUE
## 143 143 2.200166 0.03420638 0.04203837 FALSE FALSE TRUE TRUE
## 30 30 2.172800 0.03120105 0.03730329 FALSE FALSE TRUE TRUE
## 65 65 -2.140175 0.04228760 0.04965406 FALSE TRUE TRUE TRUE
Bonferroni outlier test (car):
car::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-03
Plot cepat untuk Cook’s distance:
plot(cookd, type="h", main="Cook's Distance", ylab="D", xlab="Observasi")
abline(h = thr_cook, lty=2)
Tindak lanjut (opsional):
Coba refit tanpa observasi berpengaruh ekstrem dan bandingkan koefisien/fit.
Pertimbangkan robust regression (M-estimator) jika banyak outlier: MASS::rlm.
# Refit tanpa titik yang sangat berpengaruh (contoh)
to_drop <- which(cookd > thr_cook | abs(studres) > thr_res)
mod_refit <- lm(y ~ X1 + X2 + X3, data = dat[-to_drop, ])
broom::tidy(mod0)
## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 2.08 0.0873 23.8 9.73e-60
## 2 X1 1.32 0.125 10.6 6.73e-21
## 3 X2 -0.758 0.123 -6.16 3.97e- 9
## 4 X3 0.577 0.0907 6.37 1.35e- 9
broom::tidy(mod_refit)
## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1.99 0.0566 35.1 2.69e-83
## 2 X1 1.44 0.0861 16.7 1.11e-38
## 3 X2 -0.676 0.0819 -8.26 2.87e-14
## 4 X3 0.501 0.0606 8.27 2.75e-14
mtcars: heteroskedastisitas (Breusch–Pagan),
multikolinearitas (VIF), pengaruh (Cook’s D), normalitas (QQ-plot).
Laporkan temuan & rekomendasi.mpg ~ wt dan bandingkan AIC/BIC + 5-fold CV.x dengan sinyal musiman (mis. sin(2*pi*t/12)),
interpretasikan \(\hat\rho\);
bandingkan SE OLS vs GLS.Regresi linear mengasumsikan hubungan antara variabel respon \(Y\) dan prediktor \(X\) berbentuk linear:
\[ Y = \beta_0 + \beta_1 X + \varepsilon, \quad \varepsilon \overset{iid}{\sim} N(0, \sigma^2_{\varepsilon}) \] Namun, dalam banyak fenomena nyata (biologi, farmasi, pertumbuhan, epidemiologi), hubungan tidak linear. Contoh:
Model
Model pertumbuhan bakteri seringkali mengikuti bentuk kurva logistik. Model matematisnya dapat dituliskan sebagai:
\[ Y(t) = \frac{K}{1 + \exp\{-(\alpha + \beta t)\}} + \varepsilon, \quad \varepsilon \sim N(0, \sigma^2) \]
Keterangan:
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:
Model
Model farmakokinetik Michaelis–Menten sering digunakan untuk menjelaskan hubungan antara konsentrasi obat (\(X\)) dan respon biologis (\(Y\)). Modelnya adalah:
\[ Y = \frac{V_{\text{max}} \, X}{K_m + X} + \varepsilon, \quad \varepsilon \sim N(0, \sigma^2) \]
Keterangan:
# 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:
Model
Contoh model sederhana dengan bentuk konkaf dapat dituliskan sebagai:
\[ Y = A X^{\alpha} + \varepsilon, \quad 0 < \alpha < 1, \quad \varepsilon \sim N(0, \sigma^2) \]
Alternatif lain:
\[ Y = a \, \log(1 + bX) + \varepsilon, \quad \varepsilon \sim N(0, \sigma^2) \]
Keterangan:
# 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:
Merujuk pada contoh diatas maka dibutuhkan regresi nonlinear.
Secara umum, model regresi nonlinear ditulis sebagai:
\[ Y_i = f(X_i, \boldsymbol{\theta}) + \varepsilon_i, \quad \varepsilon_i \sim N(0, \sigma^2) \]
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, \quad \varepsilon \sim N(0, \sigma^2) \]
2. Error Multiplikatif \[ Y = \beta_0 e^{\beta_1 X} \cdot \varepsilon, \quad \varepsilon > 0 \]
Ambil logaritma kedua sisi:
\[ \ln Y = \ln \beta_0 + \beta_1 X + \ln \varepsilon \]
Jika diasumsikan:
\[ \ln \varepsilon \sim N(0, \sigma^2), \]
maka model menjadi linear dalam parameter:
\[ \ln Y = \alpha_0 + \beta_1 X + u, \quad u \sim N(0, \sigma^2) \]
dengan definisi:
\[ \alpha_0 = \ln \beta_0 \]
Estimasi
\[ \hat{\beta}_0 = e^{\hat{\alpha}_0} \]
Ringkasan
Kita ingin menaksir parameter \(\boldsymbol\theta \in \mathbb{R}^p\) pada model
\[ Y_i = f(X_i,\boldsymbol\theta) + \varepsilon_i,\qquad \varepsilon_i \overset{iid}{\sim} \mathcal{N}(0,\sigma^2),\quad i=1,\dots,n. \]
Estimator nonlinear least squares (NLS) meminimalkan jumlah kuadrat residu:
\[ S(\boldsymbol\theta) = \sum_{i=1}^n r_i(\boldsymbol\theta)^2, \qquad r_i(\boldsymbol\theta)=y_i - f(X_i,\boldsymbol\theta). \]
Notasi Jacobian dari \(f\) terhadap parameter:
\[ J_{ij}(\boldsymbol\theta) = \frac{\partial f(x_i,\boldsymbol\theta)}{\partial \theta_j}, \quad J(\boldsymbol\theta)\in\mathbb{R}^{n\times p}. \]
Di sekitar tebakan \(\boldsymbol\theta^{(k)}\), pendekatan Taylor orde-1 memberi:
\[ \mathbf{r}(\boldsymbol\theta^{(k)}+\Delta) \;\approx\; \mathbf{r}^{(k)} - J^{(k)}\,\Delta, \quad \mathbf{r}^{(k)}= \mathbf{y}-\mathbf{f}(\boldsymbol\theta^{(k)}). \]
Gauss–Newton (GN) — Ide Inti
Dengan linearisasi di atas, langkah perbaikan \(\Delta\) dicari dari normal equations:
\[ (J^{\top}J)\,\Delta = J^{\top}\,\mathbf{r}. \]
Lalu diperbarui \(\boldsymbol\theta^{(k+1)}=\boldsymbol\theta^{(k)}+\Delta\).
Ulangi hingga konvergen (mis. \(\|\Delta\|\) 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}\,\mathbf{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
Kita gunakan model:
\[ Y = \beta_0\,e^{\beta_1 X} + \varepsilon,\qquad \varepsilon \sim \mathcal{N}(0,\sigma^2). \]
Untuk \(x_i\), \(f(x_i,\boldsymbol\theta)=\beta_0 e^{\beta_1 x_i}\), \(\boldsymbol\theta=(\beta_0,\beta_1)^\top\).
Jacobian \(J\) (terhadap \(\beta_0,\beta_1\)):
\[ \frac{\partial f}{\partial \beta_0} = e^{\beta_1 x_i}, \qquad \frac{\partial f}{\partial \beta_1} = \beta_0 x_i e^{\beta_1 x_i}. \] Simulasi Data & Visualisasi Awal
set.seed(123)
library(ggplot2)
library(dplyr)
library(tidyr)
beta0_true <- 2.0
beta1_true <- 0.3
sigma_true <- 0.4
# Small dataset agar mudah ditelusuri langkah demi langkah
x <- c(0, 1, 2, 3, 4, 5)
f_true <- beta0_true * exp(beta1_true * x)
y <- f_true + rnorm(length(x), sd = sigma_true)
dat <- tibble(x, y, f_true)
dat
## # A tibble: 6 × 3
## x y f_true
## <dbl> <dbl> <dbl>
## 1 0 1.78 2
## 2 1 2.61 2.70
## 3 2 4.27 3.64
## 4 3 4.95 4.92
## 5 4 6.69 6.64
## 6 5 9.65 8.96
Plot
ggplot(dat, aes(x, y)) +
geom_point(size=2) +
geom_line(aes(y = f_true), linewidth = 0.9) +
labs(title="Data Eksponensial (dengan noise) & Kurva Sejati",
x="X", y="Y") +
theme_minimal()
Gauss–Newton: Estimasi Manual (Step-by-step) Fungsi Bantu (f, residu, Jacobian)
f_exp <- function(x, theta){ # theta = c(beta0, beta1)
beta0 <- theta[1]; beta1 <- theta[2]
beta0 * exp(beta1 * x)
}
resid_vec <- function(x, y, theta){
y - f_exp(x, theta)
}
J_exp <- function(x, theta){
beta0 <- theta[1]; beta1 <- theta[2]
e <- exp(beta1 * x)
cbind( d_beta0 = e,
d_beta1 = beta0 * x * e )
}
Satu langkah iterasi Gauss–Newton (GN) diperoleh dengan menyelesaikan persamaan normal:
\[ (J^\top J)\,\Delta = J^\top r \]
gn_step <- function(x, y, theta){
r <- resid_vec(x, y, theta)
J <- J_exp(x, theta)
# Normal equations
lhs <- crossprod(J, J) # J^T J
rhs <- crossprod(J, r) # J^T r
Delta <- solve(lhs, rhs) # langkah GN
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) # tebakan awal sengaja "kasar"
step1 <- gn_step(x, y, theta0)
step1
## $Delta
## [1] 0.2542352 1.4943019
##
## $theta_new
## beta0 beta1
## 1.254235 1.494302
##
## $r
## [1] 0.7758097 1.6076466 3.2677209 3.9474096 5.6919489 8.6494041
##
## $J
## d_beta0 d_beta1
## [1,] 1 0
## [2,] 1 1
## [3,] 1 2
## [4,] 1 3
## [5,] 1 4
## [6,] 1 5
##
## $lhs
## d_beta0 d_beta1
## d_beta0 6 15
## d_beta1 15 55
##
## $rhs
## [,1]
## d_beta0 23.93994
## d_beta1 86.00013
##
## $SSE
## [1] 136.6569
Iterasi Beberapa Kali Hingga Konvergen
gn_fit <- function(x, y, theta_init, max_it=20, tol=1e-8){
th <- theta_init
hist <- tibble(iter=0, beta0=th[1], beta1=th[2],
SSE = sum(resid_vec(x,y,th)^2), step_norm = NA_real_)
for(k in 1:max_it){
st <- gn_step(x, y, th)
step_norm <- sqrt(sum(st$Delta^2))
th <- st$theta_new
hist <- add_row(hist, iter=k, beta0=th[1], beta1=th[2],
SSE = sum(resid_vec(x,y,th)^2), step_norm = step_norm)
if(step_norm < tol) break
}
list(theta = th, history = hist)
}
gn_res <- gn_fit(x, y, theta_init = theta0, max_it = 20)
gn_res$history
## # A tibble: 12 × 5
## iter beta0 beta1 SSE step_norm
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 1 0 137. NA
## 2 1 1.25 1.49 5064963. 1.52e+0
## 3 2 0.0813 1.48 16167. 1.17e+0
## 4 3 0.0851 1.29 1963. 1.95e-1
## 5 4 0.178 0.905 68.0 3.93e-1
## 6 5 0.664 0.270 98.5 8.00e-1
## 7 6 1.96 0.422 62.3 1.30e+0
## 8 7 1.89 0.347 2.26 1.04e-1
## 9 8 1.97 0.317 0.468 8.45e-2
## 10 9 1.98 0.315 0.465 8.30e-3
## 11 10 1.98 0.315 0.465 2.93e-5
## 12 11 1.98 0.315 0.465 3.99e-9
gn_res$theta
## beta0 beta1
## 1.9753264 0.3148085
Kurva Hasil GN
dat$yhat_gn <- f_exp(dat$x, gn_res$theta)
ggplot(dat, aes(x, y)) +
geom_point(size=2) +
geom_line(aes(y = yhat_gn), color="steelblue", linewidth=1) +
labs(title = "Hasil Estimasi Gauss–Newton",
subtitle = paste0("beta0 ≈ ", round(gn_res$theta[1],3),
", beta1 ≈ ", round(gn_res$theta[2],3)),
x = "X", y = "Y") +
theme_minimal()
Levenberg–Marquardt: Estimasi Manual (Damped GN) Satu Langkah LM
Kita gunakan \(D =
\operatorname{diag}(J^\top J)\).
Langkah \(\Delta\) diperoleh dari:
\[ (J^\top J + \lambda D)\,\Delta = J^\top r \]
lm_trial <- function(x, y, theta, lambda){
r <- resid_vec(x, y, theta)
J <- J_exp(x, theta)
JTJ <- crossprod(J,J)
D <- diag(diag(JTJ))
rhs <- crossprod(J, r)
lhs <- JTJ + lambda * D
Delta <- solve(lhs, rhs)
list(Delta = as.numeric(Delta),
theta_new = theta + as.numeric(Delta),
SSE_old = sum(r^2),
SSE_new = sum(resid_vec(x,y,theta + as.numeric(Delta))^2))
}
Prosedur Adaptif \(\lambda\) Jika \(S(\theta_{\text{baru}}) <
S(\theta_{\text{lama}})\), kita kecilkan \(\lambda\);
jika tidak, besarkan \(\lambda\) dan
coba lagi.
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(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 <- 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 <- add_row(hist, iter=k, beta0=th[1], beta1=th[2],
lambda=lam, SSE=tried$SSE_old,
step_norm=NA_real_, accepted = FALSE)
# do not update theta; try again with bigger lambda
}
}
list(theta = th, history = hist)
}
lm_res <- lm_fit(x, y, theta_init = c(1.0, 0.0), lambda_init = 1e-2, nu=10)
lm_res$history
## # A tibble: 11 × 7
## iter beta0 beta1 lambda SSE step_norm accepted
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <lgl>
## 1 0 1 0 0.01 137. NA TRUE
## 2 1 1 0 0.1 137. NA FALSE
## 3 2 1 0 1 137. NA FALSE
## 4 3 1 0 10 137. NA FALSE
## 5 4 1.33 0.134 1 85.7 3.58e-1 TRUE
## 6 5 2.10 0.328 0.1 3.11 7.88e-1 TRUE
## 7 6 1.99 0.316 0.01 0.484 1.10e-1 TRUE
## 8 7 1.98 0.315 0.001 0.465 1.09e-2 TRUE
## 9 8 1.98 0.315 0.0001 0.465 2.59e-4 TRUE
## 10 9 1.98 0.315 0.00001 0.465 3.85e-6 TRUE
## 11 10 1.98 0.315 0.000001 0.465 8.79e-9 TRUE
lm_res$theta
## [1] 1.9753264 0.3148085
Kurva Hasil LM
dat$yhat_lm <- f_exp(dat$x, lm_res$theta)
ggplot(dat, aes(x, y)) +
geom_point(size=2) +
geom_line(aes(y = yhat_lm), color="firebrick", linewidth=1) +
labs(title = "Hasil Estimasi Levenberg–Marquardt",
subtitle = paste0("beta0 ≈ ", round(lm_res$theta[1],3),
", beta1 ≈ ", round(lm_res$theta[2],3)),
x = "X", y = "Y") +
theme_minimal()
Pembanding: nls() (Built-in R)
fit_nls <- nls(y ~ b0*exp(b1*x), data = dat,
start = list(b0 = 1, b1 = 0))
summary(fit_nls)
##
## Formula: y ~ b0 * exp(b1 * x)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b0 1.97533 0.16309 12.11 0.000267 ***
## b1 0.31481 0.01969 15.99 8.95e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3411 on 4 degrees of freedom
##
## Number of iterations to convergence: 5
## Achieved convergence tolerance: 2.846e-07
coef(fit_nls)
## b0 b1
## 1.9753263 0.3148085
Inferensi (SE & Interval Kepercayaan)
Pada solusi \(\hat{\theta}\), aproksimasi kovarians diberikan oleh:
\[ \hat{\sigma}^2 \;=\; \frac{S(\hat{\theta})}{n - p}, \]
dan
\[ \widehat{\operatorname{Var}}(\hat{\theta}) \;\approx\; \hat{\sigma}^2 \,\Big( J(\hat{\theta})^{\top} J(\hat{\theta}) \Big)^{-1}. \]
dengan:
theta_hat <- lm_res$theta # pakai hasil LM (boleh ganti GN/nls)
r_hat <- resid_vec(x, y, theta_hat)
SSE <- sum(r_hat^2)
n <- length(x); p <- 2
sigma2_hat <- SSE/(n-p)
J_hat <- J_exp(x, theta_hat)
cov_theta <- sigma2_hat * solve(crossprod(J_hat, J_hat))
se_theta <- sqrt(diag(cov_theta))
theta_hat; se_theta
## [1] 1.9753264 0.3148085
## d_beta0 d_beta1
## 0.16309142 0.01969227
Interval kepercayaan Wald 95% untuk tiap parameter \(\theta_j\) adalah:
\[ \hat{\theta}_j \;\pm\; 1.96 \; SE(\hat{\theta}_j). \]
dengan: - \(\hat{\theta}_j\): estimasi parameter ke-\(j\), - \(SE(\hat{\theta}_j)\): standar error dari estimasi parameter ke-\(j\).
ci <- cbind(
est = theta_hat,
LCL = theta_hat - 1.96*se_theta,
UCL = theta_hat + 1.96*se_theta
)
rownames(ci) <- c("beta0","beta1")
ci
## est LCL UCL
## beta0 1.9753264 1.6556672 2.2949856
## beta1 0.3148085 0.2762116 0.3534053
Soal
Ulangi GN & LM dengan tebakan awal berbeda-beda (misalnya \((0.5,-0.5)\), \((5,0.5)\)).
Catat konvergensi, jumlah iterasi, dan SSE akhir.
Generasikan data dengan error multiplikatif:
\[ Y = \beta_0 e^{\beta_1 X} \cdot \varepsilon, \qquad \ln \varepsilon \sim N(0,\tau^2) \]
Estimasi OLS pada \(\ln Y\) dan bandingkan hasilnya dengan NLS.
Untuk hasil terbaik (GN/LM), plot residu vs fitted
dan QQ-plot residu.
Apakah asumsi normalitas dan homoskedastisitas terlihat wajar?
Terapkan GN/LM untuk model Michaelis–Menten dan
logistik.
Laporkan \(\hat{\theta}\), SE, CI, dan
interpretasi substantif parameternya.
Jawaban
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) } # theta = c(beta0, beta1)
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)
}
# - Gauss-Newton yang stabil (QR + fallback ridge + backtracking) -
gn_fit_safe <- function(x, y, theta_init, max_it = 50, tol = 1e-8,
ridge = 1e-10, max_bt = 8) {
th <- theta_init
sse_prev <- sum(resid_vec(x, y, th)^2)
for (k in 1:max_it) {
r <- resid_vec(x, y, th)
J <- J_exp(x, th)
# Coba solusi LS dengan QR (lebih stabil)
Delta <- tryCatch(
qr.coef(qr(J), r),
error = function(e) rep(NA_real_, ncol(J))
)
# Jika gagal/NA, pakai ridge kecil pada normal equations
if (any(!is.finite(Delta))) {
JTJ <- crossprod(J, J)
rhs <- crossprod(J, r)
Delta <- tryCatch(
solve(JTJ + ridge * diag(ncol(J)), rhs),
error = function(e) rep(0, ncol(J))
)
}
# Backtracking line search: kurangi langkah jika SSE naik
step <- 1.0
improved <- FALSE
for (bt in 0:max_bt) {
th_new <- th + as.numeric(step * Delta)
sse_new <- sum(resid_vec(x, y, th_new)^2)
if (is.finite(sse_new) && (sse_new <= sse_prev)) {
th <- th_new
sse_prev <- sse_new
improved <- TRUE
break
}
step <- step / 2
}
# Hentikan jika tidak ada perbaikan atau langkah sangat kecil
if (!improved || sqrt(sum((step * Delta)^2)) < tol) break
}
th
}
# - Coba beberapa starting values (hindari beta0=0) -
starts <- list(c(0.5, -0.5), c(5, 0.5), c(1, 0.0))
results <- lapply(starts, function(st) gn_fit_safe(x, y, st))
names(results) <- c("start (0.5,-0.5)", "start (5,0.5)", "start (1,0)")
results
## $`start (0.5,-0.5)`
## [1] 1.9753264 0.3148085
##
## $`start (5,0.5)`
## [1] 1.9753264 0.3148085
##
## $`start (1,0)`
## [1] 1.9753264 0.3148085
Kesimpulan
Hasil akhirnya relatif konsisten (mendekati \(\beta_0 \approx 2\), \(\beta_1 \approx 0.3\))
meskipun starting value berbeda. Namun, jumlah iterasi
dan kecepatan konvergensi bisa berbeda.
Untuk data dengan error multiplikatif, model bisa dilog-kan agar menjadi linear.
set.seed(456)
y_mult <- beta0_true * exp(beta1_true*x) * exp(rnorm(length(x), sd=0.2))
df <- data.frame(x,y_mult)
# Estimasi OLS pada log(y)
fit_lm <- lm(log(y_mult) ~ x, data=df)
summary(fit_lm)
##
## Call:
## lm(formula = log(y_mult) ~ x, data = df)
##
## Residuals:
## 1 2 3 4 5 6
## -0.20616 0.19319 0.23530 -0.19636 -0.05516 0.02919
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.63060 0.15252 4.135 0.01444 *
## x 0.29371 0.05038 5.830 0.00431 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2107 on 4 degrees of freedom
## Multiple R-squared: 0.8947, Adjusted R-squared: 0.8684
## F-statistic: 33.99 on 1 and 4 DF, p-value: 0.004312
# konversi kembali
beta0_hat <- exp(coef(fit_lm)[1])
beta1_hat <- coef(fit_lm)[2]
c(beta0_hat, beta1_hat)
## (Intercept) x
## 1.8787395 0.2937094
Kesimpulan: Dengan error multiplikatif, model linear log cocok → estimasi lebih sederhana (OLS), hasilnya dekat dengan nilai parameter sebenarnya.
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.
Gunakan nls() dengan data simulasi Michaelis–Menten.
set.seed(789)
Vmax_true <- 2; Km_true <- 1.2
X <- seq(0,6,length.out=40)
Y <- (Vmax_true*X)/(Km_true+X) + rnorm(length(X),sd=0.05)
df_mm <- data.frame(X,Y)
fit_mm <- nls(Y ~ (Vmax*X)/(Km+X), data=df_mm,
start=list(Vmax=1.5,Km=1))
summary(fit_mm)
##
## Formula: Y ~ (Vmax * X)/(Km + X)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Vmax 1.9981 0.0234 85.40 <2e-16 ***
## Km 1.2415 0.0491 25.28 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03623 on 38 degrees of freedom
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 1.317e-06
confint(fit_mm)
## 2.5% 97.5%
## Vmax 1.952336 2.04651
## Km 1.146484 1.34395
Interpretasi
Ringkasan Jawaban
Setelah parameter model nonlinear \(\hat{\theta}\) diperoleh (misal dengan
nonlinear least squares), kita ingin melakukan
inferensi: 1. Apakah parameter signifikan (≠ 0)?
2. Bagaimana ketidakpastian estimasi (interval kepercayaan)?
3. Seberapa baik model cocok dengan data (goodness-of-fit)?
Hipotesis umum untuk tiap parameter \(\theta_j\):
\[ H_0 : \theta_j = 0 \quad \text{vs} \quad H_1 : \theta_j \neq 0. \]
Statistik uji Wald:
\[ t_j = \frac{\hat{\theta}_j}{SE(\hat{\theta}_j)}, \]
dengan \(SE(\hat{\theta}_j)\) diperoleh dari akar diagonal matriks kovarians \(\widehat{\mathrm{Var}}(\hat{\theta})\).
Jika \(|t_j|\) cukup besar, atau \(p\text{-value}\) < 0.05, maka \(\theta_j\) signifikan.
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.
\(R^2\) (Pseudo-\(R^2\) pada regresi nonlinear) \[
R^2 = 1 - \frac{\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 \(y\) yang dijelaskan model. Nilai mendekati 1 → fit baik.
Akaike Information Criterion (AIC)
\[
AIC = n \cdot \ln\!\left(\frac{\text{SSE}}{n}\right) + 2p,
\] dengan:
- \(n\) = jumlah observasi,
- \(p\) = jumlah parameter.
Interpretasi: lebih kecil → model lebih baik (dengan penalti kompleksitas).
Kita gunakan model: \[ y = \beta_0 e^{\beta_1 x} + \varepsilon, \qquad \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
# Ringkasan dengan uji t dan p-value
tidy(fit, conf.int = TRUE)
## # A tibble: 2 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 b0 2.01 0.0968 20.7 6.66e- 9 1.79 2.23
## 2 b1 0.300 0.00540 55.6 9.93e-13 0.288 0.313
Interpretasi: jika p-value < 0.05 → parameter signifikan berbeda dari 0.
confint(fit)
## 2.5% 97.5%
## b0 1.7917267 2.2341574
## b1 0.2879737 0.3126612
Interval kepercayaan memberi rentang nilai parameter yang konsisten dengan data.
n <- nrow(dat)
p <- length(coef(fit))
yhat <- fitted(fit)
SSE <- sum(resid(fit)^2)
SST <- sum((y - mean(y))^2)
R2 <- 1 - SSE/SST
AIC <- n*log(SSE/n) + 2*p
list(SSE=SSE, R2=R2, AIC=AIC)
## $SSE
## [1] 2.387234
##
## $R2
## [1] 0.9984577
##
## $AIC
## [1] -12.80536
Laporkan dalam format R Markdown lengkap dengan narasi, rumus, kode, output, dan visualisasi.
Model regresi eksponensial dengan error aditif: \[ Y_i \;=\; f_i(\boldsymbol\theta) + \varepsilon_i \;=\; \beta_0 e^{\beta_1 x_i} + \varepsilon_i, \qquad \varepsilon_i \overset{iid}{\sim} \mathcal N(0,\sigma^2), \quad i=1,\dots,n. \]
Vektor parameter: \[ \boldsymbol\theta \;=\; (\beta_0,\beta_1)^\top,\qquad p=2. \]
Catatan: Di bawah asumsi normal, NLS = MLE untuk \(\boldsymbol\theta\) dan \(\sigma^2\).
Residual dan fungsi objektif (jumlah kuadrat residual): \[ r_i(\boldsymbol\theta) \;=\; y_i - \beta_0 e^{\beta_1 x_i}, \qquad S(\boldsymbol\theta) \;=\; \sum_{i=1}^n \big[ r_i(\boldsymbol\theta) \big]^2 \;=\; \lVert \mathbf r(\boldsymbol\theta)\rVert_2^2. \]
Jacobian \(J(\boldsymbol\theta)\in\mathbb R^{n\times p}\) dengan baris ke-\(i\): \[ J_{i1} \;=\; \frac{\partial f_i}{\partial \beta_0} \;=\; e^{\beta_1 x_i}, \qquad J_{i2} \;=\; \frac{\partial f_i}{\partial \beta_1} \;=\; \beta_0 x_i e^{\beta_1 x_i}. \]
Gradien di \(\boldsymbol\theta\): \[ \nabla S(\boldsymbol\theta) \;=\; \mathbf g(\boldsymbol\theta) \;=\; -2\, J(\boldsymbol\theta)^\top \mathbf r(\boldsymbol\theta). \]
Aproksimasi Hessian Gauss–Newton: \[ H_{\text{GN}}(\boldsymbol\theta) \;\approx\; 2\, J(\boldsymbol\theta)^\top J(\boldsymbol\theta). \]
(Hessian eksak: \(H = 2J^\top J - 2\sum_{i=1}^n r_i\,H_i\), dengan \(H_i\) Hessian \(f_i\).)
Pada iterasi ke-\(k\), langkah perbaikan \(\Delta^{(k)}\) diperoleh dari: \[ \big(J^{(k)\top} J^{(k)}\big)\,\Delta^{(k)} \;=\; J^{(k)\top}\,\mathbf r^{(k)}, \] dengan \[ \mathbf r^{(k)} \;=\; \big(y_i - \beta_0^{(k)} e^{\beta_1^{(k)} x_i}\big)_{i=1}^n, \qquad J^{(k)} \;=\; J(\boldsymbol\theta^{(k)}). \]
Pembaruan parameter: \[ \boldsymbol\theta^{(k+1)} \;=\; \boldsymbol\theta^{(k)} + \Delta^{(k)}. \]
Kriteria henti (contoh): \[ \lVert \Delta^{(k)} \rVert_\infty < \varepsilon_\theta \quad \text{atau}\quad \big|S^{(k+1)} - S^{(k)}\big| < \varepsilon_S \quad \text{atau}\quad \lVert \mathbf g^{(k)} \rVert_\infty < \varepsilon_g. \]
Alternatif redaman (Levenberg–Marquardt): \[ \big(J^{(k)\top}J^{(k)} + \lambda^{(k)} D^{(k)}\big)\Delta^{(k)} \;=\; J^{(k)\top} \mathbf r^{(k)}, \quad D^{(k)}=I \ \text{atau}\ \operatorname{diag}(J^{(k)\top}J^{(k)}). \]
Setelah konvergen di \(\hat{\boldsymbol\theta}\), taksiran varians error: \[ \hat\sigma^2 \;=\; \frac{S(\hat{\boldsymbol\theta})}{n-p}. \]
Kovarians asimtotik parameter (Wald): \[ \widehat{\operatorname{Cov}}\!\big(\hat{\boldsymbol\theta}\big) \;\approx\; \hat\sigma^2 \,\big(J(\hat{\boldsymbol\theta})^\top J(\hat{\boldsymbol\theta})\big)^{-1}. \]
Galat baku: \[ \operatorname{SE}(\hat\theta_j) \;=\; \sqrt{\,\hat\sigma^2\,\big[(J^\top J)^{-1}\big]_{jj}\,}\ \Big|_{\boldsymbol\theta=\hat{\boldsymbol\theta}}. \]
(Opsional, robust sandwich—untuk heteroskedastisitas): \[ \widehat{\operatorname{Cov}}_{\text{robust}}\!\big(\hat{\boldsymbol\theta}\big) \;=\; \big(J^\top J\big)^{-1} J^\top \operatorname{diag}\!\big(r_i(\hat{\boldsymbol\theta})^2\big)\, J\,\big(J^\top J\big)^{-1}. \]
Asimtotik normalitas (Wald): \[ \hat{\boldsymbol\theta} \;\dot{\sim}\; \mathcal N\!\Big(\boldsymbol\theta,\ \sigma^2 \big(J^\top J\big)^{-1}\Big). \]
Statistik-\(z\) satu parameter untuk hipotesis \(H_0:\ \theta_j=\theta_{j,0}\): \[ z_j \;=\; \frac{\hat\theta_j - \theta_{j,0}}{\operatorname{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}\,\operatorname{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\boldsymbol\theta=r\)): \[ W \;=\; \big(R\hat{\boldsymbol\theta}-r\big)^\top \Big[ R\,\widehat{\operatorname{Cov}}(\hat{\boldsymbol\theta})\,R^\top \Big]^{-1} \big(R\hat{\boldsymbol\theta}-r\big) \;\;\dot{\sim}\;\; \chi^2_q. \]
Untuk titik baru \(x_\star\), prediksi mean: \[ \hat y_\star \;=\; f(x_\star,\hat{\boldsymbol\theta}) \;=\; \hat\beta_0\,e^{\hat\beta_1 x_\star}. \]
Linearization delta method: turunan terhadap parameter di \(\hat{\boldsymbol\theta}\), \[ \mathbf g_\star \;=\; \begin{bmatrix} \frac{\partial f(x_\star,\boldsymbol\theta)}{\partial \beta_0} \\ \frac{\partial f(x_\star,\boldsymbol\theta)}{\partial \beta_1} \end{bmatrix}_{\boldsymbol\theta=\hat{\boldsymbol\theta}} = \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}\big(\hat y_\star\big) \;\approx\; \mathbf g_\star^\top\,\widehat{\operatorname{Cov}}(\hat{\boldsymbol\theta})\,\mathbf g_\star. \]
CI untuk mean \(f(x_\star,\boldsymbol\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\,}. \]
Catatan
Contoh GN Manual (n = 5) untuk Model Eksponensial
Data Kita gunakan 5 observasi:
\[ \begin{aligned} \mathbf x &= (0,\;1,\;2,\;3,\;4)^\top,\\ \mathbf y &= (2.100000,\;2.649718,\;3.724238,\;4.819206,\;6.790234)^\top . \end{aligned} \]
Model \[ y_i \;=\; \beta_0 e^{\beta_1 x_i} + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal N(0,\sigma^2). \]
Iterasi 0 (tebakan awal dari regresi \(\ln y\) pada \(x\)) \[ \boldsymbol\theta^{(0)}=\begin{bmatrix}\beta_0^{(0)}\\ \beta_1^{(0)}\end{bmatrix} =\begin{bmatrix}2.043817\\ 0.294525\end{bmatrix}. \]
Nilai fungsi, residual, dan Jacobian di \(\boldsymbol\theta^{(0)}\): \[ \mathbf f^{(0)}=\beta_0^{(0)} e^{\beta_1^{(0)}\mathbf x} =\begin{bmatrix} 2.043817\\ 2.743801\\ 3.683523\\ 4.945088\\ 6.638725 \end{bmatrix},\qquad \mathbf r^{(0)}=\mathbf y-\mathbf 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}\mathbf 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}\mathbf r^{(0)} \;\;\Rightarrow\;\; \Delta^{(0)}= \begin{bmatrix} -0.041728\\ \;\;0.008108 \end{bmatrix}. \]
Pembaruan parameter dan SSE: \[ \boldsymbol\theta^{(1)}=\boldsymbol\theta^{(0)}+\Delta^{(0)} =\begin{bmatrix}2.002088\\ 0.302633\end{bmatrix},\qquad S(\boldsymbol\theta^{(0)})=0.052467,\; S(\boldsymbol\theta^{(1)})=0.042495. \]
Iterasi 1 Hitung ulang \(J^{(1)}\), \(\mathbf r^{(1)}\) pada \(\boldsymbol\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}, \quad \mathbf 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}\mathbf r^{(1)}= \begin{bmatrix} 0.007541\\ 0.061600 \end{bmatrix}. \]
Langkah GN berikutnya: \[ \Delta^{(1)}= \begin{bmatrix} -0.000446\\ \;\;0.000125 \end{bmatrix}, \qquad \boldsymbol\theta^{(2)}=\begin{bmatrix}2.001642\\ 0.302758\end{bmatrix}, \qquad S(\boldsymbol\theta^{(2)})=0.042491. \]
Terlihat \(S(\theta)\) menurun dan perubahan parameter makin kecil → konvergen.
Kita gunakan model regresi eksponensial dengan error aditif: \[ Y_i = f_i(\boldsymbol\theta) + \varepsilon_i = \beta_0 e^{\beta_1 x_i} + \varepsilon_i,\qquad \varepsilon_i \overset{iid}{\sim} \mathcal N(0,\sigma^2), \quad i=1,\dots,n. \]
Vektor parameter: \[ \boldsymbol\theta = (\beta_0,\beta_1)^\top, \qquad p=2. \]
LM menggabungkan Gauss–Newton dan Gradient
Descent.
Iterasi ke-\(k\) menyelesaikan sistem
modifikasi:
\[ \big(J^{(k)\top} J^{(k)} + \lambda^{(k)} D^{(k)}\big)\,\Delta^{(k)} = J^{(k)\top}\,\mathbf r^{(k)}, \]
dengan: - \(J^{(k)}\) : Jacobian di
\(\boldsymbol\theta^{(k)}\)
- \(\mathbf r^{(k)} = (y_i - \beta_0^{(k)}
e^{\beta_1^{(k)} 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.
Residual dan fungsi SSE sama dengan GN:
\[ r_i(\boldsymbol\theta) = y_i - \beta_0 e^{\beta_1 x_i}, \qquad S(\boldsymbol\theta)=\sum_{i=1}^n r_i(\boldsymbol\theta)^2. \]
Jacobian (baris ke-\(i\)): \[ J_{i1} = e^{\beta_1 x_i}, \qquad J_{i2} = \beta_0 x_i e^{\beta_1 x_i}. \]
Pembaruan parameter: \[ \boldsymbol\theta^{(k+1)} = \boldsymbol\theta^{(k)} + \Delta^{(k)}. \]
Aturan adaptasi \(\lambda\): - Jika
\(S(\boldsymbol\theta^{(k+1)}) <
S(\boldsymbol\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\)).
LM memodifikasi matriks normal: \[ H_{\text{LM}}^{(k)} = J^{(k)\top} J^{(k)} + \lambda^{(k)} D^{(k)}. \]
Gradien tetap sama: \[ \mathbf g^{(k)} = -2\,J^{(k)\top}\mathbf r^{(k)}. \]
Setelah konvergen di \(\hat{\boldsymbol\theta}\): \[ \hat\sigma^2 = \frac{S(\hat{\boldsymbol\theta})}{n-p}. \]
Kovarians asimtotik (Wald): \[ \widehat{\operatorname{Cov}}(\hat{\boldsymbol\theta}) \approx \hat\sigma^2\,(J^\top J)^{-1}\Big|_{\boldsymbol\theta=\hat{\boldsymbol\theta}}. \]
SE parameter: \[ \operatorname{SE}(\hat\theta_j) = \sqrt{\,\hat\sigma^2\,[(J^\top J)^{-1}]_{jj}\,}. \]
Distribusi asimtotik: \[ \hat{\boldsymbol\theta} \;\dot{\sim}\; \mathcal N\!\Big(\boldsymbol\theta,\ \sigma^2 (J^\top J)^{-1}\Big). \]
Statistik-\(z\): \[ z_j = \frac{\hat\theta_j - \theta_{j,0}}{\operatorname{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}\,\operatorname{SE}(\hat\theta_j). \]
Prediksi mean di \(x_\star\): \[ \hat y_\star = f(x_\star,\hat{\boldsymbol\theta}) = \hat\beta_0 e^{\hat\beta_1 x_\star}. \]
Gradien terhadap parameter: \[ \mathbf 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 \mathbf g_\star^\top \,\widehat{\operatorname{Cov}}(\hat{\boldsymbol\theta})\,\mathbf 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}. \]
Catatan
Bab ini menyajikan materi komprehensif tentang outlier dalam analisis regresi, mencakup:
Outlier adalah observasi yang menyimpang jauh dari pola umum data. Dalam regresi, pembedaan penting:
Secara umum, regresi OLS meminimalkan jumlah kuadrat residual \[ \hat{\beta} = \arg\min_{\beta}\sum_{i=1}^n (y_i - \mathbf{x}_i^\top \beta)^2, \] sehingga residual besar akan dikuadratkan dan diberi bobot tidak proporsional — ini yang membuat OLS rentan terhadap outlier.
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))
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\).
Kriteria praktis: \(|r_i|>2\) atau \(|t_i|>3\) untuk kandidat outlier.
Leverage \[ h_{ii}=\mathbf{x}_i^\top (X^\top X)^{-1}\mathbf{x}_i, \quad \text{ambang populer: } h_{ii}>\frac{2p}{n}, \] dengan \(p\) jumlah parameter (termasuk intercept), \(n\) ukuran sampel.
Cook’s Distance \[ D_i = \frac{r_i^2}{p}\cdot \frac{h_{ii}}{1-h_{ii}}. \] Aturan praktis: \(D_i>1\) perlu ditinjau.
DFFITS dan DFBETAS \[ \text{DFFITS}_i = \frac{\hat{y}_i - \hat{y}_{i(i)}}{\hat{\sigma}_{(i)}\sqrt{h_{ii}}},\quad \text{DFBETAS}_{ij} = \frac{\hat{\beta}_j - \hat{\beta}_{j(i)}}{\mathrm{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).
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))
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.
M‑estimasi meminimalkan \[ \hat{\beta} = \arg\min_\beta \sum_{i=1}^n \rho\!\left(\frac{y_i-\mathbf{x}_i^\top\beta}{\sigma}\right), \] dengan fungsi loss \(\rho(\cdot)\) tertentu dan fungsi influence \(\psi=\rho'\). Dua fungsi populer:
Huber loss dengan ambang \(c\): \[ \rho(e)= \begin{cases} \frac{1}{2}e^2, & |e|\le c\\ c|e|-\frac{1}{2}c^2, & |e|>c \end{cases} ,\quad \psi(e)= \begin{cases} e, & |e|\le c\\ c\,\mathrm{sign}(e), & |e|>c \end{cases} \]
Tukey’s bisquare (redescending): \[ \rho(e)= \begin{cases} \frac{c^2}{6}\left[1-\left(1-\left(\frac{e}{c}\right)^2\right)^3\right], & |e|\le c\\ \frac{c^2}{6}, & |e|>c \end{cases} ,\quad \psi(e)= \begin{cases} e\left(1-\left(\frac{e}{c}\right)^2\right)^2, & |e|\le c\\ 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)}\).
Nilai \(c\) adalah tuning constant yang menentukan kapan sebuah residual dianggap outlier.
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
Contoh
Efek outlier dibatasi pada \(\pm c\).
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 = \frac{\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\}, \quad c=1.345\)
Huber loss:
Observasi ke-3 (outlier) hanya dapat bobot 0.269, jauh lebih kecil dibanding observasi lain.
Iterasi IRLS
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 1
Output weights(fit_rob) akan menunjukkan:
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 \text{median}\big(|e_i - \text{median}(e)|\big)\]
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() (paket
MASS):
fit_rob$s # nilai sigma robust
## [1] 2.487055
summary(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 freedom
Catatan:
Meminimalkan \(\sum_i |y_i-\mathbf{x}_i^\top\beta|\) (regresi median). Lebih tahan terhadap outlier pada Y, namun kurang efisien jika error benar‑benar normal. Dapat diselesaikan dengan pemrograman linear atau algoritme khusus.
Menggeneralisasi LAD ke kuantil \(\tau\in(0,1)\) dengan meminimalkan check loss \[ \sum_{i=1}^n \rho_\tau(\varepsilon_i),\quad \rho_\tau(u)=u\cdot(\tau-\mathbf{1}\{u<0\}). \] Memungkinkan analisis efek heterogen di berbagai titik distribusi respons (median, P90, dst.).
library(MASS) # rlm (Huber & bisquare)
library(quantreg) # quantile regression
library(dplyr); library(broom)
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>
Diagnostik 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 95
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 OLS
Perbandingan 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")
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.5
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 × 5
## idx rstd leverage cooks flagged
## <int> <dbl> <dbl> <dbl> <lgl>
## 1 55 -5.82 0.0496 0.589 TRUE
## 2 95 -5.15 0.315 4.07 TRUE
## 3 80 -4.64 0.0656 0.504 TRUE
## 4 10 4.12 0.0151 0.0869 TRUE
## 5 58 -0.768 0.0999 0.0218 TRUE
## 6 68 0.725 0.0571 0.0106 TRUE
## 7 44 -0.157 0.0528 0.000456 TRUE
Tahapan Analisis
Template Pelaporan
Catatan Praktis
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.
Petunjuk
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 \(\psi(e)\) dalam robust regression dan bedanya pada OLS, Huber, dan Tukey’s bisquare.
Bagaimana memilih nilai ambang \(c\) pada Huber dan Tukey? Jelaskan peran \(\hat{\sigma}\) (skala robust).
Diberikan residual OLS: \(e=\{-0.5,\,1.2,\,5.0\}\) dan \(c=1.345\).
Diberikan Cook’s Distance untuk tiga observasi: \(\{0.2,\,0.8,\,1.5\}\).
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.937836
coef(fit_rob)
## (Intercept) x
## 0.9707001 1.9844546
C.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 TRUE
C.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 1
(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.
(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.
Outlier: observasi yang jauh dari pola mayoritas.
OLS meminimalkan \(\sum e_i^2\). Residual besar dikuadratkan → pengaruhnya meningkat sangat cepat (unbounded influence). Karena itu OLS sensitif terhadap outlier.
Influence function \(\psi(e)=\rho'(e)\) mengukur pengaruh residual pada estimasi.
Pemilihan \(c\):
B6. \(D_i =
\dfrac{r_i^2}{p}\cdot\dfrac{h_{ii}}{1-h_{ii}}\). Jika \(D_i>1\), observasi sangat
berpengaruh terhadap model.
B7. Contoh: Tukey’s bisquare, Least Absolute Deviation
(LAD), Quantile regression.
B8. Residual 5.0 adalah outlier; bobot
Huber \(=1.345/5=0.269 \ll 1\).
| Bagian | Item | Skor per Item | Jumlah Item | Skor Maks |
|---|---|---|---|---|
| A. Pilihan Ganda | A1–A5 | 2 | 5 | 10 |
| B. Isian Singkat | B6–B8 | 5 | 3 | 15 |
| Total | 25 |
Machine Learning (ML) telah berkembang pesat dalam berbagai bidang, termasuk pengenalan pola, pemodelan prediktif, dan pengambilan keputusan berbasis data. Model pembelajaran mesin yang kompleks, seperti regresi linier berganda, regresi logistik, serta model berbasis jaringan saraf tiruan, sering kali menghadapi permasalahan yang berkaitan dengan overfitting—yakni kondisi di mana model belajar terlalu spesifik terhadap data pelatihan sehingga kehilangan kemampuan generalisasi pada data baru (Bishop, 2006).
Salah satu pendekatan utama untuk mengatasi overfitting adalah regularisasi, sebuah teknik yang menambahkan penalti pada kompleksitas model untuk mencegah koefisien atau parameter model menjadi terlalu besar. Regularisasi sangat penting terutama dalam situasi di mana jumlah fitur lebih besar dari jumlah observasi (high-dimensional data) atau ketika terdapat multikolinearitas antar prediktor (Hastie, Tibshirani, & Friedman, 2009). Regularisasi memungkinkan model tetap sederhana dan menghindari eksploitasi kebisingan (noise) dalam data yang dapat mengarah pada kesalahan prediksi yang besar.
Overfitting dan Trade-off Bias-Variance Overfitting dalam ML terjadi ketika model sangat cocok terhadap data pelatihan, tetapi tidak mampu bekerja dengan baik pada data uji. Ini biasanya terjadi karena model terlalu kompleks dan menangkap pola acak dalam data yang seharusnya tidak digeneralisasi (Goodfellow, Bengio, & Courville, 2016). Dalam konteks bias-variance tradeoff, model dengan bias rendah dan varians tinggi cenderung mengalami overfitting, sementara model dengan bias tinggi dan varians rendah mengalami underfitting. Regularisasi bertujuan untuk mencapai keseimbangan antara bias dan varians dengan mengontrol kompleksitas model.
Peran Regularisasi dalam Mengatasi Kompleksitas Model Dalam konteks ML, generalization performance merupakan aspek krusial karena model harus mampu membuat prediksi akurat pada data yang belum pernah dilihat sebelumnya. Regularisasi memastikan bahwa model tidak hanya cocok terhadap data pelatihan tetapi juga dapat bekerja dengan baik pada data uji atau data dunia nyata. Hal ini terutama berguna dalam high-dimensional data, seperti yang ditemukan dalam bioinformatika (Ng, 2004), keuangan (Friedman et al., 2010), dan pengenalan gambar (Krizhevsky, Sutskever, & Hinton, 2012).
Gambaran Umum Metode Regularisasi Terdapat tiga metode regularasisasi yang umumnya digunakan yaitu Ridge Regression, Lasso Regression, dan Elastic Net.
Ridge Regression menambahkan penalti pada jumlah kuadrat dari koefisien regresi, yang secara matematis dinyatakan sebagai berikut:
\[\begin{equation} \min_{\beta} \sum_{i=1}^{n} (y_i - X_i\beta)^2 + \lambda \sum_{j=1}^{p} \beta_j^2 \label{eq:ridge} \end{equation}\]
Dimana \(\lambda\) adalah parameter regularisasi yang mengontrol tingkat penalti. Ridge Regression membantu menangani multikolinearitas dan memastikan koefisien tetap kecil tetapi tidak menjadi nol (Hoerl & Kennard, 1970). Lihat Persamaan \(\eqref{eq:ridge}\) untuk detail formulasi.
Lasso (Least Absolute Shrinkage and Selection Operator) menggunakan penalti berbasis norma L1:
\[\begin{equation} \min_{\beta} \sum_{i=1}^{n} (y_i - X_i\beta)^2 + \lambda \sum_{j=1}^{p} |\beta_j| \label{eq:lasso} \end{equation}\]
Regularisasi L1 menyebabkan beberapa koefisien regresi menjadi nol, sehingga Lasso dapat digunakan untuk seleksi fitur (Tibshirani, 1996). Lihat Persamaan \(\eqref{eq:lasso}\) untuk lebih jelasnya.
Elastic Net menggabungkan penalti L1 dan L2 untuk mengatasi kelemahan masing-masing metode. Regularisasi ini sangat bermanfaat dalam situasi di mana jumlah fitur lebih banyak daripada jumlah observasi (Zou & Hastie, 2005).
\[\begin{equation} \min_{\beta} \sum_{i=1}^{n} (y_i - X_i\beta)^2 + \lambda_1 \sum_{j=1}^{p} |\beta_j| + \lambda_2 \sum_{j=1}^{p} \beta_j^2 \label{eq:elasticnet} \end{equation}\]
Elastic Net sangat berguna ketika terdapat korelasi tinggi antara fitur-fitur dalam dataset. Persamaan \(\eqref{eq:elasticnet}\) menunjukkan formulasi matematisnya.
Ridge Regression adalah metode regresi yang menambahkan penalti L2 terhadap koefisien regresi (Hoerl & Kennard, 1970). Model regresi biasa didefinisikan sebagai:
\[\begin{equation} Y = X\beta + \epsilon \label{eq:regresi} \end{equation}\]
dengan \(X\) adalah matriks prediktor, \(\beta\) adalah vektor koefisien, dan \(\epsilon\) adalah error.
Ilustrasi Ridge Regression
Tujuan dari ilustrasi ini adalah untuk menunjukkan bagaimana regularisasi Ridge Regression bekerja dalam mengontrol nilai koefisien regresi dengan menambahkan penalti terhadap besar koefisien. Dengan membandingkan Sum Squared Residual (SSR) dan SSR Ridge, kita dapat melihat pengaruh nilai lambda terhadap model.
Misalkan kita punya pengamatan berikut:
library(knitr)
library(dplyr)
library(ggplot2)
X <- c(2, 3, 4, 6)
Y <- c(4, 3.5, 4.5, 5)
beta_0_hat <- 3.07
beta_1_values <- seq(-0.2, 0.6, by = 0.1)
lambda_values <- c(0, 30, 40, 50, 60, 70)
Menghitung Sum Squared Residual (SSR) untuk setiap \(\beta_1\)
ssr_values <- sapply(beta_1_values, function(beta_1) {
y_hat <- beta_0_hat + beta_1 * X # Prediksi model
sum((Y - y_hat) ^ 2) # SSR
})
ssr_df <- data.frame(beta_1 = beta_1_values, SSR = ssr_values)
kable(
ssr_df,
col.names = c("$\\beta_1$", "SSR"),
caption = "Tabel Sum Squared Residual (SSR) untuk setiap $\\beta_1$",
escape = FALSE
)
| \(\beta_1\) | SSR |
|---|---|
| -0.2 | 17.5996 |
| -0.1 | 11.5596 |
| 0.0 | 6.8196 |
| 0.1 | 3.3796 |
| 0.2 | 1.2396 |
| 0.3 | 0.3996 |
| 0.4 | 0.8596 |
| 0.5 | 2.6196 |
| 0.6 | 5.6796 |
Menghitung SSR Ridge
ridge_results <- expand.grid(beta_1 = beta_1_values, lambda = lambda_values) %>%
mutate(SSR = rep(ssr_values, times = length(lambda_values)),
SSR_Ridge = SSR + lambda * beta_1^2)
kable(
ridge_results,
col.names = c("$\\beta_1$", "$\\lambda$", "SSR", "SSR Ridge"),
caption = "Tabel SSR Ridge untuk berbagai nilai $\\lambda$",
escape = FALSE
)
| \(\beta_1\) | \(\lambda\) | SSR | SSR Ridge |
|---|---|---|---|
| -0.2 | 0 | 17.5996 | 17.5996 |
| -0.1 | 0 | 11.5596 | 11.5596 |
| 0.0 | 0 | 6.8196 | 6.8196 |
| 0.1 | 0 | 3.3796 | 3.3796 |
| 0.2 | 0 | 1.2396 | 1.2396 |
| 0.3 | 0 | 0.3996 | 0.3996 |
| 0.4 | 0 | 0.8596 | 0.8596 |
| 0.5 | 0 | 2.6196 | 2.6196 |
| 0.6 | 0 | 5.6796 | 5.6796 |
| -0.2 | 30 | 17.5996 | 18.7996 |
| -0.1 | 30 | 11.5596 | 11.8596 |
| 0.0 | 30 | 6.8196 | 6.8196 |
| 0.1 | 30 | 3.3796 | 3.6796 |
| 0.2 | 30 | 1.2396 | 2.4396 |
| 0.3 | 30 | 0.3996 | 3.0996 |
| 0.4 | 30 | 0.8596 | 5.6596 |
| 0.5 | 30 | 2.6196 | 10.1196 |
| 0.6 | 30 | 5.6796 | 16.4796 |
| -0.2 | 40 | 17.5996 | 19.1996 |
| -0.1 | 40 | 11.5596 | 11.9596 |
| 0.0 | 40 | 6.8196 | 6.8196 |
| 0.1 | 40 | 3.3796 | 3.7796 |
| 0.2 | 40 | 1.2396 | 2.8396 |
| 0.3 | 40 | 0.3996 | 3.9996 |
| 0.4 | 40 | 0.8596 | 7.2596 |
| 0.5 | 40 | 2.6196 | 12.6196 |
| 0.6 | 40 | 5.6796 | 20.0796 |
| -0.2 | 50 | 17.5996 | 19.5996 |
| -0.1 | 50 | 11.5596 | 12.0596 |
| 0.0 | 50 | 6.8196 | 6.8196 |
| 0.1 | 50 | 3.3796 | 3.8796 |
| 0.2 | 50 | 1.2396 | 3.2396 |
| 0.3 | 50 | 0.3996 | 4.8996 |
| 0.4 | 50 | 0.8596 | 8.8596 |
| 0.5 | 50 | 2.6196 | 15.1196 |
| 0.6 | 50 | 5.6796 | 23.6796 |
| -0.2 | 60 | 17.5996 | 19.9996 |
| -0.1 | 60 | 11.5596 | 12.1596 |
| 0.0 | 60 | 6.8196 | 6.8196 |
| 0.1 | 60 | 3.3796 | 3.9796 |
| 0.2 | 60 | 1.2396 | 3.6396 |
| 0.3 | 60 | 0.3996 | 5.7996 |
| 0.4 | 60 | 0.8596 | 10.4596 |
| 0.5 | 60 | 2.6196 | 17.6196 |
| 0.6 | 60 | 5.6796 | 27.2796 |
| -0.2 | 70 | 17.5996 | 20.3996 |
| -0.1 | 70 | 11.5596 | 12.2596 |
| 0.0 | 70 | 6.8196 | 6.8196 |
| 0.1 | 70 | 3.3796 | 4.0796 |
| 0.2 | 70 | 1.2396 | 4.0396 |
| 0.3 | 70 | 0.3996 | 6.6996 |
| 0.4 | 70 | 0.8596 | 12.0596 |
| 0.5 | 70 | 2.6196 | 20.1196 |
| 0.6 | 70 | 5.6796 | 30.8796 |
Interpretasi Hasil
Plot Perubahan Garis Regresi untuk Setiap Nilai \(\beta_1\)
plot_data <- data.frame()
for (beta_1 in beta_1_values) {
y_pred <- beta_0_hat + beta_1 * X
temp_data <- data.frame(X, Y_pred = y_pred, beta_1 = as.factor(beta_1))
plot_data <- rbind(plot_data, temp_data)
}
ggplot(plot_data, aes(x = X, y = Y_pred, color = factor(beta_1))) +
geom_line(size=0.5) +
geom_point(data = data.frame(X, Y), aes(x = X, y = Y),
color = "black", size = 3) +
labs(title = expression("Perubahan Garis Regresi untuk Setiap"~beta[1]),
x = "X", y = "Y", color = expression(beta[1])) +
theme_bw()
Plot \(\beta_1\) vs SSR untuk berbagai \(\lambda\)
ggplot(ridge_results, aes(x = beta_1, y = SSR_Ridge, color = factor(lambda))) +
geom_line(size=0.5) +
geom_point() +
labs(title = expression(beta[1] ~ " vs SSR
Ridge untuk Berbagai Nilai " ~ lambda),
x = expression(beta[1]), y = "SSR", color = expression(lambda)) +
theme_bw()
Plot ini menunjukkan bahwa nilai optimal untuk \(\beta_1\) tidak pernah mencapai nilai nol untuk semua nilai \(\lambda\) .
Metode Estimasi Fungsi yang diminimalkan dalam Ridge Regression adalah:
\[\begin{equation} J(\beta) = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} \beta_j^2. \label{eq:ridge1} \end{equation}\]
Detail Metode Estimasi dalam Ridge Regression
1. Mendefinisikan Fungsi Objektif
Fungsi yang diminimalkan dalam Ridge Regression adalah:
\[\begin{equation} J(\beta) = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} \beta_j^2. \label{eq:ridge3} \end{equation}\]
Karena dalam regresi linier,
\[\begin{equation} \hat{y} = X \beta, \label{eq:ridge4} \end{equation}\]
maka fungsi objektif dapat ditulis dalam bentuk matriks sebagai:
\[\begin{equation} J(\beta) = (y - X\beta)^T (y - X\beta) + \lambda \beta^T \beta. \label{eq:ridge5} \end{equation}\]
2. Mencari Turunan Pertama \(\frac{\partial J}{\partial \beta}\)
Menggunakan aturan turunan dalam aljabar matriks:
\[\begin{equation} \frac{\partial}{\partial \beta} (y - X\beta)^T (y - X\beta) = -2X^T(y - X\beta). \label{eq:ridge6} \end{equation}\]
Untuk regularisasi:
\[\begin{equation} \frac{\partial}{\partial \beta} \lambda \beta^T \beta = 2\lambda \beta. \label{eq:ridge7} \end{equation}\]
Sehingga, turunan pertama dari \(J(\beta)\) adalah:
\[\begin{equation} \frac{\partial J}{\partial \beta} = -2X^T(y - X\beta) + 2\lambda \beta. \label{eq:ridge8} \end{equation}\]
3. Menyamakan Turunan dengan Nol
Menyusun ulang persamaan: \[\begin{equation} X^T X\beta - X^T y + \lambda \beta = 0. \label{eq:ridge9} \end{equation}\]
Maka:
\[\begin{equation} (X^T X + \lambda I)\beta = X^T y. \label{eq:ridge10} \end{equation}\]
4. Solusi Estimasi Parameter \(\hat{\beta}\)
Mengalikan kedua sisi dengan invers matriks \((X^T X + \lambda I)^{-1}\):
\[ \hat{\beta} = (X^{T} X + \lambda I)^{-1} X^{T} y \] {#eq:ridge10}
Ilustrasi Ridge Regression dalam Kehadiran Multicollinearity
Tujuan dari ilustrasi ini adalah untuk menunjukkan bagaimana Ridge Regression dapat menangani multicollinearity lebih baik dibandingkan OLS (Ordinary Least Squares) ketika terdapat korelasi tinggi antara variabel prediktor dalam dataset besar.
Simulasi Data dengan Multicollinearity
library(MASS)
library(glmnet)
library(ggplot2)
library(dplyr)
set.seed(123)
n <- 500 # Jumlah sampel
p <- 50 # Jumlah prediktor
# Membuat matriks kovarians dengan korelasi tinggi
Sigma <- matrix(0.9, nrow = p, ncol = p) + diag(0.1, p)
X <- mvrnorm(n, mu = rep(0, p), Sigma = Sigma) # Matriks prediktor
# Koefisien sebenarnya
beta_true <- runif(p, -2, 2)
y <- X %*% beta_true + rnorm(n, sd = 5) # Variabel respon
data <- data.frame(y, X)
Estimasi Model OLS dan Ridge Regression
# Model OLS
ols_model <- lm(y ~ ., data = data)
ols_coeff <- coef(ols_model)[-1] # Mengabaikan intercept
# Model Ridge Regression
x_matrix <- as.matrix(data[, -1])
y_vector <- data$y
ridge_model <- glmnet(x_matrix, y_vector, alpha = 0) # alpha = 0 untuk Ridge
cv_ridge <- cv.glmnet(x_matrix, y_vector, alpha = 0) # Cross-validation untuk lambda optimal
ridge_coeff <- coef(cv_ridge, s = "lambda.min")[-1] # Koefisien dengan lambda optimal
Visualisasi Model Ridge Regression
Plot Koefisien Ridge Regression terhadap \(\lambda\)
par(mfrow = c(1,1)) # Pastikan hanya ada satu plot dalam satu area
plot(ridge_model, xvar = "lambda", label = TRUE, cex.axis = 0.7, las = 2) # Menyesuaikan ukuran label
Interpretasi Plot Koefisien Ridge Regression terhadap \(\lambda\)
Shrinkage Koefisien saat lambda Meningkat
Multicollinearity Terkendali
lambda Besar vs. lambda Kecil
Kesimpulan
Cross-Validation Ridge Regression
plot(cv_ridge)
Interpretasi Plot Cross-Validation Ridge Regression
Makna dari Plot Plot ini merupakan hasil cross-validation untuk memilih nilai lambda optimal dalam Ridge Regression. Sumbu x menunjukkan log(lambda), sedangkan sumbu y menunjukkan Mean Squared Error (MSE) untuk setiap nilai lambda yang diuji.
Pola yang Terlihat dalam Plot
Bagian kiri grafik (lambda kecil):
MSE lebih rendah tetapi tidak stabil.
Ridge Regression mendekati OLS, yang berarti model masih memiliki variabilitas tinggi.
Bagian tengah grafik (nilai optimal lambda):
Dapat diamati nilai lambda optimal yang ditunjukkan oleh garis vertikal putus-putus.
Pada titik ini, MSE minimum tercapai, yang menunjukkan keseimbangan antara bias dan varians.
Bagian kanan grafik (lambda besar):
MSE meningkat secara signifikan.
Ridge Regression menyusutkan hampir semua koefisien menuju nol, menyebabkan underfitting.
Kesimpulan dari Plot
Plot Mean Squared Error untuk Berbagai \(\alpha\)
lambda_seq <- 10^seq(-5, 2, length.out = 100)
cv_fit_1 <- cv.glmnet(x_matrix, y_vector, alpha = 1, lambda = lambda_seq)
cv_fit_05 <- cv.glmnet(x_matrix, y_vector, alpha = 0.5, lambda = lambda_seq)
cv_fit_0 <- cv.glmnet(x_matrix, y_vector, alpha = 0, lambda = lambda_seq)
plot(log(cv_fit_1$lambda), cv_fit_1$cvm, col = "red", pch = 19,
xlab = "log(lambda)", ylab = "Mean-Squared Error")
points(log(cv_fit_05$lambda), cv_fit_05$cvm, col = "gray", pch = 19)
points(log(cv_fit_0$lambda), cv_fit_0$cvm, col = "blue", pch = 19)
legend("topright", legend = c("alpha=1", "alpha=.5", "alpha 0"),
col = c("red", "gray", "blue"), pch = 19)
Interpretasi Plot Cross-Validation Elastic Net Regression
Makna dari Plot Plot ini menunjukkan hasil cross-validation untuk memilih nilai lambda optimal dalam Elastic Net Regression. Tiga garis berbeda menunjukkan hasil untuk alpha = 1 (Lasso, merah), alpha = 0.5 (Elastic Net, abu-abu), dan alpha = 0 (Ridge, biru).
Sumbu x menunjukkan log(lambda), sedangkan sumbu y menunjukkan Mean Squared Error (MSE) untuk setiap nilai lambda yang diuji.
Pola yang Terlihat dalam Plot
Bagian kiri grafik (lambda kecil, log(lambda) negatif besar):
MSE rendah dan stabil untuk semua nilai alpha.
Model masih memiliki variabilitas tinggi karena regularisasi belum signifikan.
Bagian tengah grafik (lambda optimal):
MSE mencapai titik minimum untuk setiap nilai alpha.
Untuk alpha = 1 (Lasso, merah), MSE mulai meningkat lebih cepat dibandingkan Ridge.
Elastic Net (alpha = 0.5, abu-abu) menunjukkan keseimbangan antara Ridge dan Lasso.
Bagian kanan grafik (lambda besar, log(lambda) positif):
MSE meningkat tajam karena penalti yang berlebihan menyebabkan underfitting.
Lasso Regression (merah) paling cepat mengalami peningkatan MSE karena beberapa koefisien menjadi nol lebih awal.
Ridge Regression (biru) masih mempertahankan kestabilan lebih lama sebelum akhirnya mengalami peningkatan MSE.
Kesimpulan dari Plot
Perbandingan Koefisien OLS vs Ridge
comparison <- data.frame(
Predictor = paste0("X", 1:p),
OLS = ols_coeff,
Ridge = ridge_coeff
)
comparison_long <- comparison %>% tidyr::pivot_longer(cols = c("OLS", "Ridge"),
names_to = "Model", values_to = "Coefficient")
# Visualisasi perbandingan koefisien
ggplot(comparison_long, aes(x = Predictor, y = Coefficient, fill = Model)) +
geom_bar(stat = "identity", position = "dodge") +
theme_minimal() +
labs(title = "Perbandingan Koefisien OLS vs Ridge Regression",
x = "Predictor",
y = "Estimated Coefficient") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
Interpretasi Hasil
Kesimpulan Ketika terdapat multicollinearity, OLS menghasilkan koefisien yang tidak stabil. Sebaliknya, Ridge Regression memberikan solusi yang lebih robust dengan mengecilkan koefisien yang terlalu besar. Ini menunjukkan bagaimana regularisasi dapat meningkatkan stabilitas model dalam dataset dengan banyak prediktor yang berkorelasi tinggi.
Lasso Regression (Least Absolute Shrinkage and Selection Operator) adalah metode yang menggunakan penalti L1, yang memungkinkan beberapa koefisien regresi menjadi nol (Tibshirani, 1996):
\[\begin{equation} J(\beta) = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} |\beta_j|. \label{eq:Lasso1} \end{equation}\]
Lasso Regression (Least Absolute Shrinkage and Selection Operator) adalah metode regresi dengan regularisasi L1 yang memungkinkan beberapa koefisien regresi menjadi nol (Tibshirani, 1996). Fungsi objektif yang diminimalkan adalah:
\[\begin{equation} J(\beta) = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} |\beta_j|. \label{eq:Lasso2} \end{equation}\]
Karena dalam regresi linier,
\[\begin{equation} \hat{y} = X \beta, \label{eq:Lasso3} \end{equation}\]
maka fungsi objektif dapat ditulis ulang dalam bentuk matriks sebagai:
\[\begin{equation} J(\beta) = (y - X\beta)^T (y - X\beta) + \lambda \sum_{j=1}^{p} |\beta_j|. \label{eq:Lasso4} \end{equation}\]
Estimasi Parameter dengan Turunan Pertama
Berbeda dengan Ridge Regression, Lasso Regression menggunakan penalti L1, yang tidak memiliki turunan yang mudah dihitung karena keberadaan fungsi nilai absolut. Mari kita hitung gradiennya secara komponen:
1. Turunan Pertama \(\frac{\partial J}{\partial \beta}\)
Turunan dari bagian pertama:
\[\begin{equation} \frac{\partial}{\partial \beta} (y - X\beta)^T (y - X\beta) = -2X^T(y - X\beta). \label{eq:Lasso5} \end{equation}\]
Turunan dari penalti L1 berbeda untuk setiap elemen \(\beta_j\):
\[\begin{equation} \frac{\partial}{\partial \beta_j} \lambda |\beta_j| = \lambda \cdot \text{sign}(\beta_j), \label{eq:Lasso6} \end{equation}\]
dengan:
\[\begin{equation} \text{sign}(\beta_j) = \begin{cases} 1, & \text{jika } \beta_j > 0 \\ -1, & \text{jika } \beta_j < 0 \\ 0, & \text{jika } \beta_j = 0 \end{cases}. \label{eq:Lasso7} \end{equation}\]
Sehingga, turunan pertama dari \(J(\beta)\) adalah:
\[\begin{equation} \frac{\partial J}{\partial \beta} = -2X^T(y - X\beta) + \lambda \cdot \text{sign}(\beta). \label{eq:Lasso8} \end{equation}\]
2. Menyamakan dengan Nol untuk Solusi Optimal
Menyusun ulang:
\[ X^{T}(X\beta - y) + \lambda_{1}\,\partial \lVert \beta \rVert_{1} + \lambda_{2}\beta = 0 \] {#eq:Elastic9}
Karena fungsi nilai absolut tidak terdiferensiasi di \(\beta_j = 0\), solusi parameter harus ditemukan dengan metode numerik, seperti Coordinate Descent atau Least Angle Regression (LARS).
Algoritma Estimasi Parameter: Coordinate Descent
Algoritma Coordinate Descent sering digunakan untuk menghitung solusi Lasso. Langkah-langkahnya adalah:
\[\begin{equation} \beta_j^{(t+1)} = S\left( \frac{X_j^T (y - X\beta_{-j})}{X_j^T X_j}, \lambda \right), \end{equation}\]
di mana \(S(z, \lambda)\) adalah fungsi soft-thresholding:
\[\begin{equation} S(z, \lambda) = \text{sign}(z) \max(|z| - \lambda, 0). \end{equation}\]
Ilustrasi Lasso Regression Tujuan dari ilustrasi ini adalah untuk menunjukkan bagaimana regularisasi Lasso Regression bekerja dalam mengontrol nilai koefisien regresi dengan menambahkan penalti terhadap besar koefisien. Dengan membandingkan Sum Squared Residual (SSR) dan SSR Lasso, kita dapat melihat pengaruh nilai lambda terhadap model.
Misalkan kita punya pengamatan berikut:
X <- c(2, 3, 4, 6)
Y <- c(4, 3.5, 4.5, 5)
beta_0_hat <- 3.07
beta_1_values <- seq(-0.2, 0.6, by = 0.1)
lambda_values <- c(0, 30, 40, 50, 60, 70)
Menghitung Sum Squared Residual (SSR) untuk setiap \(\beta_1\)
ssr_values <- sapply(beta_1_values, function(beta_1) {
y_hat <- beta_0_hat + beta_1 * X
sum((Y - y_hat)^2)
})
ssr_df <- data.frame(beta_1 = beta_1_values, SSR = ssr_values)
kable(
ssr_df,
col.names = c("$\\beta_1$", "SSR"),
caption = "Tabel Sum Squared Residual (SSR) untuk setiap $\\beta_1$",
escape = FALSE
)
| \(\beta_1\) | SSR |
|---|---|
| -0.2 | 17.5996 |
| -0.1 | 11.5596 |
| 0.0 | 6.8196 |
| 0.1 | 3.3796 |
| 0.2 | 1.2396 |
| 0.3 | 0.3996 |
| 0.4 | 0.8596 |
| 0.5 | 2.6196 |
| 0.6 | 5.6796 |
Menghitung SSR Ridge
lasso_results <- expand.grid(beta1 = beta_1_values, lambda = lambda_values) %>%
dplyr::mutate(
SSR = rep(ssr_values, times = length(lambda_values)),
SSR_Lasso = SSR + lambda * beta1^2
)
knitr::kable(
lasso_results,
col.names = c("$\\beta_1$", "$\\lambda$", "SSR", "SSR Ridge"),
caption = "Tabel SSR lasso untuk berbagai nilai $\\lambda$",
escape = FALSE
)
| \(\beta_1\) | \(\lambda\) | SSR | SSR Ridge |
|---|---|---|---|
| -0.2 | 0 | 17.5996 | 17.5996 |
| -0.1 | 0 | 11.5596 | 11.5596 |
| 0.0 | 0 | 6.8196 | 6.8196 |
| 0.1 | 0 | 3.3796 | 3.3796 |
| 0.2 | 0 | 1.2396 | 1.2396 |
| 0.3 | 0 | 0.3996 | 0.3996 |
| 0.4 | 0 | 0.8596 | 0.8596 |
| 0.5 | 0 | 2.6196 | 2.6196 |
| 0.6 | 0 | 5.6796 | 5.6796 |
| -0.2 | 30 | 17.5996 | 18.7996 |
| -0.1 | 30 | 11.5596 | 11.8596 |
| 0.0 | 30 | 6.8196 | 6.8196 |
| 0.1 | 30 | 3.3796 | 3.6796 |
| 0.2 | 30 | 1.2396 | 2.4396 |
| 0.3 | 30 | 0.3996 | 3.0996 |
| 0.4 | 30 | 0.8596 | 5.6596 |
| 0.5 | 30 | 2.6196 | 10.1196 |
| 0.6 | 30 | 5.6796 | 16.4796 |
| -0.2 | 40 | 17.5996 | 19.1996 |
| -0.1 | 40 | 11.5596 | 11.9596 |
| 0.0 | 40 | 6.8196 | 6.8196 |
| 0.1 | 40 | 3.3796 | 3.7796 |
| 0.2 | 40 | 1.2396 | 2.8396 |
| 0.3 | 40 | 0.3996 | 3.9996 |
| 0.4 | 40 | 0.8596 | 7.2596 |
| 0.5 | 40 | 2.6196 | 12.6196 |
| 0.6 | 40 | 5.6796 | 20.0796 |
| -0.2 | 50 | 17.5996 | 19.5996 |
| -0.1 | 50 | 11.5596 | 12.0596 |
| 0.0 | 50 | 6.8196 | 6.8196 |
| 0.1 | 50 | 3.3796 | 3.8796 |
| 0.2 | 50 | 1.2396 | 3.2396 |
| 0.3 | 50 | 0.3996 | 4.8996 |
| 0.4 | 50 | 0.8596 | 8.8596 |
| 0.5 | 50 | 2.6196 | 15.1196 |
| 0.6 | 50 | 5.6796 | 23.6796 |
| -0.2 | 60 | 17.5996 | 19.9996 |
| -0.1 | 60 | 11.5596 | 12.1596 |
| 0.0 | 60 | 6.8196 | 6.8196 |
| 0.1 | 60 | 3.3796 | 3.9796 |
| 0.2 | 60 | 1.2396 | 3.6396 |
| 0.3 | 60 | 0.3996 | 5.7996 |
| 0.4 | 60 | 0.8596 | 10.4596 |
| 0.5 | 60 | 2.6196 | 17.6196 |
| 0.6 | 60 | 5.6796 | 27.2796 |
| -0.2 | 70 | 17.5996 | 20.3996 |
| -0.1 | 70 | 11.5596 | 12.2596 |
| 0.0 | 70 | 6.8196 | 6.8196 |
| 0.1 | 70 | 3.3796 | 4.0796 |
| 0.2 | 70 | 1.2396 | 4.0396 |
| 0.3 | 70 | 0.3996 | 6.6996 |
| 0.4 | 70 | 0.8596 | 12.0596 |
| 0.5 | 70 | 2.6196 | 20.1196 |
| 0.6 | 70 | 5.6796 | 30.8796 |
Interpretasi Hasil
Plot Perubahan Garis Regresi untuk Setiap Nilai \(\beta_1\)
plot_data <- data.frame()
for (beta_1 in beta_1_values) {
y_pred <- beta_0_hat + beta_1 * X
temp_data <- data.frame(X, Y_pred = y_pred, beta_1 = as.factor(beta_1))
plot_data <- rbind(plot_data, temp_data)
}
ggplot(plot_data, aes(x = X, y = Y_pred, color = factor(beta_1))) +
geom_line(size=0.5) +
geom_point(data = data.frame(X, Y), aes(x = X, y = Y),
color = "black", size = 3) +
labs(title = expression("Perubahan Garis Regresi untuk Setiap"~beta[1]),
x = "X", y = "Y", color = expression(beta[1])) +
theme_bw()
Plot \(\beta_1\) vs SSR untuk berbagai \(\lambda\)
ggplot(lasso_results, aes(x = beta_1, y = SSR_Lasso, color = factor(lambda))) +
geom_line(size=0.5) +
geom_point() +
labs(title = expression(beta[1] ~ " vs SSR Lasso untuk
Berbagai Nilai " ~ lambda),
x = expression(beta[1]), y = "SSR", color = expression(lambda)) +
theme_bw()
Plot ini menunjukkan bahwa nilai optimal untuk \(\beta_1\) dapat mencapai nol untuk beberapa nilai \(\lambda\).
Metode Estimasi Lasso tidak memiliki solusi analitik seperti Ridge. Estimasi parameter dilakukan menggunakan Coordinate Descent Algorithm.
Ilustrasi Lasso Regression dalam Kehadiran Multicollinearity
Tujuan dari ilustrasi ini adalah untuk menunjukkan bagaimana Lasso Regression dapat menangani multicollinearity lebih baik dibandingkan OLS (Ordinary Least Squares) ketika terdapat korelasi tinggi antara variabel prediktor dalam dataset besar.
Simulasi Data dengan Multicollinearity
library(MASS)
library(glmnet)
library(ggplot2)
library(dplyr)
set.seed(123)
n <- 500 # Jumlah sampel
p <- 50 # Jumlah prediktor
# Membuat matriks kovarians dengan korelasi tinggi
Sigma <- matrix(0.9, nrow = p, ncol = p) + diag(0.1, p)
X <- mvrnorm(n, mu = rep(0, p), Sigma = Sigma) # Matriks prediktor
# Koefisien sebenarnya
beta_true <- runif(p, -2, 2)
y <- X %*% beta_true + rnorm(n, sd = 5) # Variabel respon
data <- data.frame(y, X)
Estimasi Model OLS dan Lasso Regression
# Model OLS
ols_model <- lm(y ~ ., data = data)
ols_coeff <- coef(ols_model)[-1] # Mengabaikan intercept
# Model Lasso Regression
x_matrix <- as.matrix(data[, -1])
y_vector <- data$y
lasso_model <- glmnet(x_matrix, y_vector, alpha = 1) # alpha = 1 untuk Lasso
cv_lasso <- cv.glmnet(x_matrix, y_vector, alpha = 1) # Cross-validation untuk lambda optimal
lasso_coeff <- coef(cv_lasso, s = "lambda.min")[-1] # Koefisien dengan lambda optimal
Visualisasi Model Lasso Regression
Plot Koefisien Lasso Regression terhadap lambda
par(mfrow = c(1,1)) # Pastikan hanya ada satu plot dalam satu area
plot(lasso_model, xvar = "lambda", label = TRUE, cex.axis = 0.7, las = 2) # Menyesuaikan ukuran label
Interpretasi Plot Koefisien Lasso Regression terhadap lambda
Cross-Validation Lasso Regression
plot(cv_lasso)
Interpretasi Plot Cross-Validation Lasso Regression
Plot Mean Squared Error untuk Berbagai Alpha
lambda_seq <- 10^seq(-5, 2, length.out = 100)
cv_fit_1 <- cv.glmnet(x_matrix, y_vector, alpha = 1, lambda = lambda_seq)
cv_fit_05 <- cv.glmnet(x_matrix, y_vector, alpha = 0.5, lambda = lambda_seq)
cv_fit_0 <- cv.glmnet(x_matrix, y_vector, alpha = 0, lambda = lambda_seq)
plot(log(cv_fit_1$lambda), cv_fit_1$cvm, col = "red", pch = 19,
xlab = "log(lambda)", ylab = "Mean-Squared Error")
points(log(cv_fit_05$lambda), cv_fit_05$cvm, col = "gray", pch = 19)
points(log(cv_fit_0$lambda), cv_fit_0$cvm, col = "blue", pch = 19)
legend("topright", legend = c("alpha=1 (Lasso)", "alpha=.5 (Elastic Net)",
"alpha 0 (Ridge)"), col = c("red", "gray", "blue"), pch = 19)
Interpretasi Perbandingan Alpha dalam Cross-Validation
Perbandingan Koefisien OLS vs Lasso
comparison <- data.frame(
Predictor = paste0("X", 1:p),
OLS = ols_coeff,
Lasso = lasso_coeff
)
comparison_long <- comparison %>% tidyr::pivot_longer(cols = c("OLS", "Lasso"),
names_to = "Model", values_to = "Coefficient")
# Visualisasi perbandingan koefisien
ggplot(comparison_long, aes(x = Predictor, y = Coefficient, fill = Model)) +
geom_bar(stat = "identity", position = "dodge") +
theme_minimal() +
labs(title = "Perbandingan Koefisien OLS vs Lasso Regression",
x = "Predictor",
y = "Estimated Coefficient") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
Kesimpulan
Elastic Net menggabungkan penalti L1 dan L2 dengan parameter keseimbangan \(\alpha\) (Zou & Hastie, 2005):
\[\begin{equation} J(\beta) = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda_1 \sum_{j=1}^{p} |\beta_j| + \lambda_2 \sum_{j=1}^{p} \beta_j^2. \label{eq:Elastic1} \end{equation}\]
Elastic Net menggabungkan penalti L1 dan L2 dengan parameter keseimbangan \(\alpha\) (Zou & Hastie, 2005). Fungsi objektif yang diminimalkan adalah:
\[\begin{equation} J(\beta) = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda_1 \sum_{j=1}^{p} |\beta_j| + \lambda_2 \sum_{j=1}^{p} \beta_j^2. \label{eq:ElasticNet2} \end{equation}\]
Karena dalam regresi linier,
\[\begin{equation} \hat{y} = X \beta, \label{eq:Elastic3} \end{equation}\]
maka fungsi objektif dapat ditulis ulang dalam bentuk matriks sebagai:
\[\begin{equation} J(\beta) = (y - X\beta)^T (y - X\beta) + \lambda_1 \sum_{j=1}^{p} |\beta_j| + \lambda_2 \beta^T \beta. \label{eq:Elastic4} \end{equation}\]
Estimasi Parameter dengan Turunan Pertama
Elastic Net merupakan kombinasi dari Ridge dan Lasso, sehingga turunan pertamanya terdiri dari dua komponen.
1. Turunan Pertama \(\frac{\partial J}{\partial \beta}\)
Turunan dari bagian pertama:
\[\begin{equation} \frac{\partial}{\partial \beta} (y - X\beta)^T (y - X\beta) = -2X^T(y - X\beta). \label{eq:Elastic5} \end{equation}\]
Turunan dari penalti L1 berbeda untuk setiap elemen \(\beta_j\):
\[\begin{equation} \frac{\partial}{\partial \beta_j} \lambda_1 |\beta_j| = \lambda_1 \cdot \text{sign}(\beta_j), \label{eq:Elastic6} \end{equation}\]
dengan:
\[\begin{equation} \text{sign}(\beta_j) = \begin{cases} 1, & \text{jika } \beta_j > 0 \\ -1, & \text{jika } \beta_j < 0 \\ 0, & \text{jika } \beta_j = 0 \end{cases}. \end{equation}\]
Turunan dari penalti L2:
\[\begin{equation} \frac{\partial}{\partial \beta} \lambda_2 \beta^T \beta = 2\lambda_2 \beta. \label{eq:Elastic7} \end{equation}\]
Sehingga, turunan pertama dari \(J(\beta)\) adalah:
\[\begin{equation} \frac{\partial J}{\partial \beta} = -2X^T(y - X\beta) + \lambda_1 \cdot \text{sign}(\beta) + 2\lambda_2 \beta. \label{eq:Elastic8} \end{equation}\]
2. Menyamakan dengan Nol untuk Solusi Optimal
Menyusun ulang:
\[\begin{equation} X^{T}(X\beta - y) + \lambda_{1}\,\partial \|\beta\|_{1} + \lambda_{2}\beta = 0, \label{eq:Elastic9} \end{equation}\]
Karena fungsi nilai absolut tidak terdiferensiasi di \(\beta_j = 0\), solusi parameter harus ditemukan dengan metode numerik seperti Coordinate Descent.
Algoritma Estimasi Parameter: Coordinate Descent
Algoritma Coordinate Descent sering digunakan untuk menghitung solusi Elastic Net. Langkah-langkahnya adalah:
\[\begin{equation} \beta_j^{(t+1)} = S\left( \frac{X_j^T (y - X\beta_{-j})}{X_j^T X_j + \lambda_2}, \lambda_1 \right), \end{equation}\]
di mana \(S(z, \lambda_1)\) adalah fungsi soft-thresholding:
\[\begin{equation} S(z, \lambda_1) = \text{sign}(z) \max(|z| - \lambda_1, 0). \end{equation}\]
Ilustrasi Elastic Net dalam Kehadiran Multicollinearity
Tujuan dari ilustrasi ini adalah untuk menunjukkan bagaimana Elastic Net Regression dapat menangani multicollinearity lebih baik dibandingkan OLS (Ordinary Least Squares) dengan menggabungkan keunggulan Lasso dan Ridge Regression.
Simulasi Data dengan Multicollinearity
library(MASS)
library(glmnet)
library(ggplot2)
library(dplyr)
set.seed(123)
n <- 500 # Jumlah sampel
p <- 50 # Jumlah prediktor
# Membuat matriks kovarians dengan korelasi tinggi
Sigma <- matrix(0.9, nrow = p, ncol = p) + diag(0.1, p)
X <- mvrnorm(n, mu = rep(0, p), Sigma = Sigma) # Matriks prediktor
# Koefisien sebenarnya
beta_true <- runif(p, -2, 2)
y <- X %*% beta_true + rnorm(n, sd = 5) # Variabel respon
data <- data.frame(y, X)
Estimasi Model OLS dan Elastic Net Regression
# Model OLS
ols_model <- lm(y ~ ., data = data)
ols_coeff <- coef(ols_model)[-1] # Mengabaikan intercept
# Model Elastic Net Regression
x_matrix <- as.matrix(data[, -1])
y_vector <- data$y
elastic_net_model <- glmnet(x_matrix, y_vector, alpha = 0.5) # alpha = 0.5 untuk Elastic Net
cv_elastic_net <- cv.glmnet(x_matrix, y_vector, alpha = 0.5) # Cross-validation untuk lambda optimal
elastic_net_coeff <- coef(cv_elastic_net, s = "lambda.min")[-1] # Koefisien dengan lambda optimal
Visualisasi Model Elastic Net Regression
Plot Koefisien Elastic Net Regression terhadap lambda
par(mfrow = c(1,1)) # Pastikan hanya ada satu plot dalam satu area
plot(elastic_net_model,
xvar = "lambda", label = TRUE, cex.axis = 0.7, las = 2) # Menyesuaikan ukuran label
Interpretasi Plot Koefisien Elastic Net Regression terhadap lambda - Elastic Net Regression menggabungkan sifat shrinkage Ridge dan seleksi fitur Lasso. - Pada bagian kiri grafik (nilai lambda kecil), beberapa koefisien masih memiliki nilai besar. - Seiring bertambahnya lambda ke kanan, beberapa koefisien menjadi nol, sementara yang lain tetap diperkecil tetapi tidak dihilangkan sepenuhnya. - Elastic Net lebih fleksibel dibandingkan Lasso dan Ridge, karena dapat menyesuaikan shrinkage dan seleksi fitur berdasarkan alpha.
Cross-Validation Elastic Net Regression
plot(cv_elastic_net)
Interpretasi Plot Cross-Validation Elastic Net Regression
Plot Mean Squared Error untuk Berbagai Alpha
lambda_seq <- 10^seq(-5, 2, length.out = 100)
cv_fit_1 <- cv.glmnet(x_matrix, y_vector, alpha = 1, lambda = lambda_seq)
cv_fit_05 <- cv.glmnet(x_matrix, y_vector, alpha = 0.5, lambda = lambda_seq)
cv_fit_0 <- cv.glmnet(x_matrix, y_vector, alpha = 0, lambda = lambda_seq)
plot(log(cv_fit_1$lambda), cv_fit_1$cvm, col = "red", pch = 19,
xlab = "log(lambda)", ylab = "Mean-Squared Error")
points(log(cv_fit_05$lambda), cv_fit_05$cvm, col = "gray", pch = 19)
points(log(cv_fit_0$lambda), cv_fit_0$cvm, col = "blue", pch = 19)
legend("topright", legend = c("alpha=1 (Lasso)", "alpha=.5 (Elastic Net)",
"alpha 0 (Ridge)"),
col = c("red", "gray", "blue"), pch = 19)
Interpretasi Perbandingan Alpha dalam Cross-Validation
Perbandingan Koefisien OLS vs Elastic Net
comparison <- data.frame(
Predictor = paste0("X", 1:p),
OLS = ols_coeff,
ElasticNet = elastic_net_coeff
)
comparison_long <- comparison %>%
tidyr::pivot_longer(cols = c("OLS", "ElasticNet"),
names_to = "Model", values_to = "Coefficient")
# Visualisasi perbandingan koefisien
ggplot(comparison_long, aes(x = Predictor, y = Coefficient, fill = Model)) +
geom_bar(stat = "identity", position = "dodge") +
theme_minimal() +
labs(title = "Perbandingan Koefisien OLS vs Elastic Net Regression",
x = "Predictor",
y = "Estimated Coefficient") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
Kesimpulan
Dalam dunia machine learning, kita jarang menemukan satu model yang bekerja sempurna untuk semua jenis data. Oleh karena itu, pendekatan ensemble modeling menjadi sangat penting. Ensemble modeling adalah metode yang menggabungkan beberapa model dasar (base learners) untuk menghasilkan model yang lebih kuat dan akurat. Tujuan utamanya adalah untuk meningkatkan performa prediksi, baik dalam konteks klasifikasi maupun regresi.
Tiga teknik ensemble yang paling umum adalah Bagging, Boosting, dan Stacking. Ketiganya memiliki pendekatan yang berbeda terhadap kombinasi model dan memiliki kelebihan masing-masing.
Apa itu Ensemble Learning?
Ensemble learning adalah strategi dalam machine learning yang bertujuan menggabungkan prediksi dari beberapa model untuk meningkatkan akurasi dan generalisasi. Metode ini berfungsi dengan prinsip bahwa sekumpulan model yang berbeda (atau serupa) dapat mengompensasi kelemahan satu sama lain dan menghasilkan prediksi akhir yang lebih stabil.
Keuntungan Ensemble Learning:
Bagging adalah metode ensemble yang bertujuan untuk mengurangi varians model. Teknik ini bekerja dengan membuat beberapa subset data dari dataset pelatihan asli menggunakan bootstrap sampling (sampel dengan pengembalian), kemudian melatih model pada setiap subset tersebut secara paralel. Prediksi akhir ditentukan dengan agregasi (rata-rata untuk regresi, voting untuk klasifikasi).
Langkah-langkah:
Contoh: Random Forest
Random Forest adalah contoh terkenal dari metode Bagging yang menggunakan banyak decision tree dan rata-rata/voting sebagai agregator.
Implementasi di R
library(randomForest)
data(iris)
set.seed(123)
model_rf <- randomForest(Species ~ ., data = iris, ntree = 100)
print(model_rf)
##
## Call:
## randomForest(formula = Species ~ ., data = iris, ntree = 100)
## Type of random forest: classification
## Number of trees: 100
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 4.67%
## Confusion matrix:
## setosa versicolor virginica class.error
## setosa 50 0 0 0.00
## versicolor 0 47 3 0.06
## virginica 0 4 46 0.08
Konsep Dasar
Boosting adalah teknik ensemble yang bekerja secara berurutan. Model dibangun satu per satu, di mana setiap model baru mencoba memperbaiki kesalahan model sebelumnya. Tujuan utamanya adalah mengurangi bias model.
Jenis-Jenis Boosting:
Inti Proses Boosting:
Implementasi: Gradient Boosting
library(gbm)
data(mtcars)
set.seed(123)
boost_model <- gbm(
formula = mpg ~ .,
data = mtcars,
distribution = "gaussian",
n.trees = 100,
interaction.depth = 3,
shrinkage = 0.1,
bag.fraction = 1, # Gunakan seluruh data
n.minobsinnode = 5 # Kurangi dari default 10
)
summary(boost_model)
## var rel.inf
## cyl cyl 32.020538014
## wt wt 30.854244585
## hp hp 16.224242840
## disp disp 15.873942069
## drat drat 2.699074609
## qsec qsec 2.159667478
## gear gear 0.165308602
## carb carb 0.002981803
## vs vs 0.000000000
## am am 0.000000000
Konsep Dasar
Stacking adalah metode ensemble yang menggabungkan beberapa model prediksi (disebut base learners) dan memanfaatkan model meta untuk menggabungkan output dari base learners.
Langkah-langkah:
Keuntungan:
Implementasi di R dengan caretEnsemble
library(caret)
library(caretEnsemble)
# Subset 2-class data
iris2 <- iris[iris$Species != "setosa", ]
iris2$Species <- factor(iris2$Species)
set.seed(123)
control <- trainControl(method = "cv", number = 5, savePredictions = "final", classProbs = TRUE)
models <- caretList(
Species ~ ., data = iris2, trControl = control,
methodList = c("rpart", "knn", "nb")
)
stack_model <- caretStack(models, method = "glm")
summary(stack_model)
## The following models were ensembled: rpart, knn, nb
##
## Model Importance:
## rpart knn nb
## 0.0505 0.5201 0.4295
##
## Model accuracy:
## model_name metric value sd
## <char> <char> <num> <num>
## 1: ensemble ROC 0.984 0.008944272
## 2: rpart Accuracy 0.910 0.065192024
## 3: knn Accuracy 0.960 0.054772256
## 4: nb Accuracy 0.940 0.082158384
Perbandingan Bagging, Boosting, dan Stacking
| Aspek | Bagging | Boosting | Stacking |
|---|---|---|---|
| Tujuan | Mengurangi varians | Mengurangi bias | Kombinasi fleksibel dari model |
| Proses | Paralel | Berurutan | Gabungan + model meta |
| Toleransi outlier | Tinggi | Rendah | Sedang |
| Kompleksitas | Sedang | Tinggi | Tinggi |
| Contoh Algoritma | Random Forest | AdaBoost, XGBoost, LightGBM | Blending, stacking dengan meta |
Studi Kasus dan Evaluasi Model
Ensemble models sangat populer karena sering memberikan hasil yang lebih baik dibandingkan single model. Namun, perlu evaluasi dengan:
Ensemble learning merupakan teknik yang sangat kuat untuk meningkatkan performa model prediksi. Dengan memahami perbedaan dan penggunaan Bagging, Boosting, dan Stacking, praktisi data dapat memilih pendekatan yang paling tepat sesuai dengan masalah dan dataset yang dihadapi. Ketiganya memiliki kekuatan tersendiri dan seringkali menjadi solusi utama dalam dunia nyata maupun kompetisi data science.
Random Forest adalah algoritma machine learning berbasis decision tree yang termasuk dalam metode ensemble learning. Ia dikenal karena keandalannya, kemampuannya menangani dataset besar, dan performa prediksi yang sangat baik baik dalam tugas klasifikasi maupun regresi.
Tujuan dari dokumen ini adalah memberikan pemahaman mendalam tentang Random Forest, mencakup konsep teoretis, formulasi matematis, hingga implementasi praktis menggunakan bahasa pemrograman R.
Konsep Dasar Random Forest
Random Forest dibangun dengan prinsip bahwa kombinasi banyak decision tree (pohon keputusan) dapat menghasilkan model yang lebih stabil dan akurat. Model ini merupakan bentuk dari metode bagging (bootstrap aggregating), dengan penambahan seleksi acak fitur pada tiap node.
Kenapa Random?
Dalam Random Forest, “random” merujuk pada dua hal:
Model Matematis
Misalkan:
Klasifikasi: \[ \hat{y}(x) = \text{mode}(h_1(x), h_2(x), ..., h_B(x)) \]
Regresi: \[ \hat{y}(x) = \frac{1}{B} \sum_{b=1}^B h_b(x) \]
Algoritma Random Forest
Implementasi di R
Library dan Data
library(randomForest)
library(tidyverse)
data(iris)
str(iris)
## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
Pelatihan Model
set.seed(123)
rf_model <- randomForest(Species ~ ., data = iris, ntree = 100, mtry = 2, importance = TRUE)
print(rf_model)
##
## Call:
## randomForest(formula = Species ~ ., data = iris, ntree = 100, mtry = 2, importance = TRUE)
## Type of random forest: classification
## Number of trees: 100
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 4.67%
## Confusion matrix:
## setosa versicolor virginica class.error
## setosa 50 0 0 0.00
## versicolor 0 46 4 0.08
## virginica 0 3 47 0.06
Evaluasi Model
pred <- predict(rf_model, iris)
table(Predicted = pred, Actual = iris$Species)
## Actual
## Predicted setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 50 0
## virginica 0 0 50
mean(pred == iris$Species)
## [1] 1
Out-of-Bag (OOB) Error
plot(rf_model, main = "OOB Error Rate")
Feature Importance
importance(rf_model)
## setosa versicolor virginica MeanDecreaseAccuracy
## Sepal.Length 2.334999 3.8525315 3.909442 5.331391
## Sepal.Width 1.873228 -0.2344981 2.931587 2.299106
## Petal.Length 9.062834 14.1250709 10.997897 13.197705
## Petal.Width 11.214854 14.2553890 16.251155 16.117169
## MeanDecreaseGini
## Sepal.Length 10.133410
## Sepal.Width 2.021934
## Petal.Length 39.955902
## Petal.Width 47.133554
varImpPlot(rf_model, main = "Feature Importance (MeanDecreaseGini)")
Eksperimen Sederhana: 3 Pohon
set.seed(42)
rf_small <- randomForest(Species ~ ., data = iris, ntree = 3, mtry = 2, keep.forest = TRUE)
predict(rf_small, iris, predict.all = TRUE)$individual[1:5, ]
## [,1] [,2] [,3]
## 1 "setosa" "setosa" "setosa"
## 2 "setosa" "setosa" "setosa"
## 3 "setosa" "setosa" "setosa"
## 4 "setosa" "setosa" "setosa"
## 5 "setosa" "setosa" "setosa"
Analisis Hasil
Berdasarkan feature importance, fitur yang paling berkontribusi terhadap prediksi adalah Petal.Width dan Petal.Length. Ini konsisten dengan analisis EDA sebelumnya pada data iris.
Penggunaan OOB error membantu dalam estimasi error tanpa memerlukan data validasi terpisah.
Kelebihan dan Kekurangan Random Forest
Kelebihan:
Kekurangan:
Studi Kasus Singkat: Prediksi Spesies Bunga
Dengan data iris, kita bisa memanfaatkan Random Forest untuk memprediksi spesies bunga berdasarkan 4 fitur morfologi. Model ini memberikan akurasi tinggi bahkan tanpa tuning parameter secara mendalam.
Perbandingan dengan Algoritma Lain
| Algoritma | Akurasi | Interpretabilitas | Waktu Latih | Tahan Overfitting |
|---|---|---|---|---|
| Decision Tree | Sedang | Tinggi | Cepat | Tidak |
| Random Forest | Tinggi | Rendah | Sedang | Ya |
| Logistic Regression | Sedang | Sedang | Cepat | Tidak |
| KNN | Sedang | Rendah | Lambat (prediksi) | Tidak |
Strategi
mtry yang optimal (gunakan tuneRF
untuk mencari)ntree untuk menghindari overfitting
dan waktu komputasiimportance = TRUE untuk interpretasi fiturRandom Forest merupakan model yang sangat powerful dan fleksibel. Ia cocok untuk sebagian besar kasus di dunia nyata karena sifatnya yang non-parametrik, robust, dan relatif mudah digunakan.
AdaBoost (Adaptive Boosting) adalah algoritma ensemble learning yang sangat populer dalam dunia machine learning, terutama untuk klasifikasi. Algoritma ini bekerja dengan cara membangun model-model lemah (weak learners), biasanya decision stump, secara berurutan dan memperkuat model sebelumnya dengan memberikan bobot lebih besar pada observasi yang sulit diklasifikasi.
Dokumen ini akan membahas teori dasar, formulasi matematis, dan implementasi AdaBoost menggunakan bahasa R.
Konsep Dasar AdaBoost
Berbeda dengan Random Forest yang bersifat paralel, AdaBoost bekerja secara sekuensial. Setiap model berikutnya fokus untuk memperbaiki kesalahan yang dilakukan oleh model sebelumnya.
Intuisi Kerja AdaBoost
Awalnya, semua observasi diberi bobot sama.
Model pertama dilatih dan evaluasi error-nya.
Observasi yang salah diklasifikasikan akan mendapatkan bobot lebih tinggi.
Model selanjutnya dilatih pada data dengan bobot baru.
Proses diulang untuk sejumlah iterasi tertentu.
Tahap proses Adaboost
Penjelasan Matematis Algoritma AdaBoost
Misal kita memiliki data pelatihan \((x_1, y_1), ..., (x_n, y_n)\), dengan \(y_i \in \{-1, +1\}\).
Langkah Algoritma:
Load Libraries and Dataset
library(adabag)
library(caret)
library(ggplot2)
data(iris)
iris2 <- subset(iris, Species != "setosa")
iris2$Species <- factor(iris2$Species)
Train-Test Split
set.seed(123)
index <- createDataPartition(iris2$Species, p = 0.7, list = FALSE)
train <- iris2[index, ]
test <- iris2[-index, ]
Train AdaBoost Model
set.seed(123)
model_ada <- boosting(Species ~ ., data = train, boos = TRUE, mfinal = 10)
summary(model_ada)
## Length Class Mode
## formula 3 formula call
## trees 10 -none- list
## weights 10 -none- numeric
## votes 140 -none- numeric
## prob 140 -none- numeric
## class 70 -none- character
## importance 4 -none- numeric
## terms 3 terms call
## call 5 -none- call
Plot Decision Boundary
grid <- expand.grid(
Petal.Length = seq(min(iris2$Petal.Length), max(iris2$Petal.Length), length.out = 100),
Petal.Width = seq(min(iris2$Petal.Width), max(iris2$Petal.Width), length.out = 100)
)
grid$Sepal.Length <- mean(iris2$Sepal.Length)
grid$Sepal.Width <- mean(iris2$Sepal.Width)
# Predict on grid
grid_pred <- predict.boosting(model_ada, newdata = grid)
grid$pred <- as.factor(grid_pred$class)
# Plot
ggplot(train, aes(x = Petal.Length, y = Petal.Width, color = Species)) +
geom_point(size = 2) +
geom_tile(data = grid, aes(fill = pred), alpha = 0.2, color = NA) +
labs(title = "Decision Boundary AdaBoost", x = "Petal Length", y = "Petal Width") +
theme_minimal()
Evaluasi Model dan Interpretasi
pred_test <- predict.boosting(model_ada, newdata = test)
confusionMatrix(as.factor(pred_test$class), test$Species)
## Confusion Matrix and Statistics
##
## Reference
## Prediction versicolor virginica
## versicolor 14 1
## virginica 1 14
##
## Accuracy : 0.9333
## 95% CI : (0.7793, 0.9918)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : 4.34e-07
##
## Kappa : 0.8667
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9333
## Specificity : 0.9333
## Pos Pred Value : 0.9333
## Neg Pred Value : 0.9333
## Prevalence : 0.5000
## Detection Rate : 0.4667
## Detection Prevalence : 0.5000
## Balanced Accuracy : 0.9333
##
## 'Positive' Class : versicolor
##
Interpretasi:
Plot Error vs Iterasi
errores <- errorevol(model_ada, newdata = test)
plot(errores$error, type = "b", main = "Test Error vs Iterasi", xlab = "Iterasi", ylab = "Test Error")
Prediksi Data Baru
new_obs <- data.frame(
Sepal.Length = 6.0,
Sepal.Width = 2.8,
Petal.Length = 5.0,
Petal.Width = 1.8
)
new_pred <- predict.boosting(model_ada, newdata = new_obs)
new_pred$class # hasil klasifikasi
## [1] "virginica"
new_pred$prob # probabilitas
## [,1] [,2]
## [1,] 0.1255672 0.8744328
adabag.Untuk aplikasi lanjut, bisa digabungkan dengan metode pemilihan fitur atau digunakan dalam ensemble yang lebih besar.
Keunggulan dan Kelemahan AdaBoost
Keunggulan:
Kelemahan:
Perbandingan AdaBoost vs Random Forest
| Aspek | AdaBoost | Random Forest |
|---|---|---|
| Proses | Berurutan | Paralel |
| Fokus | Perbaiki kesalahan | Kurangi varians |
| Sensitivitas noise | Tinggi | Rendah |
| Interpretasi | Sedang | Rendah |
Tuning Parameter
iter: jumlah iterasi boostingnu: learning rate (biasanya 0 < nu ≤ 1)type: “discrete”, “real”, atau “gentle”Gunakan train() dari caret untuk tuning
otomatis.
AdaBoost adalah teknik boosting yang sangat efisien dalam meningkatkan akurasi model klasifikasi. Implementasinya sederhana, dan dapat diterapkan ke banyak jenis data. Namun, penggunaannya harus memperhatikan outlier dan distribusi noise karena sifat algoritma yang adaptif.
Dalam materi ini, kita akan membahas tiga teknik boosting dalam machine learning:
Tujuan pembelajaran:
Konsep Dasar Boosting
Boosting adalah metode ensemble learning yang membentuk model prediksi kuat dari banyak model lemah (weak learners), biasanya berupa decision tree sederhana (shallow trees).
Pendahuluan
Gradient Boosting adalah metode boosting yang sangat kuat dan fleksibel, digunakan baik untuk klasifikasi maupun regresi. Algoritma ini bekerja dengan membangun model secara berurutan, di mana setiap model baru mempelajari kesalahan residu dari model sebelumnya, menggunakan pendekatan gradient descent terhadap loss function.
Dokumen ini membahas teori dasar, formulasi matematis, serta implementasi Gradient Boosting menggunakan R.
Konsep Dasar Gradient Boosting
Gradient Boosting termasuk dalam metode boosting, yaitu teknik ensemble learning yang menggabungkan banyak model lemah untuk membentuk model prediksi yang kuat. Fokus utamanya adalah mengurangi bias dari model dengan cara iteratif.
Proses Umum:
Model Matematis
Misalkan:
Untuk \(m = 1\) hingga \(M\):
\[ r_{im} = -\left[\frac{\partial L(y_i, F(x_i))}{\partial F(x_i)}\right]_{F(x)=F_{m-1}(x)} \]
Fit model \(h_m(x)\) ke \(r_{im}\)
Temukan langkah optimal \(\gamma_m\):
\[ \gamma_m = \arg\min_{\gamma} \sum L(y_i, F_{m-1}(x_i) + \gamma h_m(x_i)) \]
\[ F_m(x) = F_{m-1}(x) + \gamma_m h_m(x) \]
Implementasi di R
Library dan Data
library(gbm)
library(tidyverse)
data(mtcars)
Latih Model Gradient Boosting
set.seed(123)
gb_model <- gbm(
formula = mpg ~ .,
data = mtcars,
distribution = "gaussian",
n.trees = 100,
interaction.depth = 3,
shrinkage = 0.1,
bag.fraction = 1, # Gunakan seluruh data, hindari subsampling
n.minobsinnode = 3, # Kurangi minimum agar bisa training
cv.folds = 0 # Matikan CV, karena data terlalu kecil
)
Evaluasi Model
best_iter <- which.min(gb_model$train.error)
pred <- predict(gb_model, mtcars, n.trees = best_iter)
rmse <- sqrt(mean((mtcars$mpg - pred)^2))
cat("RMSE:", rmse)
## RMSE: 0.3054586
Visualisasi dan Interpretasi
plot(gb_model)
Eksperimen Learning Rate
set.seed(123)
gb_small_eta <- gbm(
formula = mpg ~ .,
data = mtcars,
distribution = "gaussian",
n.trees = 100,
interaction.depth = 3,
shrinkage = 0.01, # Lebih kecil
bag.fraction = 1, # Gunakan seluruh data
n.minobsinnode = 3 # Minimum node kecil agar training jalan
)
gb_large_eta <- gbm(
formula = mpg ~ .,
data = mtcars,
distribution = "gaussian",
n.trees = 100,
interaction.depth = 3,
shrinkage = 0.3, # learning rate besar
bag.fraction = 1, # gunakan seluruh data
n.minobsinnode = 3 # agar training bisa jalan
)
Evaluasi model
set.seed(123)
best_iter_large <- which.min(gb_large_eta$train.error)
pred_large <- predict(gb_large_eta, mtcars, n.trees = best_iter_large)
rmse_large <- sqrt(mean((mtcars$mpg - pred_large)^2))
cat("RMSE (eta = 0.3):", rmse_large)
## RMSE (eta = 0.3): 0.03880949
Kelebihan dan Kelemahan Gradient Boosting
Kelebihan:
Kelemahan:
Perbandingan dengan Metode Lain
| Aspek | Gradient Boosting | Random Forest | AdaBoost |
|---|---|---|---|
| Proses | Berurutan (gradient) | Paralel | Berurutan (reweight) |
| Learning Rate | Ya | Tidak | Ya |
| Overfitting | Bisa terjadi | Lebih tahan | Bisa terjadi |
| Interpretasi | Sedang | Rendah | Sedang |
Strategi
cv.folds untuk memilih n.trees
optimalPenutup
Gradient Boosting adalah salah satu algoritma prediksi yang sangat fleksibel dan akurat, namun membutuhkan perhatian khusus dalam tuning parameter. Ia cocok untuk data dengan kompleksitas tinggi dan cocok dalam banyak kompetisi data science.
Pendahuluan
Extreme Gradient Boosting (XGBoost) adalah salah satu algoritma boosting paling populer dan efisien. Dikenal karena kecepatannya, akurasinya, dan kemampuannya menangani data besar, XGBoost telah menjadi algoritma andalan dalam banyak kompetisi data science seperti Kaggle.
Dalam dokumen ini, kita akan membahas teori dasar, formulasi matematis, dan implementasi XGBoost menggunakan bahasa R.
Konsep Dasar XGBoost
XGBoost merupakan pengembangan dari Gradient Boosting yang mengintegrasikan beberapa peningkatan penting:
Model Matematis
Tujuan utama XGBoost adalah meminimalkan fungsi objektif sebagai berikut:
\[ Obj = \sum_{i=1}^n l(y_i, \hat{y}_i) + \sum_{k=1}^K \Omega(f_k) \]
dengan:
Untuk setiap iterasi:
\[ \hat{y}_i^{(t)} = \hat{y}_i^{(t-1)} + f_t(x_i) \]
XGBoost menggunakan Taylor expansion hingga orde ke-2 untuk mengaproksimasi loss function.
Implementasi di R
Library dan Data
library(xgboost)
data(iris)
iris_x <- as.matrix(iris[, -5])
iris_y <- as.numeric(iris$Species) - 1
dtrain <- xgb.DMatrix(data = iris_x, label = iris_y)
Latih Model XGBoost
params <- list(
objective = "multi:softprob",
num_class = 3,
eta = 0.1,
max_depth = 3,
eval_metric = "mlogloss"
)
set.seed(123)
xgb_model <- xgb.train(
params = params,
data = dtrain,
nrounds = 50,
verbose = 0
)
Evaluasi Model
pred <- predict(xgb_model, iris_x)
pred_labels <- max.col(matrix(pred, ncol = 3, byrow = TRUE)) - 1
table(Predicted = pred_labels, Actual = iris_y)
## Actual
## Predicted 0 1 2
## 0 50 0 0
## 1 0 50 0
## 2 0 0 50
mean(pred_labels == iris_y)
## [1] 1
Feature Importance
importance <- xgb.importance(model = xgb_model)
xgb.plot.importance(importance_matrix = importance)
Eksperimen Parameter
params2 <- list(
objective = "multi:softprob",
num_class = 3,
eta = 0.05,
max_depth = 5,
subsample = 0.7,
colsample_bytree = 0.8
)
Kelebihan dan Kelemahan XGBoost
Kelebihan:
Kelemahan:
Perbandingan XGBoost dengan Metode Lain
| Aspek | XGBoost | Gradient Boosting | Random Forest |
|---|---|---|---|
| Paralelisasi | Ya | Tidak | Ya |
| Regularisasi | L1, L2 | Tidak eksplisit | Tidak eksplisit |
| Akurasi | Sangat tinggi | Tinggi | Tinggi |
| Waktu Latih | Cepat | Lambat | Sedang |
Strategi - Gunakan xgb.cv untuk mencari jumlah iterasi
optimal - Tuning parameter: eta, max_depth,
subsample, colsample_bytree - Gunakan
early_stopping_rounds untuk menghindari overfitting
XGBoost adalah algoritma boosting yang sangat kuat dan telah menjadi andalan dalam berbagai aplikasi. Dengan regularisasi, kontrol overfitting, dan kemampuan skalabilitas, XGBoost menjadi pilihan utama untuk tugas klasifikasi dan regresi dengan data kompleks.
Referensi - Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning. Springer. - Hoerl, A. E., & Kennard, R. W. (1970). Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1), 55-67. - Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B, 58(1), 267-288. - Zou, H., & Hastie, T. (2005). Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society: Series B, 67(2), 301-320.
Pendekatan Bayesian menggabungkan informasi sebelum melihat data (prior) dengan likelihood dari data untuk menghasilkan posterior, yaitu keyakinan tentang parameter setelah melihat data. Inti dari pendekatan ini adalah Teorema Bayes: \[ p(\theta\mid y) \propto p(y\mid \theta)\,p(\theta), \] di mana \(p(y\mid\theta)\) adalah likelihood, \(p(\theta)\) adalah prior, dan \(p(\theta\mid y)\) adalah posterior.
Keunggulan Bayesian:
Likelihood \(L(\theta; y) = p(y\mid\theta)\) mengukur seberapa mungkin data yang kita amati muncul di bawah parameter \(\theta\). Contoh:
Dalam inferensi Bayesian, kita tidak memaksimalkan likelihood, tetapi menggabungkan dengan prior untuk memperoleh posterior.
Prior mencerminkan keyakinan awal tentang parameter sebelum melihat data.
Ringkas: informative (kuat), weakly-informative (lemah namun membatasi ekor), vague (sangat difus), non-informative (minim info; bisa proper/improper).
Posterior menggabungkan prior dan likelihood: \[ p(\theta\mid y) \propto p(y\mid \theta)\,p(\theta). \]
Ketika prior konjugat terhadap likelihood, bentuk posterior tetap dalam keluarga distribusi yang sama dan memiliki bentuk tertutup (closed-form). Keuntungannya: perhitungan ringkas (mean, varians, kredibel interval) dan mudah melakukan posterior predictive.
Misal \(y_1,\dots,y_n
\stackrel{iid}{\sim}\text{Poisson}(\lambda)\) dan prior \(\lambda \sim \text{Gamma}(\alpha_0,
\beta_0)\) (bentuk shape–rate, yaitu densitas
\(\propto \lambda^{\alpha-1}
e^{-\beta\lambda}\)).
Maka: \[
p(\lambda\mid y) \sim \text{Gamma}\!\Big(\alpha_n,\ \beta_n\Big),\quad
\alpha_n = \alpha_0 + \sum_{i=1}^n y_i,\quad
\beta_n = \beta_0 + n.
\] Beberapa ringkasan: \[
\mathbb{E}[\lambda\mid y] = \frac{\alpha_n}{\beta_n},\qquad
\mathbb{V}[\lambda\mid y] = \frac{\alpha_n}{\beta_n^2}.
\]
Posterior predictive untuk pengamatan baru \(\tilde{y}\mid y\) adalah Negative Binomial (parameterisasi size–prob): \[ \tilde{y}\mid y \sim \text{NegBin}\big(\text{size}=\alpha_n,\ \text{prob}=\tfrac{\beta_n}{\beta_n+1}\big). \]
Dokumen ini menjabarkan mengapa prior Gamma untuk
laju \(\lambda\) pada model
Poisson menghasilkan posterior yang
juga Gamma (konjugat). Kita menggunakan
parameterisasi rate untuk Gamma:
\(\text{Gamma}(\alpha, \beta)\) dengan
pdf \[
p(\lambda)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}\,\lambda^{\alpha-1}
e^{-\beta \lambda},\qquad \lambda>0.
\]
Misalkan \(y_1,\dots,y_n \mid \lambda \stackrel{i.i.d.}{\sim} \text{Poisson}(\lambda)\) dengan pmf \[ p(y_i\mid \lambda)=\frac{e^{-\lambda}\lambda^{y_i}}{y_i!},\qquad y_i\in\{0,1,2,\dots\}. \] Ambil prior \(\lambda \sim \text{Gamma}(\alpha_0,\beta_0)\): \[ p(\lambda)=\frac{\beta_0^{\alpha_0}}{\Gamma(\alpha_0)}\,\lambda^{\alpha_0-1} e^{-\beta_0\lambda},\quad \lambda>0. \]
Karena i.i.d., \[ L(\lambda; \mathbf y)=\prod_{i=1}^n p(y_i\mid \lambda) =\left(\prod_{i=1}^n \frac{1}{y_i!}\right) e^{-n\lambda} \lambda^{\sum_{i=1}^n y_i}. \] Abaikan konstanta yang tidak bergantung pada \(\lambda\), definisikan \(S=\sum_{i=1}^n y_i\): \[ L(\lambda; \mathbf y)\ \propto\ e^{-n\lambda}\,\lambda^{S}. \]
\[ p(\lambda\mid \mathbf y)\ \propto\ p(\lambda)\,L(\lambda;\mathbf y) \ \propto\ \Big[\lambda^{\alpha_0-1} e^{-\beta_0\lambda}\Big]\Big[e^{-n\lambda}\lambda^{S}\Big] = \lambda^{(\alpha_0-1)+S} e^{-(\beta_0+n)\lambda}. \] Tetapkan \[ \alpha_n := \alpha_0 + S,\qquad \beta_n := \beta_0 + n. \] Maka kernel posterior \[ p(\lambda\mid \mathbf y)\ \propto\ \lambda^{\alpha_n-1} e^{-\beta_n\lambda}, \] yang identik dengan Gamma(\(\alpha_n,\beta_n\)) (parameterisasi rate). Dengan konstanta normalisasi, pdf lengkapnya: \[ p(\lambda\mid \mathbf y) = \frac{\beta_n^{\alpha_n}}{\Gamma(\alpha_n)}\,\lambda^{\alpha_n-1} e^{-\beta_n\lambda},\quad \alpha_n=\alpha_0+\textstyle\sum_i y_i,\ \beta_n=\beta_0+n. \]
Gunakan identitas \(\int_0^{\infty}
\lambda^{a-1} e^{-b\lambda}\,d\lambda = \Gamma(a)/b^a\) untuk
\(a>0,\ b>0\).
Sehingga konstanta normalisasi adalah \(\beta_n^{\alpha_n}/\Gamma(\alpha_n)\).
Jika \(\lambda\mid\mathbf y\sim\text{Gamma}(\alpha_n,\beta_n)\) (rate):
qgamma(q, shape=alpha_n, rate=beta_n).Untuk hitungan baru \(y_{new}\): \[ p(y_{new}\mid \mathbf y) = \int p(y_{new}\mid \lambda)\,p(\lambda\mid \mathbf y)\,d\lambda = \operatorname{NegBin}\!\left(r=\alpha_n,\ p=\tfrac{\beta_n}{\beta_n+1}\right). \] Ini berasal dari campuran Poisson–Gamma.
Contoh kecil untuk memeriksa \(\alpha_n,\beta_n\) dan ringkasan posterior:
set.seed(1)
# Data mainan: 30 hari, laju ~ 3
n <- 30
lambda_true <- 3
y <- rpois(n, lambda_true)
S <- sum(y)
# Prior Gamma(rate)
alpha0 <- 2; beta0 <- 1 # proper, tidak terlalu informatif
alpha_n <- alpha0 + S
beta_n <- beta0 + n
c(alpha_n = alpha_n, beta_n = beta_n,
mean_post = alpha_n / beta_n,
var_post = alpha_n / beta_n^2)
## alpha_n beta_n mean_post var_post
## 95.00000000 31.00000000 3.06451613 0.09885536
ql <- qgamma(0.025, shape = alpha_n, rate = beta_n)
qu <- qgamma(0.975, shape = alpha_n, rate = beta_n)
c(q025 = ql, q975 = qu)
## q025 q975
## 2.479377 3.710716
Catatan Parameterisasi Jika menggunakan parameterisasi scale \(\theta=1/\beta\), maka posterior tetap Gamma dengan \(\theta_n = 1/(\beta_0+n)\) dan bentuk pdf menyesuaikan. Pastikan konsisten antara rate vs scale saat memakai fungsi
dgamma/qgamma/rgamma.
set.seed(123)
# Data: hitung kejadian per hari (misal 30 hari)
n <- 30
lambda_true <- 3.0
y <- rpois(n, lambda_true)
sum_y <- sum(y)
# PRIOR: tiga skenario
priors <- list(
informative = list(alpha0 = 15, beta0 = 5), # mean 3, sd ~0.77
weakly = list(alpha0 = 3, beta0 = 1), # mean 3, sd ~1.73
vague = list(alpha0 = 0.1, beta0 = 0.1) # sangat difus (tetap proper)
)
post_tab <- data.frame(
prior = character(),
alpha0 = numeric(), beta0 = numeric(),
alpha_n = numeric(), beta_n = numeric(),
mean_post = numeric(), sd_post = numeric(),
q025 = numeric(), q975 = numeric(),
stringsAsFactors = FALSE
)
for (nm in names(priors)) {
a0 <- priors[[nm]]$alpha0
b0 <- priors[[nm]]$beta0
an <- a0 + sum_y
bn <- b0 + n
mn <- an / bn
sd <- sqrt(an) / bn
ql <- qgamma(0.025, shape = an, rate = bn)
qu <- qgamma(0.975, shape = an, rate = bn)
post_tab <- rbind(post_tab, data.frame(
prior = nm, alpha0 = a0, beta0 = b0,
alpha_n = an, beta_n = bn,
mean_post = mn, sd_post = sd,
q025 = ql, q975 = qu
))
}
post_tab
## prior alpha0 beta0 alpha_n beta_n mean_post sd_post q025 q975
## 1 informative 15.0 5.0 117.0 35.0 3.342857 0.3090473 2.764634 3.975170
## 2 weakly 3.0 1.0 105.0 31.0 3.387097 0.3305468 2.770310 4.064949
## 3 vague 0.1 0.1 102.1 30.1 3.392027 0.3356962 2.766083 4.080859
## MLE untuk λ pada model Poisson i.i.d. harian --
## Fakta: S = sum(y) ~ Poisson(n * λ)
## MLE: λ_hat = mean(y) = S / n
## Var(λ_hat) ≈ λ / n → SE ≈ sqrt(λ_hat / n)
## CI eksak untuk μ = nλ: [0.5 * χ²_{α/2, 2S}, 0.5 * χ²_{1-α/2, 2(S+1)}], lalu dibagi n
alpha <- 0.05
S <- sum_y
lambda_mle <- mean(y)
se_mle <- sqrt(lambda_mle / n)
# CI Wald (asymptotik normal)
z <- qnorm(1 - alpha/2)
ci_wald <- c(lambda_mle - z * se_mle, lambda_mle + z * se_mle)
ci_wald[ci_wald < 0] <- 0 # potong di 0 bila perlu
# CI eksak (berbasis Chi-square pada total μ = nλ)
# perlakuan tepi bila S = 0
if (S == 0) {
mu_low <- 0
mu_upp <- 0.5 * qchisq(1 - alpha/2, df = 2*(S + 1))
} else {
mu_low <- 0.5 * qchisq(alpha/2, df = 2*S)
mu_upp <- 0.5 * qchisq(1 - alpha/2, df = 2*(S + 1))
}
ci_exact <- c(mu_low / n, mu_upp / n)
# Log-likelihood pada λ_hat (abaikan konstanta lfactorial)
loglik <- sum(dpois(y, lambda_mle, log = TRUE))
k <- 1L
AIC_mle <- 2*k - 2*loglik
mle_tab <- data.frame(
method = "MLE (Poisson)",
lambda_mle = lambda_mle,
se = se_mle,
ci_wald_lo = ci_wald[1], ci_wald_hi = ci_wald[2],
ci_exact_lo = ci_exact[1], ci_exact_hi = ci_exact[2],
logLik = loglik,
AIC = AIC_mle
)
mle_tab
## method lambda_mle se ci_wald_lo ci_wald_hi ci_exact_lo
## 1 MLE (Poisson) 3.4 0.3366502 2.740178 4.059822 2.772294
## ci_exact_hi logLik AIC
## 1 4.127368 -60.35731 122.7146
# Ringkasan Bayesian (dari post_tab) + satu baris MLE
compare_tab <- rbind(
setNames(
data.frame(
prior = post_tab$prior,
estimator = "Posterior (mean)",
point = post_tab$mean_post,
lo = post_tab$q025,
hi = post_tab$q975,
stringsAsFactors = FALSE
),
c("prior","estimator","point","lo","hi")
),
data.frame(
prior = "—",
estimator = "MLE",
point = mle_tab$lambda_mle,
lo = mle_tab$ci_exact_lo, # pakai CI eksak agar konservatif
hi = mle_tab$ci_exact_hi,
stringsAsFactors = FALSE
)
)
compare_tab
## prior estimator point lo hi
## 1 informative Posterior (mean) 3.342857 2.764634 3.975170
## 2 weakly Posterior (mean) 3.387097 2.770310 4.064949
## 3 vague Posterior (mean) 3.392027 2.766083 4.080859
## 4 — MLE 3.400000 2.772294 4.127368
# Pilih skenario untuk divisualisasikan
nm <- "weakly"
a0 <- priors[[nm]]$alpha0
b0 <- priors[[nm]]$beta0
an <- a0 + sum_y
bn <- b0 + n
# Grid lambda
lam <- seq(0, max(8, lambda_true*3), length.out = 400)
# Plot prior & posterior (Gamma shape-rate)
plot(lam, dgamma(lam, a0, b0), type = "l", lwd = 2,
xlab = expression(lambda), ylab = "densitas",
main = paste("Prior vs Posterior (", nm, ")", sep = ""))
lines(lam, dgamma(lam, an, bn), lwd = 2, lty = 2)
legend("topright", legend = c("Prior (Gamma a0,b0)", "Posterior (Gamma an,bn)"),
lwd = 2, lty = c(1,2), bty = "n")
abline(v = lambda_true, col = "gray50", lty = 3)
# Parameter NB: size = an, prob = bn/(bn+1)
size <- an
prob <- bn/(bn+1)
# Distribusi prediktif
k <- 0:15
pmf <- dnbinom(k, size = size, prob = prob)
plot(k, pmf, type = "h", lwd = 3, xlab = expression(tilde(y)), ylab = "Probabilitas",
main = "Posterior Predictive: 1 hari ke depan")
points(k, pmf, pch = 16)
Catatan
Regresi Bayesian adalah pendekatan inferensi yang menggabungkan prior information dengan data observasi melalui Teorema Bayes. Fokus utama adalah pada distribusi posterior dari parameter, yang mencerminkan ketidakpastian parameter setelah melihat data.
Motivasi utama: - Integrasi informasi sebelumnya (prior) dan data saat ini (likelihood). - Menghasilkan credible interval yang memiliki interpretasi probabilistik langsung. - Fleksibel untuk memperluas ke model kompleks (hierarkis, spasial, survival, dll.).
Misalkan model regresi linear Gaussian:
\[ y_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip} + \varepsilon_i,\quad \varepsilon_i \sim \mathcal{N}(0,\sigma^2). \]
Dalam kerangka Bayesian, tentukan: - Likelihood (untuk vektor \(y\) berukuran \(n\) dan matriks desain \(X\) berukuran \(n\times p\)): \[ p(y\mid \beta,\sigma^2) \propto (\sigma^2)^{-\frac{n}{2}}\exp\Big\{-\tfrac{1}{2\sigma^2}(y-X\beta)'(y-X\beta)\Big\}. \]
Prior (Normal–Inverse-Gamma yang konjugat): \[ \beta\mid\sigma^2 \sim \mathcal{N}(\mu_0,\, \sigma^2\,\Sigma_0),\qquad \sigma^2 \sim \text{IG}(a_0,b_0). \]
Posterior (tertutup bentuknya untuk prior di atas): \[ \begin{aligned} \Sigma_n &= (\Sigma_0^{-1} + X'X)^{-1} \\ \mu_n &= \Sigma_n(\Sigma_0^{-1}\mu_0 + X'y) \\ a_n &= a_0 + \tfrac{n}{2} \\ b_n &= b_0 + \tfrac{1}{2}\big(y' y + \mu_0'\Sigma_0^{-1}\mu_0 - \mu_n'\Sigma_n^{-1}\mu_n\big) \end{aligned} \]
Distribusi bersyaratnya: \(\beta\mid\sigma^2,y\sim\mathcal{N}(\mu_n,\sigma^2\Sigma_n)\) dan \(\sigma^2\mid y\sim \text{IG}(a_n,b_n)\).
Catatan: Formulasi ini mudah diperluas ke prior ridge/lasso Bayesian (mis. prior Laplace untuk lasso) dan model hierarkis.
Markov Chain Monte Carlo (MCMC) membangkitkan sampel dari posterior ketika bentuk tertutup sulit/ tak tersedia. Dua skema umum: - Gibbs sampling: sampling bergantian dari sebaran bersyarat lengkap. - Metropolis–Hastings (MH): menerima/menolak proposal dari sebaran kandidat.
# Paket
suppressPackageStartupMessages({
library(coda)
library(rjags) # install.packages('rjags') jika perlu
})
set.seed(123)
n <- 200
x1 <- rnorm(n, 0, 1)
x2 <- rnorm(n, 0, 1)
y <- 1 + 2*x1 - 1.5*x2 + rnorm(n, 0, 1)
dat <- data.frame(y, x1, x2)
summary(dat)
## y x1 x2
## Min. :-5.6760 Min. :-2.30917 Min. :-2.46590
## 1st Qu.:-1.0130 1st Qu.:-0.62576 1st Qu.:-0.59077
## Median : 1.1613 Median :-0.05874 Median : 0.02283
## Mean : 0.9515 Mean :-0.00857 Mean : 0.04212
## 3rd Qu.: 2.8011 3rd Qu.: 0.56840 3rd Qu.: 0.71482
## Max. : 9.0915 Max. : 3.24104 Max. : 2.57146
data_jags <- list(
y = dat$y,
x1 = dat$x1,
x2 = dat$x2,
n = nrow(dat)
)
model_string <- "
model {
for (i in 1:n) {
y[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i]
}
# Prior diffuse (Normal(0, var=100))
beta0 ~ dnorm(0, 0.01)
beta1 ~ dnorm(0, 0.01)
beta2 ~ dnorm(0, 0.01)
# Prior untuk precision (tau): Gamma(0.01, 0.01)
tau ~ dgamma(0.01, 0.01)
sigma2 <- 1 / tau
}
"
jm <- jags.model(textConnection(model_string), data = data_jags, n.chains = 3, n.adapt = 1000)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 200
## Unobserved stochastic nodes: 4
## Total graph size: 1209
##
## Initializing model
update(jm, 2000) # burn-in
samps <- coda.samples(jm, variable.names = c("beta0","beta1","beta2","sigma2"), n.iter = 6000, thin = 2)
summary(samps)
##
## Iterations = 2002:8000
## Thinning interval = 2
## Number of chains = 3
## Sample size per chain = 3000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta0 1.0335 0.06997 0.0007375 0.0007478
## beta1 1.9668 0.07330 0.0007726 0.0007727
## beta2 -1.5504 0.06954 0.0007330 0.0007112
## sigma2 0.9467 0.09581 0.0010099 0.0010357
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta0 0.8960 0.9862 1.0339 1.080 1.172
## beta1 1.8240 1.9180 1.9666 2.016 2.109
## beta2 -1.6863 -1.5974 -1.5503 -1.505 -1.411
## sigma2 0.7798 0.8784 0.9405 1.008 1.151
plot(samps) # trace & density
gelman.diag(samps)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## beta0 1 1
## beta1 1 1
## beta2 1 1
## sigma2 1 1
##
## Multivariate psrf
##
## 1
autocorr.plot(samps)
sumtab <- summary(samps)
sumtab$statistics
## Mean SD Naive SE Time-series SE
## beta0 1.0334937 0.06996958 0.0007375441 0.0007478362
## beta1 1.9667796 0.07329995 0.0007726493 0.0007726678
## beta2 -1.5503510 0.06954030 0.0007330192 0.0007111687
## sigma2 0.9466601 0.09580949 0.0010099207 0.0010356970
sumtab$quantiles[c("beta0","beta1","beta2","sigma2"), c("2.5%","50%","97.5%")]
## 2.5% 50% 97.5%
## beta0 0.8960165 1.0338858 1.172175
## beta1 1.8239700 1.9666248 2.109275
## beta2 -1.6862561 -1.5503275 -1.411453
## sigma2 0.7798064 0.9404966 1.150941
# Ambil 1000 sampel posterior untuk prediksi
post <- as.matrix(samps)
ix <- sample(seq_len(nrow(post)), 1000, replace = TRUE)
# Prediksi untuk grid x1, x2 tetap pada rata-rata
newx1 <- seq(min(x1), max(x1), length.out = 100)
newx2 <- 0
Xnew <- cbind(1, newx1, newx2)
pred_draws <- sapply(ix, function(k){
b0 <- post[k, "beta0"]
b1 <- post[k, "beta1"]
b2 <- post[k, "beta2"]
s2 <- post[k, "sigma2"]
mu <- as.vector(Xnew %*% c(b0, b1, b2))
rnorm(length(mu), mu, sqrt(s2))
})
pred_mean <- rowMeans(pred_draws)
pred_low <- apply(pred_draws, 1, quantile, probs = 0.025)
pred_high <- apply(pred_draws, 1, quantile, probs = 0.975)
plot(newx1, pred_mean, type = "l", lwd = 2, ylab = "y", xlab = "x1")
lines(newx1, pred_low, lty = 2)
lines(newx1, pred_high, lty = 2)
points(x1, y, pch = 16, cex = 0.6)
INLA (Integrated Nested Laplace Approximation) adalah pendekatan deterministik untuk model Latent Gaussian. INLA mengevaluasi posterior marginal dari parameter dan latent field dengan aproksimasi Laplace bertingkat sehingga jauh lebih cepat daripada MCMC untuk banyak model umum (regresi, GLM, hierarkis, spasial, survival).
# install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/stable"))
suppressPackageStartupMessages(library(INLA))
res_inla <- inla(
y ~ x1 + x2,
data = dat,
family = "gaussian",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE),
control.fixed = list(
mean.intercept = 0, prec.intercept = 0.01, # prior Normal(0, var=100) -> prec=1/100=0.01
mean = 0, prec = 0.01
)
)
summary(res_inla)
## Time used:
## Pre = 0.994, Running = 0.704, Post = 0.0157, Total = 1.71
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 1.034 0.068 0.899 1.034 1.168 1.034 110.163
## x1 1.967 0.073 1.825 1.967 2.110 1.967 230.839
## x2 -1.551 0.069 -1.686 -1.551 -1.415 -1.551 179.623
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 1.08 0.108 0.877 1.07
## 0.975quant mode
## Precision for the Gaussian observations 1.30 1.07
##
## Deviance Information Criterion (DIC) ...............: 559.51
## Deviance Information Criterion (DIC, saturated) ....: 206.38
## Effective number of parameters .....................: 4.00
##
## Watanabe-Akaike information criterion (WAIC) ...: 559.92
## Effective number of parameters .................: 4.28
##
## Marginal log-Likelihood: -301.90
## CPO, PIT is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
res_inla$dic$dic; res_inla$waic$waic
## [1] 559.5077
## [1] 559.9199
# Simulasi data count
set.seed(42)
n2 <- 300
z1 <- rnorm(n2)
eta <- -0.5 + 0.8*z1
lam <- exp(eta)
ycount <- rpois(n2, lam)
df_pois <- data.frame(ycount, z1)
fit_pois <- inla(
ycount ~ z1,
data = df_pois,
family = "poisson",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE)
)
summary(fit_pois)
## Time used:
## Pre = 0.916, Running = 0.196, Post = 0.00655, Total = 1.12
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.546 0.083 -0.709 -0.546 -0.383 -0.546 46.521
## z1 0.841 0.068 0.708 0.841 0.975 0.841 91.432
##
## Deviance Information Criterion (DIC) ...............: 611.16
## Deviance Information Criterion (DIC, saturated) ....: 282.92
## Effective number of parameters .....................: 2.00
##
## Watanabe-Akaike information criterion (WAIC) ...: 610.80
## Effective number of parameters .................: 1.62
##
## Marginal log-Likelihood: -311.54
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
Misal model random intercept per grup (\(g=1,\ldots,G\)) untuk Gaussian: \[ y_{ig} = \beta_0 + \beta_1 x_{ig} + u_g + \varepsilon_{ig},\quad u_g\sim\mathcal{N}(0, \tau_u^{-1}),\ \varepsilon_{ig}\sim\mathcal{N}(0,\sigma^2). \]
Di INLA:
set.seed(7)
G <- 20
ng <- 20
x_h <- rnorm(G*ng)
grp <- rep(1:G, each = ng)
u <- rnorm(G, 0, 0.6)
y_h <- 1 + 1.2*x_h + u[grp] + rnorm(G*ng, 0, 1)
df_h <- data.frame(y_h, x_h, grp = as.factor(grp))
fit_h <- inla(
y_h ~ x_h + f(grp, model = "iid"),
data = df_h,
family = "gaussian",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE)
)
summary(fit_h)
## Time used:
## Pre = 0.982, Running = 0.192, Post = 0.0111, Total = 1.18
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 1.046 0.127 0.794 1.046 1.297 1.046 32.164
## x_h 1.205 0.049 1.109 1.205 1.301 1.205 255.030
##
## Random effects:
## Name Model
## grp IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 1.07 0.078 0.927 1.07
## Precision for grp 4.11 1.543 1.853 3.85
## 0.975quant mode
## Precision for the Gaussian observations 1.23 1.07
## Precision for grp 7.83 3.39
##
## Deviance Information Criterion (DIC) ...............: 1129.72
## Deviance Information Criterion (DIC, saturated) ....: 421.48
## Effective number of parameters .....................: 19.16
##
## Watanabe-Akaike information criterion (WAIC) ...: 1130.07
## Effective number of parameters .................: 18.66
##
## Marginal log-Likelihood: -601.15
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
|||| | Prinsip | Sampling posterior | Aproksimasi Laplace bertingkat | | Fleksibilitas | Sangat luas (dengan pemodelan yang tepat) | Luas untuk LGM; beberapa model khusus memerlukan trik | | Kecepatan | Lebih lambat, perlu konvergensi | Umumnya jauh lebih cepat | | Diagnostik | Perlu cek konvergensi (trace, R-hat, ESS) | Tidak perlu rantai; cek marginals & kriteria (DIC/WAIC/CPO) | | Output | Sampel penuh posterior | Marginal posterior aproksimasi |
Praktik baik: gunakan INLA untuk eksplorasi/produksi cepat pada model LGM standar; validasi silang dengan MCMC pada subset kecil ketika perlu.
Pada bagian ini kita akan melakukan simulasi data regresi linear sederhana dengan slope \(\beta_1\) yang ditentukan. Kita akan membandingkan:
set.seed(123)
n <- 20
beta0_true <- 1.0
beta1_true <- 2.0
sigma_true <- 1.2
x <- sort(runif(n, -2, 2))
X <- cbind(1, x)
y <- beta0_true + beta1_true * x + rnorm(n, 0, sigma_true)
ols_fit <- lm(y ~ x)
beta_ols <- coef(ols_fit)
beta_ols
## (Intercept) x
## 0.9323097 1.6955629
bayes_linreg_nig <- function(X, y, m0, V0, a0, b0) {
XtX <- crossprod(X)
Xty <- crossprod(X, y)
V0_inv <- tryCatch(solve(V0), error = function(e) matrix(0, nrow(V0), ncol(V0)))
Vn <- solve(V0_inv + XtX)
mn <- Vn %*% (V0_inv %*% m0 + Xty)
a_n <- a0 + nrow(X) / 2
b_n <- as.numeric(
b0 + 0.5 * ( crossprod(y, y) + t(m0) %*% V0_inv %*% m0 - t(mn) %*% solve(Vn) %*% mn )
)
list(mn = as.vector(mn), Vn = Vn, a_n = a_n, b_n = b_n)
}
beta1_target <- 0.5 # nilai slope target
tau_beta1 <- 0.2 # varian kecil -> prior kuat (bukan 'au_beta1')
# (opsional) guard agar tau_beta1 valid
stopifnot(is.numeric(tau_beta1), length(tau_beta1)==1, tau_beta1 > 0)
# Prior informatif
m0_inf <- c(0, beta1_target)
V0_inf <- diag(c(10^2, tau_beta1^2)) # pakai tau_beta1 yang benar
a0_inf <- 2; b0_inf <- 1
# Prior vague
m0_vag <- c(0, 0)
V0_vag <- diag(c(1e6, 1e6))
a0_vag <- 1e-3; b0_vag <- 1e-3
post_inf <- bayes_linreg_nig(X, y, m0_inf, V0_inf, a0_inf, b0_inf)
post_vag <- bayes_linreg_nig(X, y, m0_vag, V0_vag, a0_vag, b0_vag)
beta_bayes_inf <- post_inf$mn
beta_bayes_vag <- post_vag$mn
summary_tab <- data.frame(
model = c("True", "OLS", "Bayes (informative)", "Bayes (vague ≈ noninformative)"),
beta0 = c(beta0_true, beta_ols[1], beta_bayes_inf[1], beta_bayes_vag[1]),
beta1 = c(beta1_true, beta_ols[2], beta_bayes_inf[2], beta_bayes_vag[2]),
stringsAsFactors = FALSE
)
knitr::kable(summary_tab, digits = 3, caption = "Ringkasan estimasi beta0 dan beta1")
| model | beta0 | beta1 |
|---|---|---|
| True | 1.000 | 2.000 |
| OLS | 0.932 | 1.696 |
| Bayes (informative) | 1.042 | 1.151 |
| Bayes (vague ≈ noninformative) | 0.932 | 1.696 |
x_grid <- seq(min(x), max(x), length.out = 200)
y_true <- beta0_true + beta1_true * x_grid
y_ols <- beta_ols[1] + beta_ols[2] * x_grid
y_inf <- beta_bayes_inf[1]+ beta_bayes_inf[2]* x_grid
y_vag <- beta_bayes_vag[1]+ beta_bayes_vag[2]* x_grid
plot(x, y, pch = 19, cex = 1.1, xlab = "x", ylab = "y",
main = "OLS vs Bayesian Linear Regression")
lines(x_grid, y_true, lwd = 2, lty = 2, col = "black")
lines(x_grid, y_ols, lwd = 2, col = "blue")
lines(x_grid, y_inf, lwd = 2, lty = 3, col = "red")
lines(x_grid, y_vag, lwd = 2, lty = 4, col = "darkgreen")
legend("topleft",
legend = c(
sprintf("True: beta1=%.2f", beta1_true),
sprintf("OLS: beta1=%.2f", beta_ols[2]),
sprintf("Bayes (informative): beta1=%.2f (target=%.2f)", beta_bayes_inf[2], beta1_target),
sprintf("Bayes (vague): beta1=%.2f", beta_bayes_vag[2])
),
lty = c(2,1,3,4), lwd = 2, col = c("black","blue","red","darkgreen"), bty = "n")
Kita tinjau model regresi linear dengan notasi matriks: \[ \mathbf y \mid \boldsymbol\beta,\sigma^2 \sim \mathcal N\!\big(X\boldsymbol\beta,\ \sigma^2 I_n\big), \] dengan \(\mathbf y\in\mathbb R^n\), \(X\in\mathbb R^{n\times p}\), \(\boldsymbol\beta\in\mathbb R^p\), dan \(\sigma^2>0\).
Prior konjugat Normal–Inverse-Gamma (NIG): \[ \boldsymbol\beta\mid\sigma^2 \sim \mathcal N\!\big(m_0,\ \sigma^2 V_0\big),\qquad \sigma^2 \sim \mathrm{InvGamma}(a_0,b_0). \] Parameterisasi Inverse-Gamma yang dipakai di sini: \[ p(\sigma^2)=\frac{b_0^{a_0}}{\Gamma(a_0)}\,(\sigma^2)^{-(a_0+1)}\exp\!\Big(-\frac{b_0}{\sigma^2}\Big),\quad \sigma^2>0. \]
Likelihood (abaikan konstanta tak bergantung pada parameter): \[ p(\mathbf y\mid\boldsymbol\beta,\sigma^2)\ \propto\ (\sigma^2)^{-n/2}\exp\!\Big\{-\tfrac{1}{2\sigma^2}(\mathbf y-X\boldsymbol\beta)'(\mathbf y-X\boldsymbol\beta)\Big\}. \] Prior bersama (joint prior): \[ p(\boldsymbol\beta,\sigma^2)=p(\boldsymbol\beta\mid\sigma^2)\,p(\sigma^2) \ \propto\ (\sigma^2)^{-p/2}\exp\!\Big\{-\tfrac{1}{2\sigma^2}(\boldsymbol\beta-m_0)'V_0^{-1}(\boldsymbol\beta-m_0)\Big\}\cdot (\sigma^2)^{-(a_0+1)}\exp\!\Big(-\tfrac{b_0}{\sigma^2}\Big). \]
Gabungkan (kernel): \[ p(\boldsymbol\beta,\sigma^2\mid\mathbf y)\ \propto\ (\sigma^2)^{-\frac{n}{2}}\exp\!\Big\{-\tfrac{1}{2\sigma^2}\Vert\mathbf y-X\boldsymbol\beta\Vert^2\Big\}\cdot (\sigma^2)^{-\frac{p}{2}}\exp\!\Big\{-\tfrac{1}{2\sigma^2}(\boldsymbol\beta-m_0)'V_0^{-1}(\boldsymbol\beta-m_0)\Big\}\cdot (\sigma^2)^{-(a_0+1)}\exp\!\Big(-\tfrac{b_0}{\sigma^2}\Big). \]
Kumpulkan faktor (^2) dan selesaikan kuadrat untuk \(\boldsymbol\beta\). Definisikan \[ V_n=(V_0^{-1}+X'X)^{-1},\qquad m_n=V_n\big(V_0^{-1}m_0+X'\mathbf y\big). \] Gunakan identitas penyelesaian kuadrat: \[ \Vert\mathbf y-X\boldsymbol\beta\Vert^2+(\boldsymbol\beta-m_0)'V_0^{-1}(\boldsymbol\beta-m_0) = (\boldsymbol\beta-m_n)'V_n^{-1}(\boldsymbol\beta-m_n) + C, \] dengan konstanta (tidak bergantung pada \(\boldsymbol\beta\)) \[ C = \mathbf y'\mathbf y + m_0'V_0^{-1}m_0 - m_n'V_n^{-1}m_n. \]
Maka kernel posterior menjadi \[ p(\boldsymbol\beta,\sigma^2\mid\mathbf y)\ \propto\ (\sigma^2)^{-\frac{n+p}{2}}\,(\sigma^2)^{-(a_0+1)}\exp\!\Big\{-\tfrac{1}{2\sigma^2}(\boldsymbol\beta-m_n)'V_n^{-1}(\boldsymbol\beta-m_n)\Big\} \times \exp\!\Big\{-\tfrac{1}{2\sigma^2}C\Big\}\,\exp\!\Big(-\tfrac{b_0}{\sigma^2}\Big). \]
Kelompokkan lagi eksponen bagian \(\sigma^2\) (yang tidak mengandung \(\boldsymbol\beta\)): \[ \exp\!\Big\{-\tfrac{1}{\sigma^2}\Big(\tfrac{1}{2}C + b_0\Big)\Big\}. \]
Definisikan \[ a_n = a_0 + \tfrac{n}{2},\qquad b_n = b_0 + \tfrac{1}{2}\Big(\mathbf y'\mathbf y + m_0'V_0^{-1}m_0 - m_n'V_n^{-1}m_n\Big). \]
Dari bentuk kernel di atas, kita memperoleh faktorasi posterior konjugat: \[ \boxed{\;\boldsymbol\beta\mid\sigma^2,\mathbf y\ \sim\ \mathcal N\!\big(m_n,\ \sigma^2 V_n\big)\;}\qquad \boxed{\;\sigma^2\mid\mathbf y\ \sim\ \mathrm{InvGamma}(a_n,\ b_n)\;} \]
Dengan mengintegrasikan \(\sigma^2\), marginal \(\boldsymbol\beta\mid\mathbf y\) adalah sebaran Student-\(t\) multivariat: \[ \boldsymbol\beta\mid\mathbf y \sim \mathsf{MVT}_{\nu}\Big(m_n,\ \frac{b_n}{a_n}V_n\Big),\quad \nu=2a_n. \] Artinya, untuk vektor mana pun \(c\), \(c'\boldsymbol\beta\mid\mathbf y\) mengikuti Student-\(t_{\nu}\) dengan mean \(c'm_n\) dan varians \(\frac{b_n}{a_n}c'V_n c\).
Untuk kovariat baru \(x_*\in\mathbb R^p\) (termasuk 1 untuk intersep): \[ y_*\mid\mathbf y \sim t_{\nu}\Big(\mu_*,\ s_*^2\Big),\quad \mu_* = x_*'m_n,\ \ s_*^2 = \frac{b_n}{a_n}\big(1 + x_*'V_n x_*\big),\ \nu=2a_n. \]
set.seed(42)
n <- 40
p <- 2
beta_true <- c(1.0, 2.0) # (intercept, slope)
sigma_true <- 1.0
x <- sort(runif(n, -2, 2))
X <- cbind(1, x)
y <- as.vector(X %*% beta_true + rnorm(n, 0, sigma_true))
m0 <- c(0, 0)
V0 <- diag(c(10^2, 10^2))
a0 <- 2
b0 <- 1
XtX <- crossprod(X)
Xty <- crossprod(X, y)
V0_inv <- solve(V0)
Vn <- solve(V0_inv + XtX)
mn <- Vn %*% (V0_inv %*% m0 + Xty)
an <- a0 + n/2
bn <- as.numeric(b0 + 0.5 * ( crossprod(y,y) + t(m0) %*% V0_inv %*% m0 - t(mn) %*% solve(Vn) %*% mn ))
c(an = an, bn = bn)
## an bn
## 22.0000 22.9944
mn
## [,1]
## 0.849861
## x 2.038724
Vn
## x
## 0.027250856 -0.006255731
## x -0.006255731 0.017338219
Untuk beta|y: mean = m_n, kovarian bersyarat = E[sigma^2|y] * V_n =
(b_n/(a_n-1)) V_n (jika a_n>1).
Untuk sigma^2|y: mean = b_n/(a_n-1) (jika a_n>1).
mean_sigma2 <- bn / (an - 1) # a_n > 1
cov_beta <- as.numeric(mean_sigma2) * Vn
list(mean_sigma2 = mean_sigma2, cov_beta = cov_beta)
## $mean_sigma2
## [1] 1.094972
##
## $cov_beta
## x
## 0.029838914 -0.006849848
## x -0.006849848 0.018984857
coef(lm(y ~ x))
## (Intercept) x
## 0.849965 2.039024
x_grid <- seq(min(x), max(x), length.out = 100)
Xg <- cbind(1, x_grid)
mu_pred <- as.vector(Xg %*% mn)
scale2 <- as.numeric(bn/an) * (1 + rowSums((Xg %*% Vn) * Xg))
nu <- 2*an
# Interval kredibel prediktif 95% (Student-t)
qL <- mu_pred + sqrt(scale2) * qt(0.025, df = nu)
qU <- mu_pred + sqrt(scale2) * qt(0.975, df = nu)
plot(x, y, pch=19, xlab = "x", ylab = "y",
main = "Posterior Predictive t-interval (95%)")
lines(x_grid, mu_pred, lwd = 2)
lines(x_grid, qL, lty = 2)
lines(x_grid, qU, lty = 2)
Misalkan model populasi (berbasis-model / model-based):
\[ y_i = \mathbf{x}_i^\top \boldsymbol{\beta} + \varepsilon_i,\quad \mathbb{E}[\varepsilon_i \mid \mathbf{x}_i]=0,\quad \mathrm{Var}(\varepsilon_i\mid \mathbf{x}_i)=\sigma^2. \]
Dengan matriks desain \(\mathbf{X}\in\mathbb{R}^{n\times p}\) dan vektor respons \(\mathbf{y}\in\mathbb{R}^n\), penaksir OLS:
\[ \hat{\boldsymbol{\beta}}_{\text{OLS}}=(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}, \qquad \widehat{\mathrm{Var}}(\hat{\boldsymbol{\beta}}_{\text{OLS}})=\hat{\sigma}^2(\mathbf{X}^\top \mathbf{X})^{-1}, \]
dengan \(\hat{\sigma}^2=\frac{1}{n-p}\sum_{i=1}^{n}\hat{\varepsilon}_i^2\).
Regresi Berbobot (WLS) & Bobot Sampling
Jika setiap unit \(i\) memiliki probabilitas inklusi \(\pi_i\) dan bobot dasar \(w_i = 1/\pi_i\), maka pendekatan WLS (sering dipadankan dengan penimbang invers-probabilitas) adalah:
\[ \hat{\boldsymbol{\beta}}_{\text{WLS}}=(\mathbf{X}^\top \mathbf{W}\mathbf{X})^{-1}\mathbf{X}^\top \mathbf{W} \mathbf{y}, \]
dengan \(\mathbf{W}=\mathrm{diag}(w_1,\dots,w_n)\). Varians model-based untuk WLS (di bawah homoskedastisitas residu terukur setelah pembobotan) adalah:
\[ \widehat{\mathrm{Var}}(\hat{\boldsymbol{\beta}}_{\text{WLS}})=\hat{\sigma}_w^2(\mathbf{X}^\top \mathbf{W}\mathbf{X})^{-1}, \]
namun untuk data survei, varians yang tepat biasanya menggunakan sandwich/robust atau prosedur berbasis-desain (Taylor linearization, replicate weights).
Pada banyak studi survei, kita memiliki populasi berhingga berukuran \(N\) dan mengambil sampel ukuran \(n\) dengan desain tertentu (SRS, stratifikasi, klaster, PPS, dsb.). Ketika melakukan regresi untuk mengestimasi hubungan \(y\) terhadap kovariat \(X\), hasil dan ketidakpastian estimasi dapat berbeda jika kita mengabaikan desain sampling yang tidak-acak sederhana. Dalam regresi linear klasik, kita biasanya mengasumsikan bahwa data yang digunakan adalah sampel acak sederhana (SRS, simple random sample) dari populasi. Artinya, setiap unit dalam populasi memiliki peluang yang sama untuk terpilih, dan pengamatan diasumsikan saling independen.
Namun, pada praktik survei atau penelitian lapangan, desain sampling sering kali lebih kompleks:
Akibatnya, peluang inklusi berbeda-beda antar unit (\(\pi_i \neq \pi_j\)), dan observasi dalam satu klaster tidak independen (ada korelasi intraklaster).
Varians rataan sampel: \[ \mathrm{Var}(\bar{y})=\left(1-\frac{n}{N}\right)\frac{S_y^2}{n},\quad S_y^2=\frac{1}{N-1}\sum_{i=1}^{N}(y_i-\bar{Y})^2. \]
Menurunkan dan mengilustrasikan rumus ukuran sampel \(n\) untuk mengestimasi rata-rata populasi \(\mu\) dengan SRS tanpa pengembalian pada target margin of error \(d\) di tingkat keyakinan \(1-\alpha\). Rumus yang akan kita capai (menggunakan varian populasi berhingga dengan penyebut \(N-1\)) adalah:
\[ n \;\approx\; \frac{z_{\alpha/2}^2\, S_y^2}{d^2}\;\cdot\;\frac{N}{\,N-1 \;+\; z_{\alpha/2}^2 S_y^2 / d^2\,}. \]
Setting & Asumsi
Definisi Margin of Error (MoE) Untuk interval keyakinan
\(1-\alpha\):
\[
\bar{y} \pm
z_{\alpha/2}\,\sqrt{\left(1-\frac{n}{N}\right)\frac{S_y^2}{n}}.
\] Target margin of error \(d\) adalah:
\[
d \;=\; z_{\alpha/2}\,\sqrt{\left(1-\frac{n}{N}\right)\frac{S_y^2}{n}}.
\]
Penurunan Rumus \(n\) (Langkah demi langkah)
Mulai dari persamaan margin of error, kuadratkan kedua sisi: \[ d^2 \;=\; z_{\alpha/2}^2 \left(1-\frac{n}{N}\right)\frac{S_y^2}{n} = z_{\alpha/2}^2\,S_y^2 \cdot \frac{N-n}{n\,N}. \]
Kalikan silang: \[ d^2\, n\, N \;=\; z_{\alpha/2}^2\,S_y^2 \,(N-n) \;=\; z_{\alpha/2}^2\,S_y^2\,N \;-\; z_{\alpha/2}^2\,S_y^2\,n. \]
Kumpulkan suku yang mengandung \(n\) di sebelah kiri: \[ d^2 n N \;+\; z_{\alpha/2}^2 S_y^2 n \;=\; z_{\alpha/2}^2 S_y^2 N. \]
Faktorkan \(n\): \[ n\left( d^2 N \;+\; z_{\alpha/2}^2 S_y^2 \right) \;=\; z_{\alpha/2}^2 S_y^2 N. \]
Bagi kedua sisi: \[ \boxed{\; n \;=\; \frac{z_{\alpha/2}^2 S_y^2 N}{d^2 N + z_{\alpha/2}^2 S_y^2}\; }. \]
Dengan definisi varian populasi berhingga \(S_y^2 = \tfrac{1}{N-1}\sum (y_i-\mu)^2\), bentuk yang sering dipakai (dan yang terlihat pada slide) adalah bentuk aljabar ekuivalen: \[ \boxed{\; n \;\approx\; \frac{z_{\alpha/2}^2\, S_y^2}{d^2}\cdot\frac{N}{\,N-1 + z_{\alpha/2}^2 S_y^2 / d^2\,}\; }. \]
Catatan: Kedua bentuk identik secara praktis; perbedaan \(N\) vs \(N-1\) berasal dari konvensi definisi \(S_y^2\) pada populasi berhingga.
Bentuk Alternatif yang Sering Dipakai
Sering dikenalkan dua langkah: 1) Abaikan FPC (anggap \(N\) sangat besar): \[ n_0 = \frac{z_{\alpha/2}^2\,S_y^2}{d^2}. \] 2) Koreksi FPC (untuk SRS tanpa pengembalian): \[ n \;=\; \frac{n_0}{\,1 + \frac{n_0-1}{N}\,}. \]
Ilustrasi Numerik Misal:
N <- 10000
Sy <- 15
d <- 1.5
z <- 1.96
n0 <- (z^2 * Sy^2) / d^2 # tanpa FPC (anggap N sangat besar)
nF <- n0 / (1 + (n0 - 1)/N) # dengan FPC (SRS w/o replacement)
data.frame(n0_tanpa_FPC = ceiling(n0), n_dengan_FPC = ceiling(nF))
## n0_tanpa_FPC n_dengan_FPC
## 1 385 370
Interpretasi: \(n_0\) adalah kebutuhan sampel jika populasi dianggap sangat besar (tanpa FPC). Karena \(N\) terbatas (10.000), kita pakai \(n\) dengan FPC (biasanya lebih kecil dari \(n_0\)).
Validasi Singkat via Simulasi Kita buat populasi sintetis berukuran \(N\) dengan ragam mendekati \(S_y^2\), ambil SRS berukuran \(n\) (hasil dari perhitungan dengan FPC), lalu cek margin of error empiris kira-kira mendekati target \(d\).
set.seed(123)
N <- 10000
Sy <- 15
mu <- 50 # mean populasi (arbitrer)
d <- 1.5
z <- 1.96
# Populasi sintetis ~ Normal(mu, Sy^2)
pop <- rnorm(N, mean = mu, sd = Sy)
# Hitung n dengan FPC seperti di atas
n0 <- (z^2 * Sy^2) / d^2
n <- ceiling(n0 / (1 + (n0 - 1)/N))
# Uji simulasi: ambil banyak SRS, ukur half-width CI
B <- 1000
hw <- numeric(B) # half-width (MoE) empiris
for (b in 1:B) {
samp <- sample(pop, n, replace = FALSE)
ybar <- mean(samp)
# SE SRS w/o replacement menggunakan ragam sampel (unbiased) * FPC
s2 <- var(samp) # var sampel (penyebut n-1)
se <- sqrt( (1 - n/N) * s2 / n ) # SE design-based
hw[b] <- z * se # half-width (MoE)
}
c( target_d = d, mean_empirical_MoE = mean(hw), sd_MoE = sd(hw), n_used = n )
## target_d mean_empirical_MoE sd_MoE n_used
## 1.50000000 1.49867515 0.05475821 370.00000000
Biasanya nilai rata-rata MoE empiris akan mendekati target \(d\), menunjukkan bahwa rumus ukuran sampel bekerja dengan baik.
Catatan
Ringkasan Rumus ukuran sampel SRS untuk mean (dengan FPC) pada target margin of error \(d\) dan keyakinan \(1-\alpha\): \[ \boxed{\; n \;\approx\; \frac{z_{\alpha/2}^2\, S_y^2}{d^2}\cdot\frac{N}{\,N-1 + z_{\alpha/2}^2 S_y^2 / d^2\,} \;=\; \frac{n_0}{1 + \frac{n_0-1}{N}} \;} \] dengan \(n_0 = z_{\alpha/2}^2 S_y^2 / d^2\).
Jika kita langsung menerapkan OLS biasa (unweighted regression) tanpa memperhitungkan bobot atau struktur desain:
Bias pada koefisien \(\hat{\beta}\).
Jika probabilitas inklusi berkorelasi dengan variabel respons (\(y\)) atau kovariat (\(\mathbf{X}\)), maka rata-rata sampel tidak
lagi merepresentasikan rata-rata populasi.
Contoh: Jika rumah tangga kaya lebih sering terpilih dalam sampel,
regresi pengaruh pendidikan terhadap pendapatan bisa
overestimate.
Underestimation of standard errors.
Pada desain klaster, unit dalam klaster lebih mirip satu sama lain \(\Rightarrow\) effective sample
size lebih kecil daripada \(n\).
Jika kita abaikan, maka standar error dari koefisien regresi akan
terlalu kecil, menyebabkan uji signifikansi terlalu
optimis (type I error naik).
Salah interpretasi inferensi.
Regresi OLS biasa hanya valid untuk sampel yang diobservasi, bukan untuk
inferensi ke populasi target. Desain sampling menentukan “who
represents whom” di populasi, sehingga jika diabaikan, hasil
regresi bisa gagal digeneralisasi.
Untuk mengatasi masalah ini, digunakan dua pendekatan:
Terapkan bobot sampling \(w_i =
1/\pi_i\) pada regresi: \[
\hat{\boldsymbol{\beta}}_{\text{WLS}} \;=\; (\mathbf{X}^\top \mathbf{W}
\mathbf{X})^{-1} \mathbf{X}^\top \mathbf{W} \mathbf{y},
\] dengan \(\mathbf{W}=\mathrm{diag}(w_1,\dots,w_n)\).
Pendekatan ini menjamin bahwa unit dengan probabilitas inklusi rendah
tetap “terwakili”.
Konsep Bobot W Dalam sampling survei, setiap unit populasi (misalnya rumah tangga, individu) punya peluang untuk masuk ke dalam sampel.
Bobot Sampling Dasar (\(w_i\))
Untuk membuat sampel representatif terhadap populasi, digunakan inverse probability weighting (IPW): \[ w_i = \frac{1}{\pi_i}. \]
Mengapa Probabilitas Besar → Bobot Kecil?
Misalkan populasi ada 1000 orang, tapi kita hanya pilih 100 (10%).
Jika probabilitas inklusi lebih besar, misalnya \(\pi_i = 0.50\) (setiap orang punya peluang 50% dipilih), maka:
Jadi: Semakin mudah sebuah unit masuk sampel, semakin sedikit “berat” yang harus ia bawa dalam merepresentasikan populasi.
WLS (Weighted Least Squares) dalam Konteks Survei
Dengan bobot \(w_i\), regresi dilakukan dengan meminimalkan fungsi kuadrat berbobot: \[ \hat\beta_{WLS} = \arg\min_\beta \sum_{i \in s} w_i \,(y_i - x_i^\top \beta)^2. \]
Rumus penutupnya: \[ \hat\beta_{WLS} = (X^\top W X)^{-1} X^\top W y, \] dengan \(W = \mathrm{diag}(w_1,\ldots,w_n)\).
Interpretasi: bobot memastikan bahwa koefisien regresi memperhitungkan peluang pemilihan yang berbeda antar unit.
Ringkasan
Ketika kita menggunakan desain sampling kompleks (stratifikasi, klaster, PPS), menghitung varians estimasi (terutama pada regresi) tidak bisa langsung memakai rumus OLS biasa. Ada beberapa pendekatan umum yang dipakai dalam analisis data survei:
Rumus
Misalkan estimator dapat ditulis sebagai fungsi \(\hat{\theta} = f(\bar{y}, \mathbf{x})\). Dengan ekspansi Taylor orde pertama di sekitar nilai harapan, diperoleh: \[ \hat{\theta} \approx \theta \;+\; \frac{\partial f}{\partial \bar{y}}\left(\bar{y} - \mathbb{E}[\bar{y}]\right) \;+\; \cdots \] Sehingga variansnya (aproksimasi) dapat dipropagasikan sebagai: \[ \mathrm{Var}(\hat{\theta}) \;\approx\; \left(\frac{\partial f}{\partial \bar{y}}\right)^2 \mathrm{Var}(\bar{y}). \]
Contoh
svyglm() umumnya dihitung dengan
pendekatan linearization.Alih-alih menghitung varians dengan rumus analitik, kita mengulang estimasi pada beberapa replicate samples yang dibentuk dari sampel asli.
Metode umum - Jackknife Repeated Replication (JRR / JK1, JK2): buang satu unit (atau satu klaster) dari sampel, hitung kembali estimator, lalu lihat variasinya. - Balanced Repeated Replication (BRR): bagi sampel ke dalam half-samples (setengah data), lakukan estimasi berkali-kali dengan kombinasi setengah sampel; populer di survei besar seperti CPS. - Bootstrap untuk survei: ambil ulang (resampling) unit/klaster dengan skema yang meniru desain sampling, lalu hitung variasinya.
Rumus Sketsa Jika terdapat \(R\) replikasi, estimator pada replikasi ke-\(r\) adalah \(\hat{\theta}^{(r)}\). Varians replikasi dihitung sebagai: \[ \widehat{\mathrm{Var}}(\hat{\theta}) \;=\; \frac{1}{R}\sum_{r=1}^{R}\left(\hat{\theta}^{(r)} - \bar{\theta}\right)^2, \qquad \bar{\theta} \;=\; \frac{1}{R}\sum_{r=1}^{R} \hat{\theta}^{(r)}. \]
Kelebihan
Intuisi Jika populasi berhingga berukuran \(N\) dan kita ambil sampel tanpa pengembalian sebesar \(n\), maka ketidakpastian lebih kecil dibanding sampel dengan pengembalian, karena setiap unit yang sudah terpilih tidak bisa muncul lagi.
Rumus Untuk mean (atau estimator lain berbasis rata-rata) pada SRS tanpa pengembalian: \[ \mathrm{Var}(\bar{y}) \;=\; \left(1 - \frac{n}{N}\right)\frac{S_y^2}{n}, \] dengan \(S_y^2\) ragam populasi.
Ringkasnya, FPC sebagai pengali standard error adalah: \[ \mathrm{FPC} \;=\; \sqrt{\frac{N - n}{N - 1}}. \]
Contoh Numerik
Jika \(N=10{,}000\) dan \(n=2{,}000\), maka: \[ \mathrm{FPC} \;\approx\; \sqrt{\frac{10{,}000 - 2{,}000}{10{,}000 - 1}} \;\approx\; 0.894. \] Artinya, standard error berkurang ~10% dibanding asumsi populasi tak berhingga.
Jika \(n \ll N\) (misalnya 1% populasi), maka FPC \(\approx 1\) (efeknya tidak signifikan).
Ringkasan
Implementasi praktis di R: svyglm() dari paket
survey (Lumley, 2010).
Contoh Ilustrasi
Catatan
svyglm) digunakan,
tergantung tujuan analisis.Untuk Koefisien Regresi (Pendekatan Aproksimasi)
Jika varian kesalahan \(\sigma^2\) dan ragam kolom \(\mathbf{x}_j\) diketahui (atau diperkirakan), maka (OLS): \[ \mathrm{Var}(\hat{\beta}_j) \approx \sigma^2 \left[(\mathbf{X}^\top \mathbf{X})^{-1}\right]_{jj}. \] Aproksimasi kuat: pastikan \(n \gg p\), \(\mathrm{DEFF}\) mengembang varian efektif \(\Rightarrow n_{\text{efektif}} \approx n/\mathrm{DEFF}\).
Catatan Regresi Mengabaikan vs Memperhitungkan Desain
survey::svyglm atau
srvyr.Data Populasi Sintetik
Kita bangun populasi berhingga untuk simulasi: variabel prediktor \(x_1, x_2\), sebuah variabel strata, dan klaster.
N <- 10000
# strata: 3 strata berukuran tidak sama
strata <- sample(letters[1:3], size = N, replace = TRUE, prob = c(0.5, 0.3, 0.2))
cluster <- sample(1:500, size = N, replace = TRUE) # 500 klaster
x1 <- rnorm(N, mean = ifelse(strata=="a", 0, ifelse(strata=="b", 1, -1)), sd = 1.2)
x2 <- runif(N, -2, 2)
# Probabilitas seleksi (untuk PPS gaya sederhana): unit dg x1 besar sedikit lebih mungkin terpilih
p_select <- plogis( -1 + 0.2*x1 )
y <- 1.5 + 0.8*x1 - 0.5*x2 + rnorm(N, 0, 1)
pop <- data.frame(id = 1:N, strata, cluster, x1, x2, y, p_select)
head(pop)
## id strata cluster x1 x2 y p_select
## 1 1 c 424 1.54823092 1.9386524 2.1082381 0.3339544
## 2 2 a 18 0.05000342 -0.5331374 1.4897788 0.2709122
## 3 3 a 152 -0.32109098 1.6701588 0.1394602 0.2565042
## 4 4 a 261 -0.99588549 0.9906781 1.6725408 0.2316216
## 5 5 b 428 -1.29522503 -0.1892852 1.2128602 0.2211383
## 6 6 a 313 -0.26585865 0.7055983 -0.9148585 0.2586165
Ambil Sampel: (a) SRS, (b) Stratified, (c) PPS sederhana
n <- 800
# (a) SRS
srs_idx <- sample(pop$id, size = n, replace = FALSE)
srs <- pop[srs_idx, ]
srs$pi <- n / N
srs$w <- 1 / srs$pi # = N/n
# (b) Stratified proportional allocation
n_h <- round(prop.table(table(pop$strata)) * n)
strat_idx <- unlist(lapply(names(n_h), function(h) {
ids <- which(pop$strata == h)
sample(ids, size = n_h[h], replace = FALSE)
}))
strat <- pop[strat_idx, ]
Nh <- as.numeric(table(pop$strata))[match(strat$strata, names(table(pop$strata)))]
nh <- as.numeric(table(strat$strata))[match(strat$strata, names(table(strat$strata)))]
strat$pi <- nh / Nh
strat$w <- 1 / strat$pi
# (c) PPS-like (berbasis p_select sebagai ukuran ukuran)
pps_idx <- sample(pop$id, size = n, replace = FALSE, prob = pop$p_select / sum(pop$p_select))
pps <- pop[pps_idx, ]
pps$pi <- (pop$p_select[pps_idx] / sum(pop$p_select)) * n # aproksimasi
pps$pi <- pmin(pmax(pps$pi, 1e-6), 1) # jaga batas
pps$w <- 1 / pps$pi
Regresi OLS vs WLS (abaikan vs perhatikan bobot)
m_srs_ols <- lm(y ~ x1 + x2, data = srs)
m_srs_wls <- lm(y ~ x1 + x2, data = srs, weights = w)
m_str_ols <- lm(y ~ x1 + x2, data = strat)
m_str_wls <- lm(y ~ x1 + x2, data = strat, weights = w)
m_pps_ols <- lm(y ~ x1 + x2, data = pps)
m_pps_wls <- lm(y ~ x1 + x2, data = pps, weights = w)
summary(m_srs_ols); summary(m_srs_wls)
##
## Call:
## lm(formula = y ~ x1 + x2, data = srs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0961 -0.6639 0.0001 0.6450 3.4191
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.50965 0.03561 42.40 <2e-16 ***
## x1 0.81270 0.02561 31.73 <2e-16 ***
## x2 -0.56941 0.03075 -18.52 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.003 on 797 degrees of freedom
## Multiple R-squared: 0.6261, Adjusted R-squared: 0.6252
## F-statistic: 667.4 on 2 and 797 DF, p-value: < 2.2e-16
##
## Call:
## lm(formula = y ~ x1 + x2, data = srs, weights = w)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.9465 -2.3471 0.0003 2.2805 12.0883
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.50965 0.03561 42.40 <2e-16 ***
## x1 0.81270 0.02561 31.73 <2e-16 ***
## x2 -0.56941 0.03075 -18.52 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.547 on 797 degrees of freedom
## Multiple R-squared: 0.6261, Adjusted R-squared: 0.6252
## F-statistic: 667.4 on 2 and 797 DF, p-value: < 2.2e-16
summary(m_str_ols); summary(m_str_wls)
##
## Call:
## lm(formula = y ~ x1 + x2, data = strat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1292 -0.7031 -0.0449 0.6900 4.3023
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.46118 0.03706 39.43 <2e-16 ***
## x1 0.81544 0.02676 30.47 <2e-16 ***
## x2 -0.44771 0.03221 -13.90 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.045 on 797 degrees of freedom
## Multiple R-squared: 0.5837, Adjusted R-squared: 0.5827
## F-statistic: 558.8 on 2 and 797 DF, p-value: < 2.2e-16
##
## Call:
## lm(formula = y ~ x1 + x2, data = strat, weights = w)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -11.0648 -2.4863 -0.1587 2.4392 15.1920
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.46115 0.03706 39.43 <2e-16 ***
## x1 0.81544 0.02676 30.47 <2e-16 ***
## x2 -0.44772 0.03221 -13.90 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.695 on 797 degrees of freedom
## Multiple R-squared: 0.5837, Adjusted R-squared: 0.5827
## F-statistic: 558.7 on 2 and 797 DF, p-value: < 2.2e-16
summary(m_pps_ols); summary(m_pps_wls)
##
## Call:
## lm(formula = y ~ x1 + x2, data = pps)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9285 -0.7086 -0.0069 0.6715 4.2551
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.54669 0.03614 42.79 <2e-16 ***
## x1 0.81105 0.02469 32.85 <2e-16 ***
## x2 -0.47779 0.03020 -15.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.002 on 797 degrees of freedom
## Multiple R-squared: 0.6284, Adjusted R-squared: 0.6275
## F-statistic: 674 on 2 and 797 DF, p-value: < 2.2e-16
##
## Call:
## lm(formula = y ~ x1 + x2, data = pps, weights = w)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -11.3721 -2.4879 -0.0284 2.3852 18.3571
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.54556 0.03529 43.80 <2e-16 ***
## x1 0.81829 0.02451 33.38 <2e-16 ***
## x2 -0.46600 0.03005 -15.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.554 on 797 degrees of freedom
## Multiple R-squared: 0.6321, Adjusted R-squared: 0.6312
## F-statistic: 684.8 on 2 and 797 DF, p-value: < 2.2e-16
Catatan: Pada desain PPS (atau desain yang probabilitas inklusinya terkait dengan kovariat), OLS dapat memberikan estimasi koefisien yang bias dibanding WLS berbobot invers-probabilitas.
Varian Berbasis-Desain dengan
survey
Untuk memperoleh standar error yang benar
(memperhitungkan stratifikasi/klaster/FPC), gunakan paket
survey.
# install.packages("survey") # jika belum ada
library(survey)
## Loading required package: grid
##
## Attaching package: 'survey'
## The following object is masked from 'package:graphics':
##
## dotchart
# SRS tanpa klaster/strata (dapat memberi fpc=N untuk koreksi populasi berhingga)
des_srs <- svydesign(ids = ~1, weights = ~w, data = srs, fpc = ~I(rep(N, nrow(srs))))
# Stratified (tanpa klaster, hanya strata)
# Perlu supply ukuran populasi strata (opsional, untuk FPC). Kita beri approx via Nh_map.
Nh_map <- as.data.frame(table(pop$strata)); names(Nh_map) <- c("strata","Nh")
strat2 <- merge(strat, Nh_map, by = "strata")
des_str <- svydesign(ids = ~1, strata = ~strata, weights = ~w, data = strat2, fpc = ~Nh)
# Clustered (misal gunakan srs tapi anggap cluster benar-benar klaster sampling)
# Di sini contoh sintetik: treat 'cluster' dalam srs sebagai klaster (hanya contoh)
des_cl <- svydesign(ids = ~cluster, weights = ~w, data = srs)
# Regresi linear design-based
fit_srs <- svyglm(y ~ x1 + x2, design = des_srs)
fit_str <- svyglm(y ~ x1 + x2, design = des_str)
fit_cl <- svyglm(y ~ x1 + x2, design = des_cl)
summary(fit_srs)
##
## Call:
## svyglm(formula = y ~ x1 + x2, design = des_srs)
##
## Survey design:
## svydesign(ids = ~1, weights = ~w, data = srs, fpc = ~I(rep(N,
## nrow(srs))))
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.50965 0.03397 44.44 <2e-16 ***
## x1 0.81270 0.02298 35.36 <2e-16 ***
## x2 -0.56941 0.02936 -19.40 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.004002)
##
## Number of Fisher Scoring iterations: 2
summary(fit_str)
##
## Call:
## svyglm(formula = y ~ x1 + x2, design = des_str)
##
## Survey design:
## svydesign(ids = ~1, strata = ~strata, weights = ~w, data = strat2,
## fpc = ~Nh)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.46115 0.03565 40.99 <2e-16 ***
## x1 0.81544 0.02819 28.93 <2e-16 ***
## x2 -0.44772 0.03125 -14.32 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.089355)
##
## Number of Fisher Scoring iterations: 2
summary(fit_cl)
##
## Call:
## svyglm(formula = y ~ x1 + x2, design = des_cl)
##
## Survey design:
## svydesign(ids = ~cluster, weights = ~w, data = srs)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.50965 0.03361 44.92 <2e-16 ***
## x1 0.81270 0.02384 34.09 <2e-16 ***
## x2 -0.56941 0.03115 -18.28 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.004002)
##
## Number of Fisher Scoring iterations: 2
Interpretasi: Standar error dari svyglm
bisa lebih besar daripada OLS biasa, khususnya saat ada
klaster (intra-cluster correlation) atau
stratifikasi dengan design effect.
Perbandingan Ringkas Hasil
to_row <- function(fit, label) {
co <- coef(summary(fit))
data.frame(
model = label,
term = rownames(co),
estimate = co[, "Estimate"],
se = co[, "Std. Error"],
t_value = co[, "t value"],
p_value = co[, "Pr(>|t|)"]
)
}
tab <- rbind(
to_row(m_srs_ols, "OLS (SRS, unweighted)"),
to_row(m_srs_wls, "WLS (SRS, weights)"),
to_row(m_str_ols, "OLS (Strata, unweighted)"),
to_row(m_str_wls, "WLS (Strata, weights)"),
to_row(m_pps_ols, "OLS (PPS, unweighted)"),
to_row(m_pps_wls, "WLS (PPS, weights)"),
to_row(fit_srs, "Design-based (svyglm, SRS)"),
to_row(fit_str, "Design-based (svyglm, Strata)"),
to_row(fit_cl, "Design-based (svyglm, Cluster)")
)
tab
## model term estimate se
## (Intercept) OLS (SRS, unweighted) (Intercept) 1.5096475 0.03560627
## x1 OLS (SRS, unweighted) x1 0.8126989 0.02561251
## x2 OLS (SRS, unweighted) x2 -0.5694111 0.03074522
## (Intercept)1 WLS (SRS, weights) (Intercept) 1.5096475 0.03560627
## x11 WLS (SRS, weights) x1 0.8126989 0.02561251
## x21 WLS (SRS, weights) x2 -0.5694111 0.03074522
## (Intercept)2 OLS (Strata, unweighted) (Intercept) 1.4611849 0.03705895
## x12 OLS (Strata, unweighted) x1 0.8154373 0.02675833
## x22 OLS (Strata, unweighted) x2 -0.4477122 0.03220716
## (Intercept)3 WLS (Strata, weights) (Intercept) 1.4611507 0.03705957
## x13 WLS (Strata, weights) x1 0.8154374 0.02675979
## x23 WLS (Strata, weights) x2 -0.4477239 0.03220649
## (Intercept)4 OLS (PPS, unweighted) (Intercept) 1.5466933 0.03614242
## x14 OLS (PPS, unweighted) x1 0.8110501 0.02469273
## x24 OLS (PPS, unweighted) x2 -0.4777939 0.03020015
## (Intercept)5 WLS (PPS, weights) (Intercept) 1.5455601 0.03528804
## x15 WLS (PPS, weights) x1 0.8182904 0.02451364
## x25 WLS (PPS, weights) x2 -0.4659998 0.03004672
## (Intercept)6 Design-based (svyglm, SRS) (Intercept) 1.5096475 0.03397192
## x16 Design-based (svyglm, SRS) x1 0.8126989 0.02298328
## x26 Design-based (svyglm, SRS) x2 -0.5694111 0.02935609
## (Intercept)7 Design-based (svyglm, Strata) (Intercept) 1.4611507 0.03564741
## x17 Design-based (svyglm, Strata) x1 0.8154374 0.02818551
## x27 Design-based (svyglm, Strata) x2 -0.4477239 0.03125375
## (Intercept)8 Design-based (svyglm, Cluster) (Intercept) 1.5096475 0.03360921
## x18 Design-based (svyglm, Cluster) x1 0.8126989 0.02384313
## x28 Design-based (svyglm, Cluster) x2 -0.5694111 0.03114928
## t_value p_value
## (Intercept) 42.39837 1.791341e-206
## x1 31.73055 1.639441e-143
## x2 -18.52032 5.819122e-64
## (Intercept)1 42.39837 1.791341e-206
## x11 31.73055 1.639441e-143
## x21 -18.52032 5.819122e-64
## (Intercept)2 39.42867 1.914309e-189
## x12 30.47415 7.731180e-136
## x22 -13.90101 1.708304e-39
## (Intercept)3 39.42708 1.955369e-189
## x13 30.47249 7.914653e-136
## x23 -13.90167 1.695746e-39
## (Intercept)4 42.79440 1.037486e-208
## x14 32.84570 2.712585e-150
## x24 -15.82091 3.110001e-49
## (Intercept)5 43.79841 2.404339e-214
## x15 33.38102 1.549238e-153
## x25 -15.50917 1.320807e-47
## (Intercept)6 44.43810 6.583447e-218
## x16 35.36044 1.896689e-165
## x26 -19.39670 6.005869e-69
## (Intercept)7 40.98898 3.010079e-198
## x17 28.93109 2.722337e-126
## x27 -14.32545 1.418085e-41
## (Intercept)8 44.91767 1.014656e-157
## x18 34.08524 5.708805e-120
## x28 -18.28007 1.331520e-54
Memberikan rumus ukuran sampel untuk analisis regresi linear secara umum berbasis power analysis, termasuk:
pwr::pwr.f2.test) serta
simulasi singkat untuk validasi.Model regresi linear: \[
\mathbf{y} = \mathbf{X}\beta + \varepsilon,\qquad \varepsilon \sim
\mathcal{N}(0, \sigma^2 \mathbf{I}).
\] Estimator OLS: \(\hat{\beta} =
(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top
\mathbf{y}\).
Kovarians koefisien: \[
\mathrm{Var}(\hat{\beta}) = \sigma^2(\mathbf{X}^\top \mathbf{X})^{-1}.
\] Untuk koefisien ke-\(j\):
\(\mathrm{Var}(\hat{\beta}_j) =
\sigma^2[(\mathbf{X}^\top \mathbf{X})^{-1}]_{jj}\).
Uji hipotesis dua sisi: \(H_0:\beta_j=0\) vs \(H_1:\beta_j\neq 0\) pada taraf \(\alpha\) dan power \(1-\beta\).
Kriteria (pendekatan normal large-sample): \[
\frac{|\beta_j|}{\mathrm{SE}(\hat{\beta}_j)} \;\gtrsim\; z_{1-\alpha/2}
+ z_{1-\beta}.
\] Ini menuntun ke kebutuhan ukuran sampel
melalui \(\mathrm{SE}(\hat{\beta}_j)\)
yang menurun saat \(n\) naik.
Ukuran efek parsial untuk prediktor atau set
prediktor adalah: \[
f^2 \;=\; \frac{R^2_{\text{partial}}}{1 - R^2_{\text{partial}}}.
\] Konvensi Cohen (1988): kecil \(0.02\), sedang \(0.15\), besar \(0.35\).
Untuk uji set \(q\)
prediktor (mis. tambahan blok variabel), gunakan \(u=q\) derajat bebas uji numerik.
Dengan \(p\) prediktor total (termasuk yang diuji dan pengendali), uji \(u\) prediktor (biasanya \(u=1\) untuk satu koefisien), pendekatan praktis: \[ \boxed{ \; n \;\gtrsim\; \frac{(z_{1-\alpha/2} + z_{1-\beta})^2}{f^2} \;+\; p \;+\; 1 \; } \] Penjelasan ringkas: derivasi muncul dari pendekatan uji-F untuk peningkatan \(R^2\) dengan noncentrality parameter \(\lambda \approx n f^2\) dan pendekatan normal untuk kuantil F.
Catatan: rumus di atas aproksimasi; untuk desain kecil/eksak gunakan fungsi daya (
pwr.f2.test) atau uji-F non-sentral.
Jika data diperoleh via desain kompleks, skala \(n\) dengan DEFF (design
effect): \[
n_{\text{desain}} \;\approx\; \text{DEFF} \times n_{\text{SRS}}.
\] Contoh: jika \(n_{\text{SRS}}=100\) dan \(\text{DEFF}=2\) (klaster kuat), maka \(n_{\text{desain}}\approx 200\).
Untuk populasi berhingga dan SRS tanpa pengembalian, dapat pula
menerapkan FPC pada varian, tetapi pada tahap
perencanaan power biasanya cukup melalui DEFF.
Contoh Numerik Target: \(\alpha=0.05\), power \(1-\beta=0.80\), efek parsial sedang \(f^2=0.15\), jumlah prediktor total \(p=5\) (uji satu prediktor, \(u=1\)).
alpha <- 0.05; power <- 0.80; f2 <- 0.15; p <- 5
z_alpha <- qnorm(1 - alpha/2)
z_beta <- qnorm(power)
n_srs <- ( (z_alpha + z_beta)^2 / f2 ) + p + 1
ceiling(n_srs)
## [1] 59
DEFF <- 2.0
n_design <- ceiling(DEFF * n_srs)
data.frame(n_SRS = ceiling(n_srs), DEFF = DEFF, n_design = n_design)
## n_SRS DEFF n_design
## 1 59 2 117
pwrAlternatif yang lebih presisi untuk uji set prediktor menggunakan Cohen’s \(f^2\):
# install.packages("pwr") # jika perlu
suppressWarnings(suppressMessages(library(pwr)))
res <- pwr.f2.test(u = 1, v = NULL, f2 = f2, sig.level = alpha, power = power)
# v = n - p - 1 -> rearrange to n
n_from_pwr <- ceiling(res$v + p + 1)
list(pwr_result = res, n_from_pwr = n_from_pwr)
## $pwr_result
##
## Multiple regression power calculation
##
## u = 1
## v = 52.315
## f2 = 0.15
## sig.level = 0.05
## power = 0.8
##
##
## $n_from_pwr
## [1] 59
Kita simulasi data dengan \(p=5\) prediktor, menetapkan satu koefisien menghasilkan \(f^2\approx 0.15\), lalu cek power empiris.
set.seed(1)
sim_power <- function(n, B=500){
p <- 5
X <- matrix(rnorm(n*p), n, p)
# konstruksi beta agar partial f2 ~ 0.15 untuk x1 (aproksimasi): set beta1 dan noise
beta <- c(0.25, rep(0, p-1)) # nilai awal
y <- X %*% beta + rnorm(n, sd=1) # sigma=1
# calibrate roughly using first run's partial R2
fit <- lm(y ~ X)
# target via iteration simple (skip for speed): we accept approximate
rej <- logical(B)
for(b in 1:B){
Xb <- matrix(rnorm(n*p), n, p)
yb <- Xb %*% beta + rnorm(n, sd=1)
fitb <- lm(yb ~ Xb)
# uji koefisien X1 (kolom pertama)
sb <- summary(fitb)
tval <- coef(sb)[2,3]
pval <- 2*pt(abs(tval), df=n - p - 1, lower.tail = FALSE)
rej[b] <- (pval < 0.05)
}
mean(rej)
}
n_try <- ceiling(n_srs)
power_emp <- sim_power(n_try, B=200)
c(n_try = n_try, empirical_power = round(power_emp,3))
## n_try empirical_power
## 59.000 0.395
Catatan: simulasi di atas hanya ilustrasi kasar; untuk perancangan serius, gunakan parameter realistis (skala koefisien, korelasi antar-\(X\), dsb.) dan tingkatkan \(B\).
Ringkasan
pwr.f2.test.Dalam analisis data survei, heterogenitas responden sering muncul karena populasi terdiri dari subkelompok berbeda. Desain survei kompleks seperti stratified sampling secara eksplisit membagi populasi ke dalam strata yang relatif homogen untuk meningkatkan presisi dan representativitas.
Namun, kadang subkelompok tersebut tidak diketahui secara langsung atau ingin dieksplorasi lebih lanjut. Untuk itu, digunakan pendekatan Clustering Regression, yang menggabungkan teknik clustering dengan regresi untuk menangkap heterogenitas yang laten (tidak teramati jelas di awal).
Stratifikasi dalam Survei
Stratifikasi = membagi populasi menjadi subpopulasi (strata) yang relatif homogen (mis. jenis kelamin, provinsi, wilayah urban/rural). Analisis regresi dapat dilakukan per strata dengan bobot desain, kemudian digabung menggunakan pembobot populasi strata \(W_h\) (proporsi populasi strata \(h\)).
Estimator gabungan parameter regresi dapat ditulis secara skematis sebagai: \[ \widehat{\boldsymbol{\beta}}_{\text{strata}} = \sum_{h=1}^{H} W_h\;\widehat{\boldsymbol{\beta}}_h,\quad \text{dengan }\widehat{\boldsymbol{\beta}}_h \text{ adalah estimator regresi pada strata } h. \]
Dalam praktik, regresi sering dipasang per strata (atau menggunakan model multilevel dengan efek antar-strata) dan ketidakpastian total dihitung menggunakan sandwich variance yang sesuai dengan desain.
Clustering Regression
Clustering regression mengelompokkan data berdasarkan kemiripan karakteristik (mis. variabel penjelas atau gabungan \((\mathbf{x}, y)\)), lalu membangun model regresi untuk tiap cluster. Untuk unit \(i\) di cluster \(c\): \[ y_i = \mathbf{x}_i^{\top}\boldsymbol{\beta}_c + \varepsilon_i,\quad i\in \text{cluster }c. \]
Jika bobot survei \(w_i\) tersedia, estimator WLS per-cluster adalah: \[ \widehat{\boldsymbol{\beta}}_c = (\mathbf{X}_c^{\top}\mathbf{W}_c\mathbf{X}_c)^{-1}\mathbf{X}_c^{\top}\mathbf{W}_c\mathbf{y}_c, \] dengan \(\mathbf{W}_c = \mathrm{diag}(w_i)\) untuk \(i\in c\).
Hubungan Clustering Regression dengan Stratifikasi
| Aspek | Stratifikasi (Survei) | Clustering Regression |
|---|---|---|
| Pembentukan kelompok | Ditetapkan sebelum survei (mis. gender, provinsi) | Ditentukan dari data (data-driven) |
| Tujuan utama | Menjamin representasi & presisi | Menangkap heterogenitas yang tidak diketahui |
| Bobot sampling | Ada dari desain (\(w_i\)) | Perlu diintegrasikan ke proses (clustering & regresi) |
| Estimasi | Regresi per strata/multilevel | Regresi per cluster; bisa berbobot |
| Interpretasi | Strata punya makna kebijakan | Cluster bersifat eksploratif (bisa tak berlabel jelas) |
Intinya:
- Jika strata diketahui, analisis utama: regresi per
strata (atau multilevel) dengan bobot desain.
- Jika strata tidak diketahui, clustering
regression dapat menemukan “strata laten” yang
membantu pemodelan dan segmentasi.
Contoh Simulasi di R
Kita buat data simulasi dengan dua strata (A dan B),
lalu membandingkan: 1) Regresi per strata (benchmark desain),
2) Clustering regression tanpa informasi
strata,
3) Clustering regression berbobot (menggunakan
bobot sampling pada proses clustering dan regresi).
set.seed(123)
# Ukuran sampel per strata
nA <- 100; nB <- 100
# Strata A: kemiringan (slope) 0.5, intersep 2
xA <- rnorm(nA, mean = 5, sd = 1)
yA <- 2 + 0.5 * xA + rnorm(nA, mean = 0, sd = 1)
# Strata B: kemiringan 0.2, intersep 5
xB <- rnorm(nB, mean = 10, sd = 1)
yB <- 5 + 0.2 * xB + rnorm(nB, mean = 0, sd = 1)
# Misalkan bobot sampling berbeda antar strata (A overweighted)
data <- data.frame(
y = c(yA, yB),
x = c(xA, xB),
strata = rep(c("A","B"), c(nA, nB)),
weight = c(rep(1.5, nA), rep(0.5, nB))
)
head(data)
## y x strata weight
## 1 3.509356 4.439524 A 1.5
## 2 4.641795 4.769823 A 1.5
## 3 5.032662 6.558708 A 1.5
## 4 4.187712 5.070508 A 1.5
## 5 3.613025 5.129288 A 1.5
## 6 5.312505 6.715065 A 1.5
modA <- lm(y ~ x, data = subset(data, strata == "A"), weights = weight)
modB <- lm(y ~ x, data = subset(data, strata == "B"), weights = weight)
summary(modA)
##
## Call:
## lm(formula = y ~ x, data = subset(data, strata == "A"), weights = weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3360 -0.8371 -0.1072 0.7111 4.0299
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1596 0.5526 3.908 0.000172 ***
## x 0.4475 0.1069 4.187 6.17e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.189 on 98 degrees of freedom
## Multiple R-squared: 0.1518, Adjusted R-squared: 0.1431
## F-statistic: 17.53 on 1 and 98 DF, p-value: 6.173e-05
summary(modB)
##
## Call:
## lm(formula = y ~ x, data = subset(data, strata == "B"), weights = weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.72616 -0.46222 0.00566 0.52410 1.78655
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4603 1.1217 4.868 4.32e-06 ***
## x 0.1509 0.1104 1.368 0.175
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7375 on 98 degrees of freedom
## Multiple R-squared: 0.01873, Adjusted R-squared: 0.008717
## F-statistic: 1.871 on 1 and 98 DF, p-value: 0.1745
Sebagai pembanding gabungan sederhana (bukan pengganti desain lengkap), kita bisa melihat model yang mengizinkan intersep dan slope berbeda antar-strata melalui interaksi:
mod_interact <- lm(y ~ strata * x, data = data, weights = weight)
summary(mod_interact)
##
## Call:
## lm(formula = y ~ strata * x, data = data, weights = weight)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -2.3360 -0.6579 -0.0282 0.5852 4.0299
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.15956 0.45987 4.696 4.98e-06 ***
## strataB 3.30078 1.57336 2.098 0.0372 *
## x 0.44753 0.08894 5.032 1.09e-06 ***
## strataB:x -0.29659 0.17269 -1.717 0.0875 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9893 on 196 degrees of freedom
## Multiple R-squared: 0.5844, Adjusted R-squared: 0.578
## F-statistic: 91.85 on 3 and 196 DF, p-value: < 2.2e-16
Kita lakukan K-means pada \((x, y)\) lalu pasang regresi per cluster tanpa bobot sebagai baseline eksploratif.
set.seed(123)
km <- kmeans(data[, c("x","y")], centers = 2, nstart = 20)
data$cluster <- factor(km$cluster)
modC1 <- lm(y ~ x, data = subset(data, cluster == 1))
modC2 <- lm(y ~ x, data = subset(data, cluster == 2))
summary(modC1)
##
## Call:
## lm(formula = y ~ x, data = subset(data, cluster == 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9073 -0.6835 -0.0875 0.5806 3.2904
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1596 0.5526 3.908 0.000172 ***
## x 0.4475 0.1069 4.187 6.17e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9707 on 98 degrees of freedom
## Multiple R-squared: 0.1518, Adjusted R-squared: 0.1431
## F-statistic: 17.53 on 1 and 98 DF, p-value: 6.173e-05
summary(modC2)
##
## Call:
## lm(formula = y ~ x, data = subset(data, cluster == 2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4412 -0.6537 0.0080 0.7412 2.5266
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4603 1.1217 4.868 4.32e-06 ***
## x 0.1509 0.1104 1.368 0.175
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.043 on 98 degrees of freedom
## Multiple R-squared: 0.01873, Adjusted R-squared: 0.008717
## F-statistic: 1.871 on 1 and 98 DF, p-value: 0.1745
table(TrueStrata = data$strata, FoundCluster = data$cluster)
## FoundCluster
## TrueStrata 1 2
## A 100 0
## B 0 100
Kita gunakan weighted K-means dari paket
flexclust (fungsi cclust) untuk memasukkan
bobot sampling saat clustering. Kemudian,
regresi per cluster juga menggunakan bobot.
# install.packages("flexclust") # aktifkan baris ini jika paket belum terpasang
library(flexclust)
cl_w <- cclust(data[, c("x","y")], k = 2, weights = data$weight, save.data = TRUE)
data$cluster_w <- factor(cl_w@cluster)
modW1 <- lm(y ~ x, data = subset(data, cluster_w == 1), weights = weight)
modW2 <- lm(y ~ x, data = subset(data, cluster_w == 2), weights = weight)
summary(modW1)
##
## Call:
## lm(formula = y ~ x, data = subset(data, cluster_w == 1), weights = weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.72616 -0.46222 0.00566 0.52410 1.78655
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4603 1.1217 4.868 4.32e-06 ***
## x 0.1509 0.1104 1.368 0.175
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7375 on 98 degrees of freedom
## Multiple R-squared: 0.01873, Adjusted R-squared: 0.008717
## F-statistic: 1.871 on 1 and 98 DF, p-value: 0.1745
summary(modW2)
##
## Call:
## lm(formula = y ~ x, data = subset(data, cluster_w == 2), weights = weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3360 -0.8371 -0.1072 0.7111 4.0299
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1596 0.5526 3.908 0.000172 ***
## x 0.4475 0.1069 4.187 6.17e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.189 on 98 degrees of freedom
## Multiple R-squared: 0.1518, Adjusted R-squared: 0.1431
## F-statistic: 17.53 on 1 and 98 DF, p-value: 6.173e-05
table(TrueStrata = data$strata, WeightedCluster = data$cluster_w)
## WeightedCluster
## TrueStrata 1 2
## A 0 100
## B 100 0
Visualisasi
# install.packages("ggplot2") # jika perlu
library(ggplot2)
p1 <- ggplot(data, aes(x = x, y = y, color = strata, size = weight)) +
geom_point(alpha = 0.6) +
labs(title = "Data Simulasi per Strata (dengan Bobot)")
p2 <- ggplot(data, aes(x = x, y = y, color = cluster)) +
geom_point(alpha = 0.6) +
labs(title = "Clustering Tanpa Bobot (K-means)")
p3 <- ggplot(data, aes(x = x, y = y, color = cluster_w, size = weight)) +
geom_point(alpha = 0.6) +
labs(title = "Clustering Berbobot (cclust)")
p1
p2
p3
Diskusi
Catatan: Clustering regression tidak menggantikan penaksiran berbasis desain ketika tujuan utamanya adalah inferensi populasi. Namun, ia berguna untuk segmentasi dan pemodelan ketika struktur heterogenitas belum jelas.
Kesimpulan
Lampiran: Mixture Regression Singkat
Sebagai gambaran, pendekatan mixture regression
memodelkan bahwa setiap observasi berasal dari salah satu dari \(K\) komponen dengan probabilitas \(\pi_k\), dan regresi khusus komponen \(\boldsymbol{\beta}_k\): \[
f(y_i\mid \mathbf{x}_i) = \sum_{k=1}^{K} \pi_k\;\mathcal{N}\big(y_i\mid
\mathbf{x}_i^{\top}\boldsymbol{\beta}_k,\,\sigma_k^2\big).
\] Estimasi biasanya melalui EM. Dalam konteks survei, bobot
\(w_i\) dapat dimasukkan pada langkah
E/M (mis. weighted EM). Paket seperti flexmix
menyediakan mixture of regressions (tanpa bobot secara
default), sedangkan weighted EM memerlukan penyesuaian
manual.
Design Effect (DEFF) & Effective Sample Size
\[ \mathrm{DEFF}=\frac{\mathrm{Var}_{\text{design}}(\hat{\theta})}{\mathrm{Var}_{\text{SRS}}(\hat{\theta})},\qquad n_{\text{eff}}=\frac{n}{\mathrm{DEFF}}. \]
Koreksi Populasi Berhingga (FPC)
Untuk SRS tanpa pengembalian: \[ \sqrt{\mathrm{Var}(\bar{y})}=\sqrt{\left(1-\frac{n}{N}\right)\frac{S_y^2}{n}}. \]
Tips
svydesign() & svyglm() untuk
SE yang konsisten berbasis-desain.svyglm(family=...) tetap relevan.Ekstensi Singkat
as.svrepdesign().Referensi Singkat