Latihan dan Tugas
Melly Mustikasari 140720250011
1 Konsep dan Teori Analisis Regresi Tingkat Lanjut
1.1 Konsep Analisis Regresi dan Korelasi
Korelasi
Korelasi adalah ukuran kekuatan dan arah hubungan
linier antara dua variabel.
Nilai korelasi biasanya dinyatakan dengan koefisien r
(Pearson), yang berkisar dari:
- r = +1 : hubungan positif sempurna (jika X naik, Y selalu naik).
- r = 0 : tidak ada hubungan linier.
- r = -1 : hubungan negatif sempurna (jika X naik, Y selalu turun).
Korelasi hanya mengukur hubungan, bukan
sebab-akibat.
Contoh: tinggi badan dan berat badan biasanya berkorelasi positif,
tetapi tidak berarti tinggi menyebabkan berat.
Regresi
Regresi adalah metode statistik untuk memodelkan hubungan antara variabel:
- Variabel dependen (Y): yang ingin dijelaskan/diprediksi.
- Variabel independen (X): yang digunakan untuk menjelaskan/prediksi Y.
Model regresi linier sederhana: \[ Y = \beta_0 + \beta_1 X + \varepsilon \]
- \(\beta_0\) : intersep (nilai Y
saat X = 0)
- \(\beta_1\) : koefisien regresi,
menyatakan perubahan rata-rata Y ketika X naik satu unit
- \(\varepsilon\) : error (faktor lain yang tidak dimodelkan)
Regresi bisa dipakai untuk prediksi, sedangkan korelasi hanya menunjukkan hubungan.
Perbedaan Utama
\[ \begin{array}{|l|l|l|} \hline \textbf{Aspek} & \textbf{Korelasi} & \textbf{Regresi} \\ \hline \text{Tujuan} & \text{Mengukur keeratan hubungan} & \text{Membangun model prediksi/penjelasan} \\ \hline \text{Arah hubungan} & \text{Simetris (X dan Y setara)} & \text{Asimetris (Y sebagai respon, X sebagai prediktor)} \\ \hline \text{Output utama} & \text{Koefisien korelasi } r & \text{Persamaan regresi } \hat{Y} = \beta_0 + \beta_1 X \\ \hline \end{array} \]
Contoh Aplikasi d R
Data
Gunakan dataset mtcars untuk melihat hubungan mpg (miles per gallon) dengan wt (berat mobil)
## mpg wt
## Mazda RX4 21.0 2.620
## Mazda RX4 Wag 21.0 2.875
## Datsun 710 22.8 2.320
## Hornet 4 Drive 21.4 3.215
## Hornet Sportabout 18.7 3.440
## Valiant 18.1 3.460Korelasi
##
## Pearson's product-moment correlation
##
## data: mtcars$wt and mtcars$mpg
## t = -9.559, df = 30, p-value = 1.294e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9338264 -0.7440872
## sample estimates:
## cor
## -0.8676594Regresi
##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10Visualisasi
library(ggplot2)
ggplot(mtcars, aes(x = wt, y = mpg)) +
geom_point() +
geom_smooth(method = "lm", se = TRUE, color = "blue") +
labs(title = "Hubungan antara Berat Mobil (wt) dan Efisiensi Bahan Bakar (mpg)",
x = "Berat Mobil (1000 lbs)",
y = "Miles per Gallon (mpg)")
1.2 Model Regresi Linear
1.2.1 Estimator OLS & Geometri Proyeksi
Estimator ordinary least squares (OLS) meminimalkan jumlah kuadrat residu: \[ \hat{\beta}_{OLS} = \arg\min_{\beta} \, \| y - X\beta \|^2 = (X^{\top}X)^{-1}X^{\top}y \] \[ H = X(X^{\top}X)^{-1}X^{\top}, \quad \hat{y} = Hy, \quad e = (I - H)y \]
Penjelasan
- \(H =\) Hat matrix atau matriks proyeksi, memetakan \(y\) ke ruang kolom \(X\).
- \(\hat{y} = Hy\) = nilai prediksi dari model regresi.
- \(e = (I - H)y\) = residual (selisih antara observasi \(y\) dan prediksi \(\hat{y}\)).
1.2.2 Sifat Asimtotik & Varian Koefisien
Di bawah asumsi Gauss–Markov (linearitas, eksogenitas, homoskedastisitas, dan tidak ada multikolinearitas sempurna), ̂ 𝛽OLS adalah BLUE. Varians koefisien dan estimator ragam:
\[ \text{Var}(\hat{\beta}_{OLS}) = \sigma^2 (X^T X)^{-1}, \quad \hat{\sigma}^2 = \frac{\| e \|^2}{n - p}. \] Sifat Asimtotik : Untuk \(n\) yang besar, berlaku sifat asimtotik berikut:
\[ \sqrt{n} (\hat{\beta} - \beta) \overset{d}{\longrightarrow} \mathcal{N}(0, V), \]
sehingga inferensi berbasis normal/t berlaku ketika \(n\) besar.
1.3 Generalisasi : WLS & GLS
Jika galat tidak homoskedastik/berkorelasi, maka: \[ \text{Var}(\varepsilon) = \sigma^2 \Omega \neq \sigma^2 I. \]
Weighted least squares (WLS) untuk varians berbeda-beda (\(\Omega\) diagonal):
\[ \hat{\beta}_{WLS} = (X^T W X)^{-1} X^T W y, \quad W = \Omega^{-1}. \]
Generalized least squares (GLS) (\(\Omega\) umum):
\[ \hat{\beta}_{GLS} = (X^T \Omega^{-1} X)^{-1} X^T \Omega^{-1} y. \]
1.4 Robust Standard Errors (White Sandwich)
Tanpa menspesifikasi bentuk \(\Omega\), kovarians heteroskedastisity-consistent (HC) ala White adalah:
\[ \text{Var}(\hat{\beta}) = (X^T X)^{-1} \left( \sum_{i=1}^{n} x_i x_i^T \hat{\varepsilon}_i^2 \right) (X^T X)^{-1}. \]
Ini memberi robust SE yang valid terhadap heteroskedastisitas (HC0–HC3).
1.5 Multikolinearitas & Diagnostik
Ukuran pengaruh Cook’s distance: \[ D_i = \frac{e_i^2}{p \hat{\sigma}^2} \cdot \frac{h_{ii}}{(1 - h_{ii})^2}. \]
Gunakan plot residu, plot leverage vs residu, QQ-plot, serta inspeksi \(D_i\) untuk memeriksa asumsi dan outlier.
1.6 Bias–Variance Trade-off & Regularisasi (Ikhtisar)
Model makin kompleks (banyak parameter/basis) cenderung bias rendah namun varians tinggi. Ridge dan Lasso menambahkan penalti:
\[ \text{Ridge: } \hat{\beta} = \arg \min_{\beta} \| y - X\beta \|^2 + \lambda \|\beta\|_2^2, \qquad \text{Lasso: } \hat{\beta} = \arg \min_{\beta} \| y - X\beta \|^2 + \lambda \|\beta\|_1. \]
1.7 Nonlinearitas: Ekspansi Basis & Spline
Kita dapat menambah term polinomial atau basis spline sehingga tetap linear dalam parameter: \[ y = \beta_0 + \sum_k \beta_k b_k(x) + \varepsilon. \]
1.8 Studi Kasus 1: Regresi pada Data Nyata (mtcars)
1.8.1 Persiapan Data dan Eksplorasi
dat1 <- mtcars %>%
rownames_to_column("car") %>%
as_tibble() %>%
mutate(across(c(mpg, disp, hp, wt, qsec), as.double))
head(dat1)## # A tibble: 6 × 12
## car mpg cyl disp hp drat wt qsec vs am gear carb
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Mazda RX4 21 6 160 110 3.9 2.62 16.5 0 1 4 4
## 2 Mazda RX4 W… 21 6 160 110 3.9 2.88 17.0 0 1 4 4
## 3 Datsun 710 22.8 4 108 93 3.85 2.32 18.6 1 1 4 1
## 4 Hornet 4 Dr… 21.4 6 258 110 3.08 3.22 19.4 1 0 3 1
## 5 Hornet Spor… 18.7 8 360 175 3.15 3.44 17.0 0 0 3 2
## 6 Valiant 18.1 6 225 105 2.76 3.46 20.2 1 0 3 1Plot hubungan awal
GG1 <- ggplot(dat1, aes(wt, mpg)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = FALSE) +
labs(x = "Weight (1000 lbs)", y = "Miles per Gallon (mpg)",
title = "Hubungan mpg vs wt pada mtcars")
GG1
1.8.2 Estimasi OLS dan Ringkasan
##
## Call:
## lm(formula = mpg ~ wt + hp + disp, data = dat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.891 -1.640 -0.172 1.061 5.861
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.105505 2.110815 17.579 < 2e-16 ***
## wt -3.800891 1.066191 -3.565 0.00133 **
## hp -0.031157 0.011436 -2.724 0.01097 *
## disp -0.000937 0.010350 -0.091 0.92851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.639 on 28 degrees of freedom
## Multiple R-squared: 0.8268, Adjusted R-squared: 0.8083
## F-statistic: 44.57 on 3 and 28 DF, p-value: 8.65e-11## # A tibble: 4 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 37.1 2.11 17.6 1.16e-16 32.8 41.4
## 2 wt -3.80 1.07 -3.56 1.33e- 3 -5.98 -1.62
## 3 hp -0.0312 0.0114 -2.72 1.10e- 2 -0.0546 -0.00773
## 4 disp -0.000937 0.0103 -0.0905 9.29e- 1 -0.0221 0.02031.8.3 Diagnostik Dasar (Residu, QQ-plot, Cook)
diag_df <- tibble(
fitted = fitted(fit1), resid = resid(fit1),
stdres = rstandard(fit1), hat = hatvalues(fit1), cook = cooks.distance(fit1)
)
p1 <- ggplot(diag_df, aes(fitted, resid)) +
geom_point(alpha = 0.7) +
geom_hline(yintercept = 0, linetype = 2) +
labs(x = "Fitted", y = "Residuals")
p2 <- ggplot(diag_df, aes(sample = stdres)) +
stat_qq() + stat_qq_line() + labs(x = "Theoretical", y = "Standardized Residuals")
p1; p2
qplot(seq_along(diag_df$cook), diag_df$cook, geom = c("line","point"),
xlab = "Index", ylab = "Cook's D")
1.8.5 Robust SE (White HC3) & Interpretasi
##
## 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##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.10550527 2.71026981 13.6907 6.257e-14 ***
## wt -3.80089058 1.08086207 -3.5165 0.00151 **
## hp -0.03115655 0.01437457 -2.1675 0.03886 *
## disp -0.00093701 0.01012138 -0.0926 0.92690
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11.8.6 Nonlinearitas: Polinomial vs Spline pada wt
library(splines)
fit_poly <- lm(mpg ~ poly(wt, 2, raw = TRUE) + hp + disp, data = dat1)
fit_spl <- lm(mpg ~ bs(wt, df = 5) + hp + disp, data = dat1)
AIC(fit1, fit_poly, fit_spl) %>% arrange(AIC)## df AIC
## fit_poly 6 150.7305
## fit_spl 9 154.3323
## fit1 5 158.6430## df BIC
## fit_poly 6 159.5250
## fit1 5 165.9717
## fit_spl 9 167.5239# Kurva prediksi
library(modelr)
library(purrr)
grid_wt <- tibble(wt = seq(min(dat1$wt), max(dat1$wt), length.out = 200),
hp = mean(dat1$hp), disp = mean(dat1$disp))
pred <- grid_wt %>%
mutate(OLS = predict(fit1, newdata = grid_wt),
Poly = predict(fit_poly, newdata = grid_wt),
Spline = predict(fit_spl, newdata = grid_wt)) %>%
pivot_longer(OLS:Spline, names_to = "model", values_to = "yhat")
GG2 <- ggplot() +
geom_point(data = dat1, aes(wt, mpg), alpha = 0.3) +
geom_line(data = pred, aes(wt, yhat, linetype = model), linewidth = 0.9) +
labs(x = "wt", y = "mpg", title = "Perbandingan OLS vs Polinomial vs Spline")
GG2
1.8.7 Validasi Silang (K-fold) untuk Perbandingan Model
k <- 5
cv_res <- crossv_kfold(dat1, k = k)
rmse_fun <- function(train, test, formula) {
m <- lm(formula, data = as.data.frame(train))
sqrt(mean((as.data.frame(test)$mpg - predict(m, newdata = as.data.frame(test)))^2))
}
forms <- list(
OLS = mpg ~ wt + hp + disp,
Poly = mpg ~ poly(wt, 2, raw = TRUE) + hp + disp,
Spline = mpg ~ bs(wt, df = 5) + hp + disp
)
cv_tbl <- map_dfr(names(forms), function(nm){
fml <- forms[[nm]]
tibble(model = nm, fold = seq_len(k),
rmse = map2_dbl(cv_res$train, cv_res$test, ~ rmse_fun(.x, .y, fml)))
})
cv_tbl %>% group_by(model) %>% summarise(mean_rmse = mean(rmse), sd_rmse = sd(rmse)) %>% arrange(mean_rmse)## # A tibble: 3 × 3
## model mean_rmse sd_rmse
## <chr> <dbl> <dbl>
## 1 Poly 2.64 0.816
## 2 OLS 2.67 1.06
## 3 Spline 2.85 0.8581.9 Studi Kasus 2: Simulasi Heteroskedastisitas (WLS vs OLSvs Robust)
n <- 400
x1 <- runif(n, 0, 10)
x2 <- rnorm(n, 5, 2)
sigma_i <- 0.6 + 0.2 * x1 # var naik dengan x1
err <- rnorm(n, 0, sigma_i)
y <- 2 + 1.5*x1 - 1.2*x2 + err
sim <- tibble(y, x1, x2, w = 1/sigma_i^2)
fit_ols <- lm(y ~ x1 + x2, data = sim)
fit_wls <- lm(y ~ x1 + x2, data = sim, weights = w)
# Tabel perbandingan SE
comp_se <- tibble(
term = names(coef(fit_ols)),
OLS_SE = coef(summary(fit_ols))[, "Std. Error"],
Robust_SE = sqrt(diag(vcovHC(fit_ols, type = "HC3"))),
WLS_SE = coef(summary(fit_wls))[, "Std. Error"]
)
comp_se## # A tibble: 3 × 4
## term OLS_SE Robust_SE WLS_SE
## <chr> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.241 0.193 0.152
## 2 x1 0.0277 0.0280 0.0243
## 3 x2 0.0372 0.0333 0.0266##
## studentized Breusch-Pagan test
##
## data: fit_ols
## BP = 31.359, df = 2, p-value = 1.551e-07ggplot(data.frame(fitted=fitted(fit_ols), resid=resid(fit_ols)), aes(fitted, resid)) +
geom_point(alpha=.6) + geom_smooth(se=FALSE, method="loess") +
geom_hline(yintercept = 0, linetype=2) + labs(x="Fitted", y="Residuals")
1.10 Catatan
- Mulai dari OLS + diagnostik; bila heteroskedastik, pilih robust SE atau WLS/GLS sesuai struktur ragam.
- Periksa multikolinearitas (VIF), leverage, Cook’s D, dan normalitas residu.
- Gunakan ekspansi basis (polinomial/spline) saat ada indikasi nonlinearitas.
- Validasi model (CV) dan bandingkan AIC/BIC untuk kompleksitas vs akurasi.
- Untuk dimensi tinggi, gunakan regularisasi (Ridge/Lasso) dan pipeline (standarisasi, CV).
1.11 GLS dengan Korelasi AR(1): Contoh dan Simulasi
Model galat mengikuti proses AR(1):
\[
\varepsilon_t = \rho \varepsilon_{t-1} + u_t, \quad \text{dengan } u_t
\sim \mathcal{N}(0, \sigma_u^2).
\]
Kovarians residual memiliki bentuk Toeplitz:
\[
\Omega_{ts} = \sigma_u^2 \rho^{|t-s|}/(1-\rho^2).
\]
GLS menimbang data dengan \(\Omega^{-1}\) (atau transformasi ekuivalen) agar estimasi efisien.
1.11.1 Simulasi data dengan galat AR(1)
set.seed(42)
T <- 400
x <- runif(T, -2, 2)
rho <- 0.6
u <- rnorm(T, 0, 1)
err <- numeric(T)
for (t in 2:T) err[t] <- rho*err[t-1] + u[t]
y <- 1 + 2*x + err
dar1 <- tibble::tibble(t = 1:T, x, y)
head(dar1)## # A tibble: 6 × 3
## t x y
## <int> <dbl> <dbl>
## 1 1 1.66 4.32
## 2 2 1.75 4.83
## 3 3 -0.855 0.661
## 4 4 1.32 6.53
## 5 5 0.567 2.49
## 6 6 0.0764 0.2131.11.2 OLS vs GLS (nlme::gls dengan corAR1)
library(nlme)
# OLS biasa
fit_ols_ar1 <- lm(y ~ x, data = dar1)
# GLS dengan korelasi AR(1) di residual
fit_gls_ar1 <- gls(y ~ x, data = dar1,
correlation = corAR1(form = ~ t))
summary(fit_ols_ar1)##
## Call:
## lm(formula = y ~ x, data = dar1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3360 -0.8231 0.0026 0.8211 3.9013
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.95441 0.06043 15.79 <2e-16 ***
## x 2.06045 0.05162 39.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.209 on 398 degrees of freedom
## Multiple R-squared: 0.8002, Adjusted R-squared: 0.7997
## F-statistic: 1594 on 1 and 398 DF, p-value: < 2.2e-16## Generalized least squares fit by REML
## Model: y ~ x
## Data: dar1
## AIC BIC logLik
## 1139.541 1155.487 -565.7706
##
## Correlation Structure: AR(1)
## Formula: ~t
## Parameter estimate(s):
## Phi
## 0.5772332
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.9552097 0.11651869 8.19791 0
## x 2.0449090 0.03743358 54.62766 0
##
## Correlation:
## (Intr)
## x 0.004
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.752477543 -0.674750151 -0.007333444 0.679546648 3.225627947
##
## Residual standard error: 1.210605
## Degrees of freedom: 400 total; 398 residual1.11.3 Diagnostik serial correlation
par(mfrow = c(1,2))
acf(resid(fit_ols_ar1), main = "ACF Residuals (OLS)")
acf(resid(fit_gls_ar1, type = "normalized"), main = "ACF Residuals (GLS)")
##
## Durbin-Watson test
##
## data: fit_ols_ar1
## DW = 0.84993, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 01.11.4 Alternatif: Quasi-differencing (Cochrane–Orcutt) manual
# Estimasi rho dari regresi resid_t ~ resid_{t-1}
res <- resid(fit_ols_ar1)
rho_hat <- coef(lm(res[-1] ~ 0 + res[-length(res)]))[1]
# Quasi-differencing (transformasi data)
y_star <- dar1$y[-1] - rho_hat*dar1$y[-T]
x_star <- dar1$x[-1] - rho_hat*dar1$x[-T]
fit_qd <- lm(y_star ~ x_star)
summary(fit_qd)##
## Call:
## lm(formula = y_star ~ x_star)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6691 -0.6723 -0.0232 0.6970 3.2363
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.40603 0.04955 8.194 3.52e-15 ***
## x_star 2.04499 0.03757 54.439 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9898 on 397 degrees of freedom
## Multiple R-squared: 0.8819, Adjusted R-squared: 0.8816
## F-statistic: 2964 on 1 and 397 DF, p-value: < 2.2e-16Catatan: gls() memperkirakan 𝜌 dan ragam secara simultan dalam kerangka likelihood sehingga sering lebih stabil dibanding pendekatan dua tahap.
1.12 Latihan
# Paket yang diperlukan
pkgs <- c("dplyr","ggplot2","car","lmtest","sandwich","MASS","boot","broom")
new <- pkgs[!(pkgs %in% installed.packages()[,"Package"])]
if(length(new)) install.packages(new, dependencies = TRUE)
lapply(pkgs, library, character.only = TRUE)## [[1]]
## [1] "nlme" "purrr" "modelr" "tidyr" "splines" "sandwich"
## [7] "lmtest" "zoo" "broom" "tibble" "dplyr" "ggplot2"
## [13] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [19] "base"
##
## [[2]]
## [1] "nlme" "purrr" "modelr" "tidyr" "splines" "sandwich"
## [7] "lmtest" "zoo" "broom" "tibble" "dplyr" "ggplot2"
## [13] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [19] "base"
##
## [[3]]
## [1] "car" "carData" "nlme" "purrr" "modelr" "tidyr"
## [7] "splines" "sandwich" "lmtest" "zoo" "broom" "tibble"
## [13] "dplyr" "ggplot2" "stats" "graphics" "grDevices" "utils"
## [19] "datasets" "methods" "base"
##
## [[4]]
## [1] "car" "carData" "nlme" "purrr" "modelr" "tidyr"
## [7] "splines" "sandwich" "lmtest" "zoo" "broom" "tibble"
## [13] "dplyr" "ggplot2" "stats" "graphics" "grDevices" "utils"
## [19] "datasets" "methods" "base"
##
## [[5]]
## [1] "car" "carData" "nlme" "purrr" "modelr" "tidyr"
## [7] "splines" "sandwich" "lmtest" "zoo" "broom" "tibble"
## [13] "dplyr" "ggplot2" "stats" "graphics" "grDevices" "utils"
## [19] "datasets" "methods" "base"
##
## [[6]]
## [1] "MASS" "car" "carData" "nlme" "purrr" "modelr"
## [7] "tidyr" "splines" "sandwich" "lmtest" "zoo" "broom"
## [13] "tibble" "dplyr" "ggplot2" "stats" "graphics" "grDevices"
## [19] "utils" "datasets" "methods" "base"
##
## [[7]]
## [1] "boot" "MASS" "car" "carData" "nlme" "purrr"
## [7] "modelr" "tidyr" "splines" "sandwich" "lmtest" "zoo"
## [13] "broom" "tibble" "dplyr" "ggplot2" "stats" "graphics"
## [19] "grDevices" "utils" "datasets" "methods" "base"
##
## [[8]]
## [1] "boot" "MASS" "car" "carData" "nlme" "purrr"
## [7] "modelr" "tidyr" "splines" "sandwich" "lmtest" "zoo"
## [13] "broom" "tibble" "dplyr" "ggplot2" "stats" "graphics"
## [19] "grDevices" "utils" "datasets" "methods" "base"1.Simulasi Data
Tujuan: membuat data dengan 3 prediktor (X1, X2, X3), disuntik multikolinearitas ringan (X1 & X2 berkorelasi) dan heteroskedastisitas (ragam error meningkat saat |X1| besar), serta sedikit outlier.
n <- 200
X1 <- rnorm(n, 0, 1)
# X2 berkorelasi ~0.7 dengan X1
rho <- 0.7
X2 <- rho*X1 + sqrt(1-rho^2)*rnorm(n)
X3 <- rnorm(n)
# Heteroskedastisitas: sd error tergantung |X1|
eps <- rnorm(n, sd = 0.5 + 0.5*abs(X1))
# Model benar (ground truth): y = 2 + 1.5*X1 - 0.8*X2 + 0.5*X3 + eps
y <- 2 + 1.5*X1 - 0.8*X2 + 0.5*X3 + eps
# Tambah beberapa outlier pada respon
idx_out <- sample(1:n, 3)
y[idx_out] <- y[idx_out] + rnorm(3, mean = 6, sd = 1)
dat <- data.frame(y, X1, X2, X3)
summary(dat)## y X1 X2 X3
## Min. :-2.3696 Min. :-3.01793 Min. :-2.65816 Min. :-2.928477
## 1st Qu.: 0.8631 1st Qu.:-0.70921 1st Qu.:-0.75951 1st Qu.:-0.743733
## Median : 1.9085 Median :-0.06983 Median : 0.06149 Median : 0.009401
## Mean : 1.8694 Mean :-0.12802 Mean :-0.03760 Mean :-0.044587
## 3rd Qu.: 2.6646 3rd Qu.: 0.56642 3rd Qu.: 0.71289 3rd Qu.: 0.582471
## Max. : 7.7179 Max. : 2.22353 Max. : 2.57526 Max. : 2.6337102. Estimasi Parameter (OLS)
##
## Call:
## lm(formula = y ~ X1 + X2 + X3, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3740 -0.5558 -0.0325 0.4862 5.7915
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.04203 0.07744 26.369 < 2e-16 ***
## X1 1.40441 0.10917 12.864 < 2e-16 ***
## X2 -0.79319 0.09733 -8.150 4.23e-14 ***
## X3 0.50870 0.07732 6.579 4.22e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.081 on 196 degrees of freedom
## Multiple R-squared: 0.5249, Adjusted R-squared: 0.5177
## F-statistic: 72.19 on 3 and 196 DF, p-value: < 2.2e-16## # 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.525 0.518 1.08 72.2 1.71e-31 3 -297. 605. 621.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>Catatan interpretasi cepat: Koefisien menunjukkan perubahan rata-rata y untuk kenaikan 1 unit pada prediktor (dengan prediktor lain konstan). Lihat Pr(>|t|) untuk signifikansi.
3.Uji Hipotesis
3.1 Uji serentak (overall F-test) Sudah tersedia pada output summary(mod0) (Signifikansi model secara keseluruhan).
3.2 Uji parsial (t-test per koefisien) Juga ada di summary(mod0).
3.3 Uji hipotesis gabungan (contoh:
##
## 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 370.42
## 2 196 229.13 2 141.29 60.432 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.4 Uji dengan SE robust (menghadapi heteroskedastisitas)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.042033 0.074088 27.5623 < 2.2e-16 ***
## X1 1.404415 0.130454 10.7656 < 2.2e-16 ***
## X2 -0.793190 0.089124 -8.8999 3.756e-16 ***
## X3 0.508697 0.068733 7.4011 3.883e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14. Uji Asumsi 4.1 Multikolinearitas
## X1 X2 X3
## 1.751257 1.763584 1.009529## [1] 2.399881Rule of thumb: VIF > 10 (atau > 5) → indikasi kolinearitas meningkat.
4.2 Heteroskedastisitas
##
## studentized Breusch-Pagan test
##
## data: mod0
## BP = 4.3812, df = 3, p-value = 0.2231# White test (via bptest dengan kuadrat fitted)
lmtest::bptest(mod0, ~ fitted(mod0) + I(fitted(mod0)^2))##
## studentized Breusch-Pagan test
##
## data: mod0
## BP = 3.299, df = 2, p-value = 0.1921Jika signifikan, gunakan SE robust (lihat 3.4) untuk inferensi yang lebih andal.
4.3 Normalitas residual
##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.89843, p-value = 1.999e-10
4.4 Autokorelasi residual (opsional) (Biasanya untuk data runtun waktu)
##
## Durbin-Watson test
##
## data: mod0
## DW = 2.0477, p-value = 0.629
## alternative hypothesis: true autocorrelation is greater than 04.5 Uji spesifikasi (linearitas) – Ramsey RESET
##
## RESET test
##
## data: mod0
## RESET = 1.9032, df1 = 2, df2 = 194, p-value = 0.15194.6 Plot Diagnostik Standar
5. Seleksi Model Strategi: mulai dari model kandidat (termasuk transformasi sederhana) lalu gunakan AIC (stepwise), bandingkan dengan model awal, dan validasi via CV.
dat2 <- dat %>%
mutate(X1_sq = X1^2, X2_sq = X2^2, X3_sq = X3^2,
X1X2 = X1*X2, X1X3 = X1*X3, X2X3 = X2*X3)
mod_full <- lm(y ~ X1 + X2 + X3 + X1_sq + X2_sq + X3_sq + X1X2 + X1X3 + X2X3, data = dat2)
AIC(mod0); AIC(mod_full)## [1] 604.768## [1] 609.6585# 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 + X1X3, data = dat2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3372 -0.5195 -0.0579 0.4668 5.7224
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.04365 0.07704 26.528 < 2e-16 ***
## X1 1.39614 0.10870 12.844 < 2e-16 ***
## X2 -0.78368 0.09696 -8.082 6.55e-14 ***
## X3 0.50795 0.07692 6.604 3.71e-10 ***
## X1X3 0.14131 0.08050 1.755 0.0808 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.076 on 195 degrees of freedom
## Multiple R-squared: 0.5323, Adjusted R-squared: 0.5227
## F-statistic: 55.49 on 4 and 195 DF, p-value: < 2.2e-16## # A tibble: 1 × 4
## r.squared adj.r.squared AIC BIC
## <dbl> <dbl> <dbl> <dbl>
## 1 0.532 0.523 604. 623.Validasi silang (10-fold CV) untuk membandingkan RMSE:
# Fungsi CV RMSE menggunakan boot::cv.glm
cv_rmse <- function(model, data, K=10){
# cv.glm mengembalikan delta: [1] raw CV, [2] adjusted CV
set.seed(123)
res <- boot::cv.glm(data = data, glmfit = model, K = K)
sqrt(res$delta[1])
}
rmse_mod0 <- cv_rmse(mod0, dat)
rmse_modstep<- cv_rmse(mod_step, dat2)
c(RMSE_mod0 = rmse_mod0, RMSE_modStep = rmse_modstep)## RMSE_mod0 RMSE_modStep
## NaN NaNPilih model dengan trade-off terbaik antara kesederhanaan dan kinerja (RMSE, AIC/BIC, Adj, R-squared.
6. Evaluasi Outlier & Titik Berpengaruh Gunakan studentized residuals, leverage (hat), dan Cook’s distance. Tambahkan Bonferroni outlier test.
## 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
## 7 -0.18 0.22 0.09 -0.03 -0.45_* 0.85_* 0.05 0.02
## 20 0.01 -0.02 0.04 0.00 0.04 1.07_* 0.00 0.05
## 21 0.05 0.08 -0.08 0.04 0.11 1.07_* 0.00 0.05
## 25 -0.14 0.06 0.06 0.08 -0.23 0.93_* 0.01 0.01
## 35 0.00 -0.01 0.01 0.00 -0.01 1.06_* 0.00 0.04
## 44 -0.03 -0.03 0.02 0.07 -0.08 1.06_* 0.00 0.04
## 69 0.28 -0.26 -0.16 0.17 0.64_* 0.69_* 0.09 0.02
## 82 0.29 -0.83 0.64 -0.05 0.95_* 0.55_* 0.19 0.03
## 100 0.02 -0.03 0.06 0.07 0.09 1.10_* 0.00 0.08_*
## 106 -0.01 -0.02 0.03 0.00 -0.04 1.07_* 0.00 0.04
## 152 -0.02 0.08 -0.05 -0.04 -0.09 1.07_* 0.00 0.05
## 155 0.03 -0.05 -0.02 0.09 0.13 1.09_* 0.00 0.07_*
## 180 -0.03 -0.04 0.06 -0.02 -0.07 1.09_* 0.00 0.06_*
## 188 -0.02 0.05 -0.09 0.00 -0.10 1.07_* 0.00 0.05
## 190 -0.19 -0.20 -0.10 0.32 -0.53_* 0.96 0.07 0.05
## 192 0.37 0.21 -0.19 0.25 0.48_* 0.67_* 0.05 0.01# 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
## 82 82 5.871112 0.02546700 0.19235005 TRUE FALSE TRUE TRUE
## 192 192 4.714648 0.01024445 0.05189670 TRUE FALSE TRUE TRUE
## 69 69 4.569573 0.01899437 0.09176694 TRUE FALSE TRUE TRUE
## 7 7 -3.225228 0.01892940 0.04787938 TRUE FALSE TRUE TRUE
## 190 190 -2.314024 0.04956792 0.06829860 FALSE TRUE TRUE TRUE
## 122 122 2.080521 0.03502345 0.03862006 FALSE FALSE TRUE TRUE
## 147 147 -1.906798 0.02851264 0.02632377 FALSE FALSE TRUE TRUE
## 18 18 -1.822566 0.02571029 0.02165765 FALSE FALSE TRUE TRUE
## 197 197 -1.726085 0.04077432 0.03134486 FALSE TRUE TRUE TRUE
## 47 47 1.725732 0.05615937 0.04385802 FALSE TRUE TRUE TRUEBonferroni outlier test (car):
## rstudent unadjusted p-value Bonferroni p
## 82 5.871112 1.8362e-08 3.6724e-06
## 192 4.714648 4.6016e-06 9.2032e-04
## 69 4.569573 8.6590e-06 1.7318e-03Plot cepat untuk Cook’s distance:
plot(cookd, type="h", main="Cook's Distance", ylab="D", xlab="Observasi")
abline(h = thr_cook, lty=2)
Tindak lanjut (opsional): Coba refit tanpa observasi berpengaruh ekstrem
dan bandingkan koefisien/fit. Pertimbangkan robust regression
(M-estimator) jika banyak outlier: MASS::rlm.
# Refit tanpa titik yang sangat berpengaruh (contoh)
to_drop <- which(cookd > thr_cook | abs(studres) > thr_res)
mod_refit <- lm(y ~ X1 + X2 + X3, data = dat[-to_drop, ])
broom::tidy(mod0)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 2.04 0.0774 26.4 2.22e-66
## 2 X1 1.40 0.109 12.9 7.42e-28
## 3 X2 -0.793 0.0973 -8.15 4.23e-14
## 4 X3 0.509 0.0773 6.58 4.22e-10## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 2.01 0.0560 35.9 5.72e-85
## 2 X1 1.49 0.0806 18.4 1.09e-43
## 3 X2 -0.847 0.0704 -12.0 5.51e-25
## 4 X3 0.430 0.0572 7.51 2.52e-127. Ringkasan & Rekomendasi
• Cek signifikansi koefisien (t-test) dan kecocokan model (F-test, R²).
• Perhatikan multikolinearitas (VIF); bila tinggi, pertimbangkan transformasi, pengurangan variabel, atau ridge/lasso.
• Jika heteroskedastisitas terdeteksi, gunakan SE robust untuk inferensi atau model varian (mis. WLS).
• Lakukan seleksi model (AIC/BIC, CV) untuk keseimbangan bias–varian.
• Evaluasi outlier/pengaruh; pertimbangkan robust regression bila perlu.
• Validasi hasil secara substantif (masuk akal secara domain) dan replikasi (seed, pipeline).
1.13 Tugas/Praktikum dan Bank Soal
1.13.1 Praktikum
1.Diagnostik lengkap OLS pada data mtcars: heteroskedastisitas (Breusch–Pagan), multikolinearitas (VIF), pengaruh (Cook’s D), normalitas (QQ-plot). Laporkan temuan & rekomendasi.
jawab
Analisis ini menggunakan dataset mtcars untuk memodelkan konsumsi bahan bakar (mpg) dengan semua variabel lain sebagai prediktor.
# Load libraries
library(car)
library(lmtest)
library(sandwich)
library(ggplot2)
library(broom)
library(MASS)
library(glmnet)
library(corrplot)
# Data & model
data(mtcars)
mod <- lm(mpg ~ ., data = mtcars)
summary(mod)##
## Call:
## lm(formula = mpg ~ ., data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4506 -1.6044 -0.1196 1.2193 4.6271
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.30337 18.71788 0.657 0.5181
## cyl -0.11144 1.04502 -0.107 0.9161
## disp 0.01334 0.01786 0.747 0.4635
## hp -0.02148 0.02177 -0.987 0.3350
## drat 0.78711 1.63537 0.481 0.6353
## wt -3.71530 1.89441 -1.961 0.0633 .
## qsec 0.82104 0.73084 1.123 0.2739
## vs 0.31776 2.10451 0.151 0.8814
## am 2.52023 2.05665 1.225 0.2340
## gear 0.65541 1.49326 0.439 0.6652
## carb -0.19942 0.82875 -0.241 0.8122
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.65 on 21 degrees of freedom
## Multiple R-squared: 0.869, Adjusted R-squared: 0.8066
## F-statistic: 13.93 on 10 and 21 DF, p-value: 3.793e-07##
## studentized Breusch-Pagan test
##
## data: mod
## BP = 14.914, df = 10, p-value = 0.1352Interpretasi:
H0: Varians residual homogen (homoskedastisitas)
H1: Varians residual tidak homogen (heteroskedastisitas)
Jika p-value < 0.05, terdapat heteroskedastisitas. Jika heteroskedastisitas terdeteksi, gunakan robust standard errors.
## cyl disp hp drat wt qsec vs am
## 15.373833 21.620241 9.832037 3.374620 15.164887 7.527958 4.965873 4.648487
## gear carb
## 5.357452 7.908747Interpretasi
VIF > 5 atau 10 menunjukkan multikolinearitas tinggi.
Jika ada variabel dengan VIF tinggi, pertimbangkan untuk menghapus atau menggabungkan variabel tersebut.
# Pengaruh Observasi (Cook's Distance)
cooksd <- cooks.distance(model)
plot(cooksd, type = "h", main = "Cook's Distance", ylab = "Cook's distance")
abline(h = 4/(nrow(mtcars)-length(model$coefficients)), col = "red", lty = 2)
Interpretasi
Observasi dengan Cook’s D > threshold (4/(n-k-1)) dianggap berpengaruh tinggi. Periksa observasi tersebut dan pertimbangkan analisis lebih lanjut.
# Normalitas Residual (QQ Plot)
qqnorm(residuals(mod), main = "QQ-Plot Residuals")
qqline(residuals(mod), col = "red")
Interpretasi
-Titik-titik yang mengikuti garis lurus menunjukkan residual mendekati normal.
-Penyimpangan di ekor dapat menunjukkan pelanggaran asumsi normalitas.
2.Robust SE (HC3) vs WLS pada data simulasi heteroskedastik, bandingkan standar error dan p-value
jawab
#Simulasi Data Heteroskedastik
set.seed(123)
n <- 100
x <- runif(n, 0, 10)
# Varians error meningkat dengan x (heteroskedastisitas)
sigma <- 0.5 + 0.5 * x
eps <- rnorm(n, 0, sigma)
y <- 2 + 3 * x + eps
sim_data <- data.frame(x, y)##
## Call:
## lm(formula = y ~ x, data = sim_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.9096 -1.7710 0.0479 1.3846 11.1165
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9238 0.6514 2.953 0.00393 **
## x 2.9758 0.1136 26.202 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.22 on 98 degrees of freedom
## Multiple R-squared: 0.8751, Adjusted R-squared: 0.8738
## F-statistic: 686.6 on 1 and 98 DF, p-value: < 2.2e-16library(sandwich)
library(lmtest)
robust_se <- coeftest(model_sim, vcov = vcovHC(model_sim, type = "HC3"))
robust_se##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.92380 0.45972 4.1847 6.231e-05 ***
## x 2.97578 0.12682 23.4652 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Model WLS (Weighted Least Squares)
weights <- 1 / (sigma^2)
model_wls <- lm(y ~ x, data = sim_data, weights = weights)
summary(model_wls)##
## Call:
## lm(formula = y ~ x, data = sim_data, weights = weights)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -2.24077 -0.61176 -0.00507 0.59445 2.23207
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.02519 0.22423 9.032 1.53e-14 ***
## x 2.95723 0.08258 35.810 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9694 on 98 degrees of freedom
## Multiple R-squared: 0.929, Adjusted R-squared: 0.9283
## F-statistic: 1282 on 1 and 98 DF, p-value: < 2.2e-16Interpretasi
- Robust SE dan WLS memberikan standar error yang lebih realistis dibandingkan OLS biasa saat heteroskedastisitas ada.
- P-value pada Robust SE dan WLS biasanya lebih konservatif (lebih besar) dibanding OLS biasa.
3.Spline vs polinomial: modelkan mpg ~ wt dan bandingkan AIC/BIC + 5-fold CV
jawab
library(splines)
library(caret)
# Model polinomial
poly_mod <- lm(mpg ~ poly(wt, 3), data=mtcars)
# Model spline
spline_mod <- lm(mpg ~ bs(wt, df=4), data=mtcars)
# AIC & BIC
AIC(poly_mod, spline_mod)## df AIC
## poly_mod 5 160.0365
## spline_mod 6 161.7642## df BIC
## poly_mod 5 167.3652
## spline_mod 6 170.5586# Cross-validation 5-fold
train_control <- trainControl(method="cv", number=5)
poly_cv <- train(mpg ~ poly(wt, 3), data=mtcars, method="lm",
trControl=train_control)
spline_cv <- train(mpg ~ bs(wt, df=4), data=mtcars, method="lm",
trControl=train_control)
poly_cv## Linear Regression
##
## 32 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 24, 27, 26, 26, 25
## Resampling results:
##
## RMSE Rsquared MAE
## 3.170709 0.8688092 2.729868
##
## Tuning parameter 'intercept' was held constant at a value of TRUE## Linear Regression
##
## 32 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 26, 26, 24, 26, 26
## Resampling results:
##
## RMSE Rsquared MAE
## 2.962085 0.7361487 2.471983
##
## Tuning parameter 'intercept' was held constant at a value of TRUEInterpretasi:
-Model dengan AIC/BIC lebih kecil lebih baik.
-CV error yang lebih kecil lebih baik.
4.AR(1) GLS: gunakan kode GLS di atas, ganti x dengan sinyal musiman (mis. sin(2pit/12)), interpretasikan ; bandingkan SE OLS vs GLS.
jawab
library(nlme)
# Buat sinyal musiman
t <- 1:100
x <- sin(2*pi*t/12)
y <- 1 + 0.5*x + arima.sim(list(ar=0.5), n=100)
# OLS
ols_mod <- lm(y ~ x)
# GLS dengan AR(1)
gls_mod <- gls(y ~ x, correlation = corAR1(form = ~1))
summary(ols_mod)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3006 -0.7296 -0.0317 0.5553 3.2250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.5483 0.1106 14.00 < 2e-16 ***
## x 0.5651 0.1553 3.64 0.000438 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.105 on 98 degrees of freedom
## Multiple R-squared: 0.1191, Adjusted R-squared: 0.1101
## F-statistic: 13.25 on 1 and 98 DF, p-value: 0.000438## Generalized least squares fit by REML
## Model: y ~ x
## Data: NULL
## AIC BIC logLik
## 308.9965 319.3363 -150.4982
##
## Correlation Structure: AR(1)
## Formula: ~1
## Parameter estimate(s):
## Phi
## 0.2298899
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.5449714 0.1400980 11.027794 0.0000
## x 0.5623764 0.1876046 2.997668 0.0034
##
## Correlation:
## (Intr)
## x -0.048
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.06733415 -0.65611450 -0.02345069 0.50390198 2.90474978
##
## Residual standard error: 1.110599
## Degrees of freedom: 100 total; 98 residualInterpretasi:
-Bandingkan estimasi koefisien x dari OLS vs GLS.
-GLS memperhitungkan autokorelasi residual, sehingga SE lebih akurat.
Kesimpulan
-Diagnostik OLS menunjukkan apakah ada heteroskedastisitas, multikolinearitas, observasi berpengaruh, atau non-normalitas residual.
-Robust SE (HC3) dan WLS sama-sama solusi heteroskedastisitas; WLS lebih efisien jika bobot tepat.
-Perbandingan spline vs polinomial menunjukkan trade-off antara fleksibilitas dan kompleksitas.
-GLS lebih baik dibanding OLS pada data dengan autokorelasi.
1.13.2 Soal Latihan
Soal 1
Mengapa OLS tetap unbiased di bawah heteroskedastisitas, namun tidak efisien? Bagaimana robust SE mengatasi masalah inferensi?
Penjelasan
Unbiasedness (tak bias): OLS estimator ditulis sebagai:
\(\hat{\beta}_{OLS} = (X^{\top}X)^{-1}X^{\top}y = (X^{\top}X)^{-1}X^{\top}(X\beta + \varepsilon) = \beta + (X^{\top}X)^{-1}X^{\top}\varepsilon.\)
Dengan asumsi klasik bahwa \(E[\varepsilon\mid X]=0\), didapat:
\(E[\hat{\beta}_{OLS}\mid X] = \beta + (X^{\top}X)^{-1}X^{\top}E[\varepsilon\mid X] = \beta.\)
Perhatikan bahwa tidak diperlukan asumsi homoskedastisitas — cukup mean nol. Oleh sebab itu, OLS tetap unbiased walau residu heteroskedastik.
Efisiensi (BLUE): Menurut teorema Gauss–Markov, OLS adalah BLUE jika \(Var(\varepsilon|X) = \sigma^2 I\). Jika heteroskedastik, varians OLS menjadi:
\(Var(\hat{\beta}_{OLS}\mid X) = (X^{\top}X)^{-1}X^{\top}\Omega X (X^{\top}X)^{-1},\)
yang umumnya lebih besar dari \(\sigma^2 (X^{\top}X)^{-1}\). Jadi OLS tidak lagi efisien.
Robust SE (sandwich estimator)
Estimator kovarian robust heteroskedastisitas adalah:
\(\widehat{Var}(\hat{\beta}) = (X^{\top}X)^{-1}\Bigg(\sum_{i=1}^n x_i x_i^{\top}\hat{u}_i^2\Bigg)(X^{\top}X)^{-1}.\)
Varian HC3 khususnya:
\(\hat{u}_i^2/(1-h_{ii})^2,\)
sehingga lebih konservatif untuk sampel kecil. Robust SE menjaga koefisien OLS tak berubah, hanya menghitung ulang standar error agar inferensi valid.
Soal 2
Turunkan estimator WLS saat varian residual \(\propto |x_i|\)
Penjelasan
Jika \(Var(\varepsilon_i) = \sigma^2 |x_i|\), maka bobot WLS:
\(w_i = \frac{1}{|x_i|}.\)
Estimator WLS adalah:
\(\hat{\beta}_{WLS} = (X^{\top}WX)^{-1}X^{\top}Wy,\)
dengan \(W = diag(1/|x_1|, 1/|x_2|, \dots, 1/|x_n|).\)
Soal 3
Bandingkan Ridge dan Lasso dari sudut pandang bias–variance dan seleksi variabel.
Penjelasan
Ridge
\(\hat{\beta}^{ridge} = \arg\min_{\beta}\left\{\sum_{i=1}^n (y_i - x_i^{\top}\beta)^2 + \lambda \sum_{j=1}^p \beta_j^2\right\}.\)
Lasso
\(\hat{\beta}^{lasso} = \arg\min_{\beta}\left\{\sum_{i=1}^n (y_i - x_i^{\top}\beta)^2 + \lambda \sum_{j=1}^p |\beta_j|\right\}.\)
- Ridge: shrinkage halus, tidak membuat koefisien persis nol.
- Lasso: menghasilkan sparse solution (beberapa \(\beta_j=0\)) → seleksi variabel.
- Bias–variance trade-off: keduanya menambah bias, mengurangi varians.
Soal 4
Hitungan: Diberikan matriks 𝑋 dan vektor 𝑦, hitung ̂ 𝛽 OLS serta 𝐻; identifikasi leverage maksimum.
Penjelasan
Contoh dengan data sederhana:
X <- cbind(1, c(1,2,3,4,5,20)) # intercept + 1 regressor
Y <- c(2,3,4,5,6,30)
beta_hat <- solve(t(X)%*%X) %*% t(X) %*% Y
H <- X %*% solve(t(X)%*%X) %*% t(X)
leverages <- diag(H)
list(beta_hat=beta_hat, H=H, leverage=leverages, max_leverage=max(leverages))## $beta_hat
## [,1]
## [1,] -0.4651163
## [2,] 1.5083056
##
## $H
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.2598007 0.24053156 0.221262458 0.20199336 0.1827243 -0.106312292
## [2,] 0.2405316 0.22524917 0.209966777 0.19468439 0.1794020 -0.049833887
## [3,] 0.2212625 0.20996678 0.198671096 0.18737542 0.1760797 0.006644518
## [4,] 0.2019934 0.19468439 0.187375415 0.18006645 0.1727575 0.063122924
## [5,] 0.1827243 0.17940199 0.176079734 0.17275748 0.1694352 0.119601329
## [6,] -0.1063123 -0.04983389 0.006644518 0.06312292 0.1196013 0.966777409
##
## $leverage
## [1] 0.2598007 0.2252492 0.1986711 0.1800664 0.1694352 0.9667774
##
## $max_leverage
## [1] 0.9667774Secara umum:
- Estimator OLS: \(\hat{\beta}_{OLS} = (X^{\top}X)^{-1}X^{\top}y.\)
- Matriks hat: \(H = X(X^{\top}X)^{-1}X^{\top}.\)
- Leverage: diagonal dari H (\(h_{ii}\)). Rata-rata leverage = \(p/n\).
1.13.3 Pilihan Ganda
(A) VIF tinggi → multikolinearitas\(VIF_j = \frac{1}{1-R_j^2}\).
(B) Cook’s D → mendeteksi observasi berpengaruh.
(C) Dalam AR(1): \(Cov(\varepsilon_t, \varepsilon_{t-k}) = \sigma_{\varepsilon}^2 \rho^k,\) dengan \(\sigma_{\varepsilon}^2 = \frac{\sigma_u^2}{1-\rho^2}.\)
Ringkasan
OLS tetap unbiased di bawah heteroskedastisitas, tapi tidak efisien. Robust SE memperbaiki inferensi.
WLS optimal jika bobot proporsional terhadap varians residual.
Ridge & Lasso: trade-off bias–variance, Lasso juga melakukan seleksi variabel.
Leverage dan Cook’s D penting untuk mendeteksi observasi berpengaruh.
2 Regresi Nonlinier: Model, Estimasi, dan Aplikasi
2.1 Tujuan Pembelajaran Analisis Regresi Nonlinear
- Memperkenalkan berbagai permasalahan yang relevan dengan pemodelan dan pemfittingan fungsi regresi nonlinear.
- Menyajikan representasi grafis untuk menilai kualitas dari interval kepercayaan aproksimasi.
- Mengenalkan beberapa bagian dari perangkat lunak statistik R yang dapat membantu dalam menyelesaikan masalah konkret.
2.2 Motivasi
Regresi linear mengasumsikan hubungan antara variabel respon 𝑌 dan
prediktor 𝑋 berbentuk linear:
\[Y = \beta_0 + \beta_1 X + \varepsilon,\
\varepsilon \ \text{iid} \sim N(0,\ \sigma^2_{\varepsilon})\]
Namun, dalam banyak fenomena nyata (biologi, farmasi, pertumbuhan,
epidemiologi), hubungan tidak linear.
Contoh: ### Pertumbuhan bakteri mengikuti kurva logistik.
Model
Model pertumbuhan bakteri seringkali mengikuti bentuk kurva logistik. Model matematisnya dapat dituliskan sebagai:
\[Y (t) = \dfrac{K}{1 + \exp\{-(\alpha + \beta t)\}} + \varepsilon,\ \varepsilon \sim N(0,\ \sigma^2)\] Keterangan:
- \(Y (t)\) : jumlah populasi bakteri pada waktu \(t\)
- \(t\) : waktu pengamatan (misalnya jam, hari)
- \(K\) : daya dukung (carrying capacity), yaitu ukuran maksimum populasi bakteri yang dapat dicapai dalam kondisi lingkungan tertentu
- \(\alpha\) : konstanta yang berhubungan dengan kondisi awal (menentukan titik tengah pertumbuhan)
- \(\beta\) : laju pertumbuhan bakteri (growth rate), semakin besar nilainya semakin cepat bakteri berkembang
- \(\exp\{-(\alpha + \beta t)\}\) : komponen eksponensial yang mengatur bentuk kurva pertumbuhan logistik
- \(\varepsilon\) : error atau deviasi acak dari model ideal, diasumsikan berdistribusi normal
- \(\varepsilon \sim N(0,\ \sigma^2)\) : menyatakan bahwa error memiliki rata-rata nol dan varians \(\sigma^2\)
set.seed(42)
library(ggplot2)
library(dplyr)
library(tidyr)
# Parameter
K <- 1.0 # kapasitas maksimum (misal proporsi)
alpha <--3.0 # intercept logit
beta <- 0.9 # laju pertumbuhan
sigma <- 0.03 # noise
# Data
t <- seq(0, 8, length.out = 120)
mu <- K/(1 + exp(-(alpha + beta*t)))
y <- mu + rnorm(length(t), sd = sigma)
logistic_df <- tibble(t = t, y = y, mu = mu)
# Plot
ggplot(logistic_df, aes(t, y)) +
geom_point(alpha = 0.6) +
geom_line(aes(y = mu), linewidth = 1) +
labs(title = "Pertumbuhan Bakteri (Kurva Logistik)",
x = "Waktu (t)", y = "Ukuran/Proporsi Populasi") +
theme_minimal()
Interpretasi:
- Pada awal waktu (\(t\) kecil), populasi \(Y (t)\) masih relatif rendah.
- Seiring waktu, laju pertumbuhan meningkat pesat hingga mendekati maksimum (\(K\)).
- Setelah mendekati kapasitas lingkungan, pertumbuhan melambat dan akhirnya stabil di sekitar \(K\).
2.2.1 Respon obat mengikuti fungsi Michaelis-Menten.
Model
Model farmakokinetik Michaelis–Menten sering digunakan untuk menjelaskan
hubungan antara konsentrasi obat (𝑋) dan respon biologis (𝑌). Modelnya
adalah:
\[Y = \dfrac{V_{\text{max}} X}{K_m + X} + \varepsilon,\ \varepsilon \sim N(0,\ \sigma^2)\]
Keterangan:
- \(Y\) : respon obat (misalnya, laju reaksi enzim atau efek farmakologis yang diamati)
- \(X\) : konsentrasi obat atau substrat
- \(V_{\max}\) : respon maksimum (laju reaksi atau efek maksimum yang dapat dicapai bila semua reseptor/enzim sudah jenuh)
- \(K_m\) : konstanta Michaelis–Menten, yaitu konsentrasi \(X\) saat respon mencapai setengah dari \(V_{\max}\)
- \(\varepsilon\) : error atau deviasi acak dari model, diasumsikan berdistribusi normal
- \(\varepsilon \sim N(0,\ \sigma^2)\) : menyatakan bahwa error memiliki rata-rata nol dan varians \(\sigma^2\)
# Parameter
Vmax <- 2.0
Km <- 1.2
sigma <- 0.06
# Data
X <- seq(0, 6, length.out = 80)
mu <- (Vmax * X) / (Km + X)
Y <- mu + rnorm(length(X), sd = sigma)
mm_df <- tibble(X = X, Y = Y, mu = mu)
ggplot(mm_df, aes(X, Y)) +
geom_point(alpha = 0.6) +
geom_line(aes(y = mu), linewidth = 1) +
labs(
title = "Respon Obat (Michaelis–Menten)",
x = "Konsentrasi Obat (X)",
y = "Respon (Y)"
) +
theme_minimal()
Interpretasi:
- Pada konsentrasi obat rendah (\(X \ll K_m\)), respon \(Y\) meningkat hampir linear terhadap \(X\).
- Saat konsentrasi obat mendekati \(K_m\), respon mulai melambat, menandakan jenuh parsial.
- Pada konsentrasi obat sangat tinggi (\(X \gg K_m\)), respon mendekati \(V_{\max}\), sehingga tidak ada peningkatan berarti meskipun dosis obat ditambah.
2.2.2 Ekonomi: diminishing return mengikuti bentuk kurva nonlinear.
Model Contoh model sederhana dengan bentuk konkaf dapat dituliskan sebagai:
\[Y = A X^{\alpha} + \varepsilon,\ \ 0 < \alpha < 1,\ \ \varepsilon \sim N(0,\ \sigma^2)\]
Alternatif lain:
\[Y = a \log(1 + bX) + \varepsilon,\ \ \varepsilon \sim N(0,\ \sigma^2)\]
Keterangan:
- \(Y\) : output ekonomi (misalnya
hasil produksi, keuntungan, atau utilitas)
- \(X\) : input (misalnya modal,
tenaga kerja, atau konsumsi)
- \(A\) : konstanta skala (menentukan
besarnya output secara keseluruhan)
- \(\alpha\) : parameter elastisitas
(dengan \(0 < \alpha < 1\) untuk
menunjukkan sifat diminishing return)
- \(a, b\) : parameter dalam model
logaritmik yang mengatur skala (\(a\))
dan sensitivitas input (\(b\))
- \(\varepsilon\) : error atau
deviasi acak dari model ideal
- \(\varepsilon \sim N(0,\ \sigma^2)\) : error diasumsikan berdistribusi normal dengan rata-rata nol dan varians \(\sigma^2\)
# Parameter
A <- 5
alpha <- 0.5 # < 1 → diminishing marginal return
sigma <- 0.3
# Data
X <- seq(0, 100, length.out = 120)
mu <- A * X^alpha
Y <- mu + rnorm(length(X), sd = sigma)
dr_df <- tibble(X = X, Y = Y, mu = mu)
ggplot(dr_df, aes(X, Y)) +
geom_point(alpha = 0.6) +
geom_line(aes(y = mu), linewidth = 1) +
labs(title = "Diminishing Return: Y = A * X^{alpha}, 0<alpha<1",
x = "Input/Modal (X)", y = "Output (Y)") +
theme_minimal()
Interpretasi:
- Pada model \(Y = A X^{\alpha}\), karena \(\alpha < 1\), tambahan output dari setiap unit tambahan input semakin menurun. Contohnya, tambahan 1 pekerja pertama menambah output banyak, tapi tambahan pekerja ke-100 hanya menambah sedikit output.
- Pada model \(Y = a \log(1 + bX)\), kenaikan output juga berbentuk konkaf: besar pada awalnya (saat \(X\) kecil) lalu menurun dan mendekati asimtot saat \(X\) besar.
Merujuk pada contoh diatas maka dibutuhkan regresi nonlinear.
2.3 Model Regresi Nonlinear
Secara umum, model regresi nonlinear ditulis sebagai:
\[Y_i = f(X_i, \theta) + \varepsilon_i,\ \ \varepsilon_i \sim N(0,\ \sigma^2)\]
- \(f(\cdot)\) adalah fungsi
nonlinear terhadap parameter \(\theta\).
- Estimasi parameter biasanya menggunakan Nonlinear Least Squares (NLS).
Contoh model eksponensial:
\[Y = \beta_0 e^{\beta_1 X} + \varepsilon\]
Catatan: Model Eksponensial
Model eksponensial dapat dituliskan sebagai: \[Y = \beta_0 e^{\beta_1 X} + \varepsilon\]
Model ini pada dasarnya nonlinear dalam parameter jika ditulis langsung. Namun, jika asumsi error diubah (misalnya error bersifat multiplikatif, bukan aditif), maka model dapat dipandang sebagai linear setelah transformasi logaritma.
Dua Pendekatan Umum
1. Error Aditif
\[Y = \beta_0 e^{\beta_1 X} + \varepsilon,\ \varepsilon \sim N(0,\ \sigma^2)\] - Model ini benar-benar nonlinear terhadap parameternya. - Estimasi memerlukan metode Nonlinear Least Squares (NLS).
2. Error Multiplikatif
\[Y = \beta_0 e^{\beta_1 X} \cdot \varepsilon,\ \varepsilon > 0\]
Ambil logaritma kedua sisi:
\[\ln Y = \ln \beta_0 + \beta_1 X + \ln \varepsilon\]
Jika diasumsikan:
\[\ln \varepsilon \sim N(0,\ \tau^2),\] maka model menjadi linear dalam parameter:
\[\ln Y = \alpha_0 + \beta_1 X + u,\ u \sim N(0,\ \sigma^2)\] dengan definisi:
\[\alpha_0 = \ln \beta_0\] Estimasi
- Model dapat diestimasi menggunakan regresi linear biasa (OLS).
- Setelah estimasi, kita bisa kembali ke skala asli:
\[\hat{\beta}_0 = e^{\hat{\alpha}_0}\] Ringkasan
- Error aditif → model nonlinear → estimasi dengan NLS.
- Error multiplikatif → bisa ditransformasi log → estimasi dengan OLS.
2.4 Estimasi & Inferensi
2.4.1 Estimasi parameter menggunakan metode iteratif (Gauss-Newton,Levenberg-Marquardt).
Kita ingin menaksir parameter \(\theta \in \mathbb{R}^p\) pada model
\[Y_i = f(X_i, \theta) + \varepsilon_i,\ \varepsilon_i \ \text{iid} \sim \mathcal{N}(0, \sigma^2),\ i = 1, \ldots, n.\] Estimator nonlinear least squares (NLS) meminimalkan jumlah kuadrat residu:
\[S(\theta) = \sum_{i=1}^{n} r_i(\theta) = y_i - f(X_i, \theta).\] Notasi Jacobian dari \(f\) terhadap parameter: \[J_{ij}(\theta) = \dfrac{\partial f(x_i, \theta)}{\partial \theta_j}, \quad J(\theta) \in \mathbb{R}^{n\times p}.\] Di sekitar tebakan \(\theta^{(k)}\), pendekatan Taylor orde-1 memberi: \[r(\theta^{(k)} + \Delta) \approx r^{(k)} - J^{(k)} \Delta,\ \ r^{(k)} = y- f(\theta^{(k)}).\]
Gauss–Newton (GN) — Ide Inti
Dengan linearisasi di atas, langkah perbaikan \(\Delta\) dicari dari normal equations: \[(J^\top J) \ \Delta = J^\top r.\] Lalu diperbarui \(\theta^{(k+1)} = \theta^{(k)} + \Delta.\)
Ulangi hingga konvergen (mis. \(\lVert \Delta \rVert\) kecil atau \(|S^{(k+1)} - S^{(k)}|\) kecil).
Levenberg–Marquardt (LM) — Ide Inti (Damped GN / Trust-Region)
LM menambahkan redaman agar lebih stabil ketika \(J^\top J\) buruk-kondisi / tebakan awal jauh:
\[(J^\top J + \lambda D) \ \Delta = J^\top r,\] dengan \(\lambda > 0\) (faktor redaman) dan \(D = I\) atau \(D = \operatorname{diag}(J^\top J).\)
Aturan sederhana:
Jika langkah \(\Delta\) menurunkan \(S(\theta)\), kecilkan \(\lambda\) (lebih agresif).
Jika tidak, besarkan \(\lambda\) (lebih konservatif) dan coba lagi. Kasus Spesifik: Model Eksponensial
Kasus Spesifik: Model Eksponensial Kita gunakan model:
\[Y = \beta_0 e^{\beta_1 X} + \varepsilon\], \[\varepsilon \sim \mathcal{N}(0, \sigma^2)\].
Untuk \(x_i\), \(f(x_i, \theta) = \beta_0 e^{\beta_1 x_i}\) , \(\theta = (\beta_0, \beta_1)^\top\) .
Jacobian \(J\) (terhadap \(\beta_0, \beta_1\)):
\[\dfrac{\partial f}{\partial \beta_0} \dfrac{\partial f}{\partial \beta_1} = e^{\beta_1 x_i} , = \beta_0 x_i e^{\beta_1 x_i} . \]
Simulasi Data & Visualisasi Awal
set.seed(123)
library(ggplot2)
library(dplyr)
library(tidyr)
beta0_true <- 2.0
beta1_true <- 0.3
sigma_true <- 0.4
# Small dataset agar mudah ditelusuri langkah demi langkah
x <- c(0, 1, 2, 3, 4, 5)
f_true <- beta0_true * exp(beta1_true * x)
y <- f_true + rnorm(length(x), sd = sigma_true)
dat <- tibble(x, y, f_true)
dat## # A tibble: 6 × 3
## x y f_true
## <dbl> <dbl> <dbl>
## 1 0 1.78 2
## 2 1 2.61 2.70
## 3 2 4.27 3.64
## 4 3 4.95 4.92
## 5 4 6.69 6.64
## 6 5 9.65 8.96Plot
ggplot(dat, aes(x, y)) +
geom_point(size=2) +
geom_line(aes(y = f_true), linewidth = 0.9) +
labs(title="Data Eksponensial (dengan noise) & Kurva Sejati",
x="X", y="Y") +
theme_minimal()
Gauss–Newton: Estimasi Manual (Step-by-step) Fungsi Bantu (f, residu, Jacobian)
f_exp <- function(x, theta){ # theta = c(beta0, beta1)
beta0 <- theta[1]; beta1 <- theta[2]
beta0 * exp(beta1 * x)
resid_vec <- function(x, y, theta){
y- f_exp(x, theta)
}
}
J_exp <- function(x, theta){
beta0 <- theta[1]; beta1 <- theta[2]
e <- exp(beta1 * x)
cbind( d_beta0 = e,
d_beta1 = beta0 * x * e )
}Satu langkah iterasi Gauss–Newton (GN) diperoleh dengan menyelesaikan persamaan normal:
\[(J^\top J)\ \Delta = J^\top r\]
# --- Definisi model eksponensial: Y = beta0 * exp(beta1 * X) + e
f_exp <- function(x, theta) theta[1] * exp(theta[2] * x)
J_exp <- function(x, theta){
e <- exp(theta[2] * x)
cbind(d_beta0 = e, d_beta1 = theta[1] * x * e)
}
resid_vec <- function(x, y, theta) y - f_exp(x, theta)
# Jika x/y belum ada, buat data contoh (opsional, aman di-skip jika sudah ada)
if(!exists("x") || !exists("y")){
set.seed(1); n <- 120; x <- runif(n, 0, 3)
beta_true <- c(2, 0.8); y <- beta_true[1]*exp(beta_true[2]*x) + rnorm(n, sd = 0.5)
}
# --- Langkah GN
gn_step <- function(x, y, theta){
r <- resid_vec(x, y, theta)
J <- J_exp(x, theta)
lhs <- crossprod(J, J) # J^T J
rhs <- crossprod(J, r) # J^T r
Delta <- solve(lhs, rhs) # Δ
list(Delta = as.numeric(Delta),
theta_new = theta + as.numeric(Delta),
r = r, J = J, lhs = lhs, rhs = rhs,
SSE = sum(r^2))
}
theta0 <- c(beta0 = 1.0, beta1 = 0.0)
step1 <- gn_step(x, y, theta0)
step1## $Delta
## [1] 0.2542352 1.4943019
##
## $theta_new
## beta0 beta1
## 1.254235 1.494302
##
## $r
## [1] 0.7758097 1.6076466 3.2677209 3.9474096 5.6919489 8.6494041
##
## $J
## d_beta0 d_beta1
## [1,] 1 0
## [2,] 1 1
## [3,] 1 2
## [4,] 1 3
## [5,] 1 4
## [6,] 1 5
##
## $lhs
## d_beta0 d_beta1
## d_beta0 6 15
## d_beta1 15 55
##
## $rhs
## [,1]
## d_beta0 23.93994
## d_beta1 86.00013
##
## $SSE
## [1] 136.6569Iterasi Beberapa Kali Hingga Konvergen
gn_fit <- function(x, y, theta_init, max_it = 20, tol = 1e-8) {
th <- theta_init
hist <- data.frame(
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 <- rbind(hist, data.frame(
iter = k,
beta0 = th[1],
beta1 = th[2],
SSE = sum(resid_vec(x, y, th)^2),
step_norm = step_norm
))
if (!is.na(step_norm) && step_norm < tol) break
}
list(theta = th, history = hist)
}
gn_res <- gn_fit(x, y, theta_init = theta0, max_it = 20)
# tampilkan aman di knit
knitr::kable(gn_res$history)| iter | beta0 | beta1 | SSE | step_norm | |
|---|---|---|---|---|---|
| beta0 | 0 | 1.0000000 | 0.0000000 | 1.366569e+02 | NA |
| beta01 | 1 | 1.2542352 | 1.4943019 | 5.064963e+06 | 1.5157750 |
| beta02 | 2 | 0.0812940 | 1.4821976 | 1.616697e+04 | 1.1730036 |
| beta03 | 3 | 0.0851451 | 1.2867953 | 1.962659e+03 | 0.1954403 |
| beta04 | 4 | 0.1776123 | 0.9050502 | 6.804454e+01 | 0.3927843 |
| beta05 | 5 | 0.6642144 | 0.2703978 | 9.850920e+01 | 0.7997283 |
| beta06 | 6 | 1.9594441 | 0.4222512 | 6.231191e+01 | 1.3041009 |
| beta07 | 7 | 1.8882580 | 0.3468196 | 2.261074e+00 | 0.1037178 |
| beta08 | 8 | 1.9672064 | 0.3166507 | 4.683862e-01 | 0.0845164 |
| beta09 | 9 | 1.9752973 | 0.3148120 | 4.653472e-01 | 0.0082971 |
| beta010 | 10 | 1.9753264 | 0.3148085 | 4.653472e-01 | 0.0000293 |
| beta011 | 11 | 1.9753264 | 0.3148085 | 4.653472e-01 | 0.0000000 |
## beta0 beta1
## 1.9753264 0.3148085Kurva Hasil GN
dat$yhat_gn <- f_exp(dat$x, gn_res$theta)
ggplot(dat, aes(x, y)) +
geom_point(size = 2, color = "steelblue") +
geom_line(aes(y = yhat_gn), color = "darkred", linewidth = 1) +
stat_ellipse(level = 0.95, color = "blue") +
labs(
title = "Hasil Estimasi Gauss–Newton",
subtitle = paste0(
"beta0 = ", round(gn_res$theta[1], 3),
", beta1 = ", round(gn_res$theta[2], 3)
),
x = "X", y = "Y"
) +
theme_minimal()
Levenberg–Marquardt: Estimasi Manual (Damped GN) Satu Langkah
LM
Kita gunakan \(D = \operatorname{diag}(J^\top
J)\).
Langkah \(\Delta\) diperoleh dari:
\[(J^\top J + \lambda D)\,\Delta = J^\top r\]
# --- Definisi model eksponensial: Y = beta0 * exp(beta1 * X) + e
f_exp <- function(x, theta) theta[1] * exp(theta[2] * x)
J_exp <- function(x, theta){
e <- exp(theta[2] * x)
cbind(d_beta0 = e, d_beta1 = theta[1] * x * e)
}
resid_vec <- function(x, y, theta) y - f_exp(x, theta)
# --- Satu langkah LM (trial) untuk lambda tertentu
lm_trial <- function(x, y, theta, lambda){
r_old <- resid_vec(x, y, theta)
J <- J_exp(x, theta)
lhs0 <- crossprod(J, J) # J^T J
D <- diag(diag(lhs0)) # diag(J^T J)
rhs <- crossprod(J, r_old) # J^T r
Delta <- solve(lhs0 + lambda * D, rhs)
theta_new <- theta + as.numeric(Delta)
r_new <- resid_vec(x, y, theta_new)
list(Delta = as.numeric(Delta),
theta_new = theta_new,
SSE_old = sum(r_old^2),
SSE_new = sum(r_new^2))
}
# --- Levenberg–Marquardt: Estimasi manual (damped GN)
lm_fit <- function(x, y, theta_init, lambda_init=1e-2, nu=10, max_it=50, tol=1e-8){
th <- theta_init
lam <- lambda_init
hist <- tibble::tibble(iter=0, beta0=th[1], beta1=th[2],
lambda=lam, SSE = sum(resid_vec(x,y,th)^2),
step_norm = NA_real_, accepted = TRUE)
for(k in 1:max_it){
tried <- lm_trial(x, y, th, lam)
if(tried$SSE_new < tried$SSE_old){ # accept
th_new <- tried$theta_new
step_norm <- sqrt(sum(tried$Delta^2))
lam <- lam/nu
hist <- tibble::add_row(hist, iter=k, beta0=th_new[1], beta1=th_new[2],
lambda=lam, SSE=tried$SSE_new,
step_norm=step_norm, accepted = TRUE)
th <- th_new
if(step_norm < tol) break
} else { # reject and increase lambda
lam <- lam*nu
hist <- tibble::add_row(hist, iter=k, beta0=th[1], beta1=th[2],
lambda=lam, SSE=tried$SSE_old,
step_norm=NA_real_, accepted = FALSE)
# do not update theta; try again with bigger lambda
}
}
list(theta = th, history = hist)
}
lm_res <- lm_fit(x, y, theta_init = c(1.0, 0.0), lambda_init = 1e-2, nu=10)
lm_res$history## # A tibble: 11 × 7
## iter beta0 beta1 lambda SSE step_norm accepted
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <lgl>
## 1 0 1 0 0.01 137. NA TRUE
## 2 1 1 0 0.1 137. NA FALSE
## 3 2 1 0 1 137. NA FALSE
## 4 3 1 0 10 137. NA FALSE
## 5 4 1.33 0.134 1 85.7 3.58e-1 TRUE
## 6 5 2.10 0.328 0.1 3.11 7.88e-1 TRUE
## 7 6 1.99 0.316 0.01 0.484 1.10e-1 TRUE
## 8 7 1.98 0.315 0.001 0.465 1.09e-2 TRUE
## 9 8 1.98 0.315 0.0001 0.465 2.59e-4 TRUE
## 10 9 1.98 0.315 0.00001 0.465 3.85e-6 TRUE
## 11 10 1.98 0.315 0.000001 0.465 8.79e-9 TRUE## [1] 1.9753264 0.3148085Kurva Hasil LM
dat$yhat_lm <- f_exp(dat$x, lm_res$theta)
ggplot(dat, aes(x, y)) +
geom_point(size=2) +
geom_line(aes(y = yhat_lm), color="firebrick", linewidth=1) +
labs(title = "Hasil Estimasi Levenberg–Marquardt",
subtitle = paste0("beta0 = ", round(lm_res$theta[1],3),
", beta1 = ", round(lm_res$theta[2],3)),
x = "X", y = "Y") +
theme_minimal()
Pembanding: nls() (Built-in R)
##
## 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.828e-07## b0 b1
## 1.9753263 0.3148085Inferensi (SE & Interval Kepercayaan)
Pada solusî \(\theta\), aproksimasi kovarians diberikan oleh: \[\hat{\sigma}^2 = \dfrac{S(\hat{\theta})}{n - p}\] dan \[\operatorname{Var}(\hat{\theta}) \approx \hat{\sigma}^2 \left(J(\hat{\theta})^\top J(\hat{\theta})\right)^{-1}.\] dengan:
- \(S(\hat{\theta}) = \sum_{i=1}^n r_i(\hat{\theta})^2\) adalah jumlah kuadrat residu pada solusi,
- \(n\) jumlah observasi,
- \(p\) jumlah parameter,
- \(J(\hat{\theta})\) matriks Jacobian dievaluasi padâ \(\theta\).
theta_hat <- lm_res$theta # pakai hasil LM (boleh ganti GN/nls)
r_hat <- resid_vec(x, y, theta_hat)
SSE <- sum(r_hat^2)
n <- length(x); p <- 2
sigma2_hat <- SSE/(n-p)
J_hat <- J_exp(x, theta_hat)
cov_theta <- sigma2_hat * solve(crossprod(J_hat, J_hat))
se_theta <- sqrt(diag(cov_theta))
theta_hat; se_theta## [1] 1.9753264 0.3148085## d_beta0 d_beta1
## 0.16309142 0.01969227Interval kepercayaan Wald 95% untuk tiap parameter \(\theta_j\) adalah:
\[\hat{\theta}_j \pm 1.96\, SE(\hat{\theta}_j).\]
dengan:
- \(\hat{\theta}_j\): estimasi parameter ke-\(j\),
- \(SE(\hat{\theta}_j)\): standar error dari estimasi parameter ke-\(j\).
SOAL
Sensitivity to Start Ulangi GN & LM dengan tebakan awal berbeda-beda (misalnya (0.5, −0.5), (5, 0.5)). Catat konvergensi, jumlah iterasi, dan SSE akhir.
Bandingkan dengan Transformasi Log Generasikan data dengan error multiplikatif: \(Y = \beta_0 e^{\beta_1 X} \cdot \varepsilon,\ \ln \varepsilon \sim N(0,\ \tau^2)\) Estimasi OLS pada \(\ln Y\) dan bandingkan hasilnya dengan NLS.
Diagnostik Residu Untuk hasil terbaik (GN/LM), plot residu vs fitted dan QQ-plot residu. Apakah asumsi normalitas dan homoskedastisitas terlihat wajar?
Tambahan Terapkan GN/LM untuk model Michaelis–Menten dan logistik. Laporkan \(\hat{\theta}\), SE, CI, dan interpretasi substantif parameternya.
jawaban
- Sensitivity to Start (Pengaruh Tebakan Awal) Kita coba jalankan Gauss–Newton (GN) dan Levenberg–Marquardt (LM) dengan beberapa tebakan awal berbeda.
# Data simulasi tetap sama
set.seed(123)
beta0_true <- 2.0
beta1_true <- 0.3
x <- c(0,1,2,3,4,5)
y <- beta0_true * exp(beta1_true * x) + rnorm(length(x), sd = 0.4)
# - Fungsi bantu model eksponensial -
f_exp <- function(x, theta) {
theta[1] * exp(theta[2] * x)
}
resid_vec <- function(x, y, theta) {
y- f_exp(x, theta)
}
J_exp <- function(x, theta) {
b0 <- theta[1]; b1 <- theta[2]
e <- exp(b1 * x)
cbind(
d_beta0 = e,
d_beta1 = b0 * x * e
)
# theta = c(beta0, beta1)
}
# - Gauss-Newton yang stabil (QR + fallback ridge + backtracking) -
gn_fit_safe <- function(x, y, theta_init, max_it = 50, tol = 1e-8,
ridge = 1e-10, max_bt = 8) {
th <- theta_init
sse_prev <- sum(resid_vec(x, y, th)^2)
for (k in 1:max_it) {
r <- resid_vec(x, y, th)
J <- J_exp(x, th)
# Coba solusi LS dengan QR (lebih stabil)
Delta <- tryCatch(
qr.coef(qr(J), r),
error = function(e) rep(NA_real_, ncol(J))
)
# Jika gagal/NA, pakai ridge kecil pada normal equations
if (any(!is.finite(Delta))) {
JTJ <- crossprod(J, J)
rhs <- crossprod(J, r)
Delta <- tryCatch(
solve(JTJ + ridge * diag(ncol(J)), rhs),
error = function(e) rep(0, ncol(J))
)
}
# Backtracking line search: kurangi langkah jika SSE naik
step <- 1.0
improved <- FALSE
for (bt in 0:max_bt) {
th_new <- th + as.numeric(step * Delta)
sse_new <- sum(resid_vec(x, y, th_new)^2)
if (is.finite(sse_new) && (sse_new <= sse_prev)) {
th <- th_new
sse_prev <- sse_new
improved <- TRUE
break
}
step <- step / 2
}
# Hentikan jika tidak ada perbaikan atau langkah sangat kecil
if (!improved || sqrt(sum((step * Delta)^2)) < tol) break
}
th
}
# - Coba beberapa starting values (hindari beta0=0) -
starts <- list(c(0.5,-0.5), c(5, 0.5), c(1, 0.0))
results <- lapply(starts, function(st) gn_fit_safe(x, y, st))
names(results) <- c("start (0.5,-0.5)", "start (5,0.5)", "start (1,0)")
results## $`start (0.5,-0.5)`
## [1] 1.9753264 0.3148085
##
## $`start (5,0.5)`
## [1] 1.9753264 0.3148085
##
## $`start (1,0)`
## [1] 1.9753264 0.3148085Kesimpulan
Hasil akhirnya relatif konsisten (mendekati \(\beta_0 \approx 2, \beta_1 \approx 0.3\)) meskipun starting value berbeda. Namun, jumlah iterasi dan kecepatan konvergensi bisa berbeda.
- Bandingkan dengan Transformasi Log
Untuk data dengan error multiplikatif, model bisa dilog-kan agar menjadi linear.
set.seed(456)
y_mult <- beta0_true * exp(beta1_true * x) * exp(rnorm(length(x), sd = 0.2))
df <- data.frame(x, y_mult)
# Estimasi OLS pada log(y)
fit_lm <- lm(log(y_mult) ~ x, data = df)
summary(fit_lm)##
## Call:
## lm(formula = log(y_mult) ~ x, data = df)
##
## Residuals:
## 1 2 3 4 5 6
## -0.20616 0.19319 0.23530 -0.19636 -0.05516 0.02919
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.63060 0.15252 4.135 0.01444 *
## x 0.29371 0.05038 5.830 0.00431 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2107 on 4 degrees of freedom
## Multiple R-squared: 0.8947, Adjusted R-squared: 0.8684
## F-statistic: 33.99 on 1 and 4 DF, p-value: 0.004312# konversi kembali
beta0_hat <- exp(coef(fit_lm)[1])
beta1_hat <- coef(fit_lm)[2]
c(beta0_hat, beta1_hat)## (Intercept) x
## 1.8787395 0.2937094Kesimpulan:
Dengan error multiplikatif, model linear log cocok → estimasi lebih sederhana (OLS), hasilnya dekat dengan nilai parameter sebenarnya.
- Diagnostik Residu
Setelah fitting, cek residu.
fit_nls <- nls(
y ~ b0 * exp(b1 * x),
start = list(b0 = 1, b1 = 0),
trace = FALSE
)
resid_fit <- resid(fit_nls)
fitted_fit <- fitted(fit_nls)
par(mfrow = c(1, 2))
plot(
fitted_fit, resid_fit,
main = "Residu vs Fitted",
xlab = "Fitted", ylab = "Residuals"
)
abline(h = 0, lty = 2, col = "red")
qqnorm(resid_fit)
qqline(resid_fit, col = "blue")
Interpretasi:
Plot residu vs fitted → apakah ada pola sistematis? Jika tidak, asumsi homoskedastisitas cukup terpenuhi. QQ-plot → memeriksa apakah distribusi residu mendekati normal.
- Tambahan: Model Michaelis–Menten
Gunakan nls() dengan data simulasi Michaelis–Menten.
set.seed(789)
Vmax_true <- 2
Km_true <- 1.2
X <- seq(0, 6, length.out = 40)
Y <- (Vmax_true * X) / (Km_true + X) + rnorm(length(X), sd = 0.05)
df_mm <- data.frame(X, Y)
fit_mm <- nls(
Y ~ (Vmax * X) / (Km + X),
data = df_mm,
start = list(Vmax = 1.5, Km = 1)
)
summary(fit_mm)##
## Formula: Y ~ (Vmax * X)/(Km + X)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Vmax 1.9981 0.0234 85.40 <2e-16 ***
## Km 1.2415 0.0491 25.28 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03623 on 38 degrees of freedom
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 1.317e-06## 2.5% 97.5%
## Vmax 1.952336 2.04651
## Km 1.146484 1.34395Interpretasi
- \(\hat{V}_{\max}\) mendekati 2 = kapasitas maksimum.
- \(\hat{K}_m\) mendekati 1.2 = konsentrasi yang menghasilkan setengah \(V_{\max}\).
- CI (interval kepercayaan) memberi gambaran ketidakpastian estimasi.
Ringkasan Jawaban
- GN & LM sensitif terhadap tebakan awal, tetapi dengan tebakan masuk akal tetap konvergen.
- Untuk error multiplikatif, model bisa dilog-kan menjadi linear → estimasi mudah dengan OLS.
- Diagnostik residu penting untuk validasi asumsi normalitas & homoskedastisitas.
- Pada Michaelis–Menten, interpretasi parameter \(V_{\max}\) dan \(K_m\) sangat substantif (kapasitas & efisiensi sistem).
2.4.2 Inferensi: uji signifikansi parameter, interval kepercayaan,goodness-of-fit (AIC, R²).
Setelah parameter model nonlinear̂ 𝜃 diperoleh (misal dengan nonlinear least squares), kita ingin melakukan inferensi:
- Apakah parameter signifikan (= 0)?
- Bagaimana ketidakpastian estimasi (interval kepercayaan)?
- Seberapa baik model cocok dengan data (goodness-of-fit)?
2.4.2.1 Uji Signifikansi Parameter
Hipotesis umum untuk tiap parameter \(\theta_j\):
\(H_0 \colon \theta_j = 0\) vs \(H_1 \colon \theta_j \ne 0\).
Statistik uji Wald: \[t_j = \frac{\hat{\theta}_j}{SE(\hat{\theta}_j)},\] dengan \(SE(\hat{\theta}_j)\) diperoleh dari akar diagonal matriks kovarianŝ \(\operatorname{Var}(\hat{\theta})\).
Jika \(|t_j|\) cukup besar, atau \(p\)-value < 0.05, maka \(\theta_j\) signifikan.
2.4.2.2 Interval Kepercayaan
Wald confidence interval 95% untuk tiap parameter: \[\hat{\theta}_j \pm 1.96 \cdot SE(\hat{\theta}_j).\] Lebih umum, untuk taraf kepercayaan \((1 - \alpha)\): \[\hat{\theta}_j \pm z_{1-\alpha/2} \cdot SE(\hat{\theta}_j).\] Dengan \(z_{1-\alpha/2}\) = kuantil distribusi normal standar.
2.4.2.3 Goodness-of-Fit
𝑅² (Pseudo-𝑅² pada regresi nonlinear)
\[R^2 = 1 -
\dfrac{\text{SSE}}{\text{SST}},\] dengan:
- SSE = \(\sum (y_i - \hat{y}_i)^2\)
(residual sum of squares),
- SST = \(\sum (y_i - \bar{y})^2\)
(total sum of squares).
Interpretasi: proporsi variasi 𝑦 yang dijelaskan model. Nilai mendekati 1 → fit baik.
Akaike Information Criterion (AIC)
\[AIC = n \cdot
\ln\!\left(\dfrac{\text{SSE}}{n}\right) + 2p,\] dengan:
- \(n\) = jumlah observasi,
- \(p\) = jumlah parameter.
Interpretasi: lebih kecil → model lebih baik (dengan penalti kompleksitas).
2.4.2.4 Studi Kasus: Model Eksponensial
Kita gunakan model: \[y = \beta_0 e^{\beta_1 x} + \varepsilon,\ \varepsilon \sim N(0,\ \sigma^2).\]
set.seed(123)
library(broom)
library(ggplot2)
# Simulasi data
beta0_true <- 2; beta1_true <- 0.3
x <- 0:10
y <- beta0_true * exp(beta1_true * x) + rnorm(length(x), sd = 0.5)
dat <- data.frame(x, y)
# Estimasi dengan nls()
fit <- nls(
y ~ b0 * exp(b1 * x),
data = dat,
start = list(b0 = 1, b1 = 0.1)
)
summary(fit)##
## Formula: y ~ b0 * exp(b1 * x)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b0 2.005728 0.096821 20.72 6.66e-09 ***
## b1 0.300091 0.005399 55.59 9.93e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.515 on 9 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 4.467e-062.4.2.5 Uji Signifikansi Parameter
## # A tibble: 2 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 b0 2.01 0.0968 20.7 6.66e- 9 1.79 2.23
## 2 b1 0.300 0.00540 55.6 9.93e-13 0.288 0.313Interpretasi: jika p-value < 0.05 → parameter signifikan berbeda dari 0.
2.4.3 Goodness-of-Fit
n <- nrow(dat)
p <- length(coef(fit))
yhat <- fitted(fit)
SSE <- sum(resid(fit)^2)
SST <- sum((y - mean(y))^2)
R2 <- 1 - SSE / SST
AIC <- n * log(SSE / n) + 2 * p
list(SSE = SSE, R2 = R2, AIC = AIC)## $SSE
## [1] 2.387234
##
## $R2
## [1] 0.9984577
##
## $AIC
## [1] -12.805362.4.3.1 Studi Kasus: Model Eksponensial
Kita gunakan model: \[y = \beta_0 e^{\beta_1 x} + \varepsilon,\ \varepsilon \sim N(0,\ \sigma^2).\]
set.seed(123)
library(broom)
library(ggplot2)
# Simulasi data
beta0_true <- 2; beta1_true <- 0.3
x <- 0:10
y <- beta0_true * exp(beta1_true * x) + rnorm(length(x), sd = 0.5)
dat <- data.frame(x, y)
# Estimasi dengan nls()
fit <- nls(
y ~ b0 * exp(b1 * x),
data = dat,
start = list(b0 = 1, b1 = 0.1)
)
summary(fit)##
## Formula: y ~ b0 * exp(b1 * x)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b0 2.005728 0.096821 20.72 6.66e-09 ***
## b1 0.300091 0.005399 55.59 9.93e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.515 on 9 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 4.467e-06- Uji Signifikansi Parameter
## # A tibble: 2 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 b0 2.01 0.0968 20.7 6.66e- 9 1.79 2.23
## 2 b1 0.300 0.00540 55.6 9.93e-13 0.288 0.313Interpretasi: jika p-value < 0.05 → parameter signifikan berbeda dari 0.
- Interval Kepercayaan 95%
## 2.5% 97.5%
## b0 1.7917267 2.2341574
## b1 0.2879737 0.3126612Interval kepercayaan memberi rentang nilai parameter yang konsisten dengan data.
- Goodness-of-Fit
n <- nrow(dat)
p <- length(coef(fit))
yhat <- fitted(fit)
SSE <- sum(resid(fit)^2)
SST <- sum((y - mean(y))^2)
R2 <- 1 - SSE / SST
AIC <- n * log(SSE / n) + 2 * p
list(SSE = SSE, R2 = R2, AIC = AIC)## $SSE
## [1] 2.387234
##
## $R2
## [1] 0.9984577
##
## $AIC
## [1] -12.805362.5 Tugas
2.5.1 Studi Data Nyata: Ambil dataset pertumbuhan tanaman/hewan (atau data ekonomi) yang menunjukkan pola kurva nonlinear.
• Fit-kan model nonlinear dengan nls().
• Laporkan estimasi parameter, interpretasi biologis/ekonomisnya, serta goodness-of-fit.
jawab
2.5.1.1 Studi Data Nyata
Pada bagian ini, kita akan menggunakan dataset pertumbuhan tanaman (misalnya data pertumbuhan jagung). Dataset diambil dari paket datasets atau dibuat simulasi sederhana yang menunjukkan pola kurva nonlinear (logistik).
Data
set.seed(123)
time <- 1:10
# Simulasi pertumbuhan tanaman (logistik)
growth <- 100 / (1 + exp(-0.5*(time-5))) + rnorm(10, 0, 2)
data <- data.frame(time, growth)
head(data)## time growth
## 1 1 10.79934
## 2 2 17.78220
## 3 3 30.01156
## 4 4 37.89508
## 5 5 50.25858
## 6 6 65.676062.5.1.2 Model Nonlinear (nls)
Kita fit model nonlinear logistik dengan fungsi
nls():
\[ y = \frac{A}{1 + e^{-k(t - t_0)}} \]
dengan parameter: - \(A\): Asimtot maksimum (kapasitas pertumbuhan maksimum) - \(k\): Laju pertumbuhan - \(t_0\): Titik tengah pertumbuhan
nls_model <- nls(growth ~ A / (1 + exp(-k*(time - t0))),
data = data,
start = list(A = 100, k = 0.5, t0 = 5))
summary(nls_model)##
## Formula: growth ~ A/(1 + exp(-k * (time - t0)))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 95.96414 2.30668 41.60 1.21e-09 ***
## k 0.52700 0.02892 18.22 3.71e-07 ***
## t0 4.72209 0.13782 34.26 4.68e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.735 on 7 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 6.867e-06Interpretasi parameter: - \(A\) ≈ kapasitas maksimum (bobot maksimum tanaman). - \(k\) ≈ laju pertumbuhan. - \(t_0\) ≈ waktu ketika pertumbuhan mencapai setengah maksimum.
2.5.1.3 Goodness-of-fit
# Prediksi
data$pred_nls <- predict(nls_model)
# R-squared manual
SSE <- sum((data$growth - data$pred_nls)^2)
SST <- sum((data$growth - mean(data$growth))^2)
R2_nls <- 1 - SSE/SST
R2_nls## [1] 0.9972567Visualisasi:
plot(data$time, data$growth, pch=19, main="Nonlinear Fit vs Data")
lines(data$time, data$pred_nls, col="blue", lwd=2)
2.5.2 Bandingkan dengan Regresi Linear
• Fit model linear sederhana pada data yang sama. • Bandingkan SSE, R², dan visualisasi dengan model nonlinear. • Diskusikan mengapa model nonlinear lebih sesuai.
jawab
2.5.2.1 Fit Model Linear
##
## Call:
## lm(formula = growth ~ time, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.9280 -1.9006 -0.4940 0.8489 6.5009
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.7616 2.6760 0.658 0.529
## time 9.5689 0.4313 22.187 1.8e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.917 on 8 degrees of freedom
## Multiple R-squared: 0.984, Adjusted R-squared: 0.982
## F-statistic: 492.3 on 1 and 8 DF, p-value: 1.8e-082.5.2.2 Bandingkan SSE dan R²
# Prediksi linear
data$pred_lm <- predict(lm_model)
SSE_lm <- sum((data$growth - data$pred_lm)^2)
R2_lm <- summary(lm_model)$r.squared
c(SSE_nls = SSE, R2_nls = R2_nls, SSE_lm = SSE_lm, R2_lm = R2_lm)## SSE_nls R2_nls SSE_lm R2_lm
## 21.0598976 0.9972567 122.7610884 0.98400892.5.2.3 Diskusi
- Model nonlinear lebih baik karena sesuai dengan pola biologis pertumbuhan (logistik).
- Regresi linear hanya cocok untuk fase awal pertumbuhan, tapi gagal menangkap saturasi.
- Nilai R² lebih tinggi pada model nonlinear, menandakan kecocokan model lebih baik.
Visualisasi perbandingan:
plot(data$time, data$growth, pch=19, main="Linear vs Nonlinear Fit")
lines(data$time, data$pred_nls, col="blue", lwd=2)
lines(data$time, data$pred_lm, col="red", lwd=2, lty=2)
legend("bottomright", legend=c("Data","Nonlinear","Linear"),
col=c("black","blue","red"), pch=c(19,NA,NA), lty=c(NA,1,2))
2.5.3 Simulasi: Pertumbuhan Eksponensial (Bakteri)
• Gunakan nls() untuk estimasi parameter. • Laporkan hasil dan interpretasi.
jawab
Kita buat simulasi data pertumbuhan bakteri mengikuti model eksponensial:
\[ y = y_0 \, e^{rt} \]
dengan parameter: - \(y_0\): jumlah awal bakteri - \(r\): laju pertumbuhan - \(t\): waktu
2.5.3.1 Simulasi Data
2.5.3.2 Fit Model Nonlinear
nls_bac <- nls(growth_bac ~ y0 * exp(r * time_bac),
data = data_bac,
start = list(y0=5, r=0.2))
summary(nls_bac)##
## Formula: growth_bac ~ y0 * exp(r * time_bac)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## y0 10.45162 0.50386 20.74 3.06e-08 ***
## r 0.39470 0.00521 75.75 1.03e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.652 on 8 degrees of freedom
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 3.322e-062.5.3.3 Visualisasi
data_bac$pred <- predict(nls_bac)
plot(data_bac$time_bac, data_bac$growth_bac, pch=19,
main="Eksponensial Fit Pertumbuhan Bakteri")
lines(data_bac$time_bac, data_bac$pred, col="blue", lwd=2)
Interpretasi:
- Estimasi \(y_0\) ≈ jumlah awal
bakteri. - Estimasi \(r\) ≈ laju
pertumbuhan.
Model eksponensial cocok untuk pertumbuhan awal populasi bakteri.
2.5.3.4 Interpretasi
- Model nonlinear (logistik/eksponensial) memberikan hasil lebih baik
dibanding model linear.
- Parameter dapat diinterpretasikan secara biologis/ekonomis.
- Simulasi eksponensial menunjukkan bahwa
nls()dapat mengestimasi parameter dengan baik.
2.5.4 Kesimpulan
- Regresi nonlinear diperlukan ketika hubungan X–Y tidak bisa dijelaskan dengan garis lurus.
- Parameter nonlinear sering punya interpretasi substantif (misalnya: laju pertumbuhan, kapasitas maksimum).
- Estimasi dilakukan dengan algoritma iteratif (nls() di R).
- Visualisasi sangat penting untuk membandingkan data vs kurva model.
2.6 Appendix: A. Gauss–Newton (GN) untuk Model Eksponensial
2.6.1 Model dan Asumsi
Model regresi eksponensial dengan error aditif: \[Y_i = f_i(\theta) + \varepsilon_i = \beta_0 e^{\beta_1 x_i} + \varepsilon_i,\ \varepsilon_i\ \text{iid} \sim \mathcal{N}(0,\ \sigma^2),\ i = 1, \ldots, n.\]
Vektor parameter:
\[\theta = (\beta_0, \beta_1)^\top,\ p = 2.\]
Catatan: Di bawah asumsi normal, NLS = MLE untuk \(\theta\) dan \(\sigma^2\).
2.6.2 Residual, Fungsi Tujuan, dan Jacobian
Residual dan fungsi objektif (jumlah kuadrat residual): \[r_i(\theta) = y_i - \beta_0 e^{\beta_1 x_i}, \quad S(\theta) = \sum_{i=1}^{n} [r_i(\theta)]^2 = \lVert r(\theta) \rVert_2^2.\]
Jacobian \(J(\theta) \in \mathbb{R}^{n \times p}\) dengan baris ke-\(i\):
\[J_{i1} = \dfrac{\partial f_i}{\partial \beta_0} = e^{\beta_1 x_i}, \quad J_{i2} = \dfrac{\partial f_i}{\partial \beta_1} = \beta_0 x_i e^{\beta_1 x_i}.\]
2.6.3 Gradien dan (Aproksimasi) Hessian dari 𝑆(𝜃)
Gradien di \(\theta\): \[\nabla S(\theta) = g(\theta) = -2\, J(\theta)^\top r(\theta).\]
Aproksimasi Hessian Gauss–Newton: \[H_{\text{GN}}(\theta) \approx 2\, J(\theta)^\top J(\theta).\]
(Hessian eksak: \(H = 2 J^\top J - 2 \sum_{i=1}^{n} r_i H_i\), dengan \(H_i\) Hessian \(f_i\).)
2.6.4 Satu Langkah Iterasi GN (Linearized Least Squares)
Pada iterasi ke-\(k\), langkah perbaikan \(\Delta^{(k)}\) diperoleh dari: \[(J^{(k)\top} J^{(k)})\,\Delta^{(k)} = J^{(k)\top} r^{(k)},\] dengan \[r^{(k)} = \bigl(y_i - \beta_0^{(k)} e^{\beta_1^{(k)} x_i}\bigr)_{i=1}^n,\quad J^{(k)} = J(\theta^{(k)}).\]
Pembaruan parameter: \[\theta^{(k+1)} = \theta^{(k)} + \Delta^{(k)}.\]
Kriteria henti (contoh): \[\lVert \Delta^{(k)} \rVert_\infty < \varepsilon_\theta \quad \text{atau} \quad \bigl|S^{(k+1)} - S^{(k)}\bigr| < \varepsilon_S \quad \text{atau} \quad \lVert g^{(k)} \rVert_\infty < \varepsilon_g.\]
Alternatif redaman (Levenberg–Marquardt): \[(J^{(k)\top} J^{(k)} + \lambda^{(k)} D^{(k)})\,\Delta^{(k)} = J^{(k)\top} r^{(k)},\quad D^{(k)} = I \ \text{atau}\ \operatorname{diag}(J^{(k)\top} J^{(k)}).\]
2.6.5 Estimasi Varians Error dan Kovarians Parameter
Setelah konvergen dî \(\theta\), taksiran varians error: \[\hat{\sigma}^2 = \dfrac{S(\hat{\theta})}{n - p}.\]
Kovarians asimtotik parameter (Wald): \[\widehat{\mathrm{Cov}}(\hat{\theta}) \approx \hat{\sigma}^2 \left(J(\hat{\theta})^\top J(\hat{\theta})\right)^{-1}.\]
Galat baku: \[SE(\hat{\theta}_j) = \sqrt{\hat{\sigma}^2 \left[\left(J^\top J\right)^{-1}\right]_{jj}}\ \big|_{\theta=\hat{\theta}}.\]
(Opsional, robust sandwich—untuk heteroskedastisitas): \[\widehat{\text{Covrobust}}(\hat{\theta}) = \left(J^\top J\right)^{-1} J^\top \operatorname{diag}\!\big(r_i(\hat{\theta})^2\big)\, J \left(J^\top J\right)^{-1}.\]
2.6.6 Inferensi: Statistik-𝑧, CI Wald, dan Uji Wald Bersama
Asimtotik normalitas (Wald): \[\hat{\theta} \ \dot{\sim}\ \mathcal{N}\!\big(\theta,\ \sigma^2 (J^\top J)^{-1}\big).\]
Statistik-\(z\) satu parameter untuk hipotesis \(H_0 \colon \theta_j = \theta_{j,0}\): \[z_j = \frac{\hat{\theta}_j - \theta_{j,0}}{SE(\hat{\theta}_j)} \ \dot{\sim}\ \mathcal{N}(0,1).\]
CI Wald \(100(1-\alpha)\%\) untuk \(\theta_j\): \[\hat{\theta}_j \pm z_{1-\alpha/2}\, SE(\hat{\theta}_j),\] dengan \(z_{1-\alpha/2}\) kuantil standar normal (mis. \(z_{0.975} \approx 1.96\)).
(Pada sampel kecil, dapat memakai \(t_{n-p,\ 1-\alpha/2}\).)
Uji Wald bersama (\(q\) kendala linear \(R\theta = r\)):
\[W = (R\hat{\theta} - r)^\top \big[\, R\, \widehat{\mathrm{Cov}}(\hat{\theta})\, R^\top \big]^{-1} (R\hat{\theta} - r)\ \dot{\sim}\ \chi^2_q.\]
2.6.7 Interval untuk Rata-rata Kondisional dan Prediksi
Untuk titik baru \(x^\star\), prediksi mean: \[\hat{y}^\star = f(x^\star, \hat{\theta}) = \hat{\beta}_0 \, e^{\hat{\beta}_1 x^\star}.\]
Linearization delta method: turunan terhadap parameter dî \(\theta\), \[g^\star = \left[\frac{\partial f(x^\star,\theta)}{\partial \beta_0}\ \ \frac{\partial f(x^\star,\theta)}{\partial \beta_1}\right]_{\theta=\hat{\theta}} = \left[ e^{\hat{\beta}_1 x^\star}\ \ \hat{\beta}_0 x^\star e^{\hat{\beta}_1 x^\star} \right].\]
Varian prediksi mean: \[\operatorname{Var}(\hat{y}^\star) \approx {g^\star}^\top\, \widehat{\operatorname{Cov}}(\hat{\theta})\, g^\star.\]
CI untuk mean \(f(x^\star,\theta)\): \[\hat{y}^\star \pm z_{1-\alpha/2}\, \sqrt{\operatorname{Var}(\hat{y}^\star)}.\]
PI (prediction interval) untuk \(Y^\star\): \[\hat{y}^\star \pm z_{1-\alpha/2}\, \sqrt{\operatorname{Var}(\hat{y}^\star) + \hat{\sigma}^2}.\]
2.6.8 Ringkas Algoritma GN (Praktis)
Pilih tebakan awal \(\theta^{(0)}\) (mis. via regresi \(\ln y\) pada \(x\) untuk inisialisasi).
Ulangi hingga konvergen:
- Hitung \(r^{(k)}\) dan \(J^{(k)}\).
- Pecahkan \((J^{(k)\top} J^{(k)})\,\Delta^{(k)} = J^{(k)\top} r^{(k)}\).
- Set \(\theta^{(k+1)} = \theta^{(k)} + \Delta^{(k)}\).
- Hentikan jika kriteria henti terpenuhi; jika \(S\) tidak turun, gunakan line-search/LM.
- Padâ \(\theta\), hitung \(\hat{\sigma}^2\), \(\widehat{\text{Cov}}(\hat{\theta})\), SE, \(z\)-stat, dan CI Wald.
Catatan
- Sensitivitas starting value: GN/LM dapat sensitif; gunakan beberapa tebakan dan bandingkan SSE akhir.
- Nonlinearitas kuat: CI Wald bisa kurang akurat; pertimbangkan profil likelihood.
- Skala: Standarisasi \(x\) dapat memperbaiki kondisioning \(J^\top J\).
2.6.9 Perhitungan Manual
Contoh GN Manual (n = 5) untuk Model Eksponensial
Data
Kita gunakan 5 observasi:
\[x = (0, 1, 2, 3, 4)^\top,\]
\[y = (2.100000, 2.649718, 3.724238,
4.819206, 6.790234)^\top.\]
Model
\[y_i = \beta_0 e^{\beta_1 x_i} +
\varepsilon_i,\ \ \varepsilon_i \sim \mathcal N(0,\
\sigma^2).\]
Iterasi 0 (tebakan awal dari regresi \(\ln y\) pada \(x\))
\[
\theta^{(0)} =
\begin{bmatrix}
\beta^{(0)}_0\\[2pt]
\beta^{(0)}_1
\end{bmatrix}
=
\begin{bmatrix}
2.043817\\[2pt]
0.294525
\end{bmatrix}.
\]
Nilai fungsi, residual, dan Jacobian di \(\theta^{(0)}\):
\[ f^{(0)} = \beta^{(0)}_0 e^{\beta^{(0)}_1 x} = \begin{bmatrix} 2.043817\\ 2.743801\\ 3.683523\\ 4.945088\\ 6.638725 \end{bmatrix},\qquad r^{(0)} = y - f^{(0)} = \begin{bmatrix} \ \ 0.056183\\ -0.094083\\ \ \ 0.040715\\ -0.125881\\ \ \ 0.151509 \end{bmatrix}. \]
Jacobian (baris ke-\(i\): \([\partial f_i/\partial \beta_0,\ \partial f_i/\partial \beta_1]\)): \[ J^{(0)} = \begin{bmatrix} 1.000000 & \ \ 0.000000\\ 1.342489 & \ \ 2.743801\\ 1.802276 & \ \ 7.367045\\ 2.419536 & 14.835263\\ 3.248200 & 26.554899 \end{bmatrix}. \]
Matriks normal dan vektor kanan: \[ J^{(0)\top}J^{(0)} = \begin{bmatrix} \ 22.455430 & 139.111033\\ 139.111033 & 987.049457 \end{bmatrix},\qquad J^{(0)\top} r^{(0)} = \begin{bmatrix} 0.190815\\ 2.197633 \end{bmatrix}. \]
Langkah GN (persamaan normal): \[ \big(J^{(0)\top}J^{(0)}\big)\,\Delta^{(0)} = J^{(0)\top} r^{(0)} \ \ \Rightarrow\ \ \Delta^{(0)} = \begin{bmatrix} -0.041728\\ \ \ 0.008108 \end{bmatrix}. \]
Pembaruan parameter dan SSE: \[ \theta^{(1)} = \theta^{(0)} + \Delta^{(0)} = \begin{bmatrix} 2.002088\\ 0.302633 \end{bmatrix},\qquad S\big(\theta^{(0)}\big)=0.052467,\ \ S\big(\theta^{(1)}\big)=0.042495. \]
Iterasi 1
Hitung ulang \(J^{(1)},\ r^{(1)}\) pada
\(\theta^{(1)}\): \[
J^{(1)} =
\begin{bmatrix}
1.000000 & \ \ 0.000000\\
1.353417 & \ \ 2.709661\\
1.831738 & \ \ 7.334604\\
2.479106 & 14.890170\\
3.355265 & 26.870152
\end{bmatrix},\qquad
r^{(1)} =
\begin{bmatrix}
\ \ 0.097912\\
-0.059943\\
\ \ 0.056936\\
-0.144184\\
\ \ 0.072696
\end{bmatrix}.
\]
\[ J^{(1)\top}J^{(1)} = \begin{bmatrix} \ 23.590779 & 144.173187\\ 144.173187 & 1004.860917 \end{bmatrix},\qquad J^{(1)\top} r^{(1)} = \begin{bmatrix} 0.007541\\ 0.061600 \end{bmatrix}. \]
Langkah GN berikutnya: \[ \Delta^{(1)} = \begin{bmatrix} -0.000446\\ \ \ 0.000125 \end{bmatrix},\qquad \theta^{(2)} = \begin{bmatrix} 2.001642\\ 0.302758 \end{bmatrix},\qquad S\big(\theta^{(2)}\big)=0.042491. \]
Terlihat 𝑆(𝜃) menurun dan perubahan parameter makin kecil → konvergen.
2.7 Appendix: Levenberg–Marquardt (LM) untuk Model Eksponensial
2.7.1 Model dan Asumsi
Kita gunakan model regresi eksponensial dengan error aditif: \[Y_i = f_i(\theta) + \varepsilon_i = \beta_0 e^{\beta_1 x_i} + \varepsilon_i,\ \varepsilon_i\ \text{iid} \sim \mathcal{N}(0,\ \sigma^2),\ i = 1, \ldots, n.\]
Vektor parameter: \[\theta = (\beta_0, \beta_1)^\top,\ p = 2.\]
2.7.2 Ide Inti LM (Damped Gauss–Newton)
LM menggabungkan Gauss–Newton dan Gradient Descent.
Iterasi ke-\(k\) menyelesaikan sistem modifikasi:
\[(J^{(k)\top} J^{(k)} + \lambda^{(k)} D^{(k)})\, \Delta^{(k)} = J^{(k)\top} r^{(k)},\]
dengan:
- \(J^{(k)}\) : Jacobian di \(\theta^{(k)}\)
- \(r^{(k)} = (y_i - \beta^{(k)}_0 e^{\beta^{(k)}_1 x_i})_{i=1}^n\)
- \(D^{(k)} = I\) atau \(\operatorname{diag}(J^{(k)\top} J^{(k)})\)
- \(\lambda^{(k)} > 0\) : faktor redaman
Jika \(\lambda^{(k)} \to 0\), metode mendekati GN; jika \(\lambda^{(k)}\) besar, mendekati Gradient Descent.
2.7.3 Residual, Fungsi Tujuan, dan Jacobian
Residual dan fungsi SSE sama dengan GN: \[r_i(\theta) = y_i - \beta_0 e^{\beta_1 x_i}, \quad S(\theta) = \sum_{i=1}^{n} r_i(\theta)^2.\]
Jacobian (baris ke-\(i\)): \[J_{i1} = e^{\beta_1 x_i}, \quad J_{i2} = \beta_0 x_i e^{\beta_1 x_i}.\]
2.7.4 Update Parameter dengan Redaman
Pembaruan parameter: \[\theta^{(k+1)} = \theta^{(k)} + \Delta^{(k)}.\]
Aturan adaptasi \(\lambda\):
- Jika \(S(\theta^{(k+1})) < S(\theta^{(k)})\), kecilkan \(\lambda\) (lebih agresif).
- Jika tidak, besarkan \(\lambda\) (lebih konservatif) dan coba ulangi.
Sehingga: \[ \lambda^{(k+1)} = \begin{cases} \lambda^{(k)}/\nu, & \text{jika } S(\theta^{(k+1})) < S(\theta^{(k)}),\\ \nu\,\lambda^{(k)}, & \text{jika tidak.} \end{cases} \]
(dengan \(\nu > 1\), misalnya \(\nu = 10\)).
2.7.5 Gradien dan (Aproksimasi) Hessian dalam LM
LM memodifikasi matriks normal: \[H_{\text{LM}}^{(k)} = J^{(k)\top} J^{(k)} + \lambda^{(k)} D^{(k)}.\]
Gradien tetap sama: \[g^{(k)} = -2\, J^{(k)\top} r^{(k)}.\]
2.7.6 Estimasi Varians Error dan Kovarians Parameter
Setelah konvergen dî \(\theta\): \[\hat{\sigma}^2 = \dfrac{S(\hat{\theta})}{n - p}.\]
Kovarians asimtotik (Wald): \[\widehat{\mathrm{Cov}}(\hat{\theta}) \approx \hat{\sigma}^2 (J^\top J)^{-1}\big|_{\theta=\hat{\theta}}.\]
SE parameter: \[SE(\hat{\theta}_j) = \sqrt{\hat{\sigma}^2 \,\big[(J^\top J)^{-1}\big]_{jj}}.\]
2.7.7 Inferensi: Statistik-𝑧 dan CI
Distribusi asimtotik: \[\hat{\theta} \ \dot{\sim}\ \mathcal{N}\!\big(\theta,\ \sigma^2 (J^\top J)^{-1}\big).\]
Statistik-\(z\): \[z_j = \frac{\hat{\theta}_j - \theta_{j,0}}{SE(\hat{\theta}_j)} \ \dot{\sim}\ \mathcal{N}(0,1).\]
CI Wald \(100(1 - \alpha)\%\) untuk \(\theta_j\): \[\hat{\theta}_j \pm z_{1-\alpha/2}\, SE(\hat{\theta}_j).\]
2.7.8 Interval Prediksi dan Mean
Prediksi mean di \(x^\star\): \[\hat{y}^\star = f(x^\star, \hat{\theta}) = \hat{\beta}_0 e^{\hat{\beta}_1 x^\star}.\]
Gradien terhadap parameter: \[ g^\star = \begin{bmatrix} e^{\hat{\beta}_1 x^\star} & \hat{\beta}_0 x^\star e^{\hat{\beta}_1 x^\star} \end{bmatrix}. \]
Varian prediksi mean: \[\operatorname{Var}(\hat{y}^\star) \approx {g^\star}^\top\, \widehat{\operatorname{Cov}}(\hat{\theta})\, g^\star.\]
CI mean: \[\hat{y}^\star \pm z_{1-\alpha/2}\, \sqrt{\operatorname{Var}(\hat{y}^\star)}.\]
PI (prediction interval): \[\hat{y}^\star \pm z_{1-\alpha/2}\, \sqrt{\operatorname{Var}(\hat{y}^\star) + \hat{\sigma}^2}.\]
2.7.9 Ringkas Algoritma LM
- Pilih tebakan awal \(\theta^{(0)}\), pilih \(\lambda^{(0)} > 0\).
- Iterasi:
- Hitung \(r^{(k)}\) dan \(J^{(k)}\).
- Selesaikan \((J^{(k)\top}J^{(k)} + \lambda^{(k)}D^{(k)})\Delta^{(k)} = J^{(k)\top}r^{(k)}\).
- Update \(\theta^{(k+1)} = \theta^{(k)} + \Delta^{(k)}\).
- Adaptasi \(\lambda\) sesuai aturan sukses/gagal.
- Hentikan bila kriteria konvergensi tercapai.
- Setelah konvergen: hitung \(\hat{\sigma}^2\), SE, \(z\)-statistik, dan CI.
Catatan - LM lebih stabil daripada GN ketika \(J^\top J\) hampir singular atau starting value jauh. - \(\lambda\) adaptif menjaga keseimbangan eksplorasi (besar \(\lambda\)) dan efisiensi (kecil \(\lambda\)). - Interval Wald tetap bisa kurang akurat → gunakan profil likelihood bila perlu.
3 Robust Regression
Bab ini menyajikan materi komprehensif tentang outlier dalam analisis regresi, mencakup:
- (i) definisi dan klasifikasi,
- (ii) teknik identifikasi grafis dan statistik,
- (iii) strategi penanganan, serta
- (iv) ide dasar dan praktik robust regression berikut contoh kode R yang siap dijalankan.
3.1 Konsep Dasar Outlier
Outlier adalah observasi yang menyimpang jauh dari pola umum data. Dalam regresi, pembedaan penting:
- Outlier di ruang respons (Y-space): residual sangat besar.
- Outlier di ruang prediktor (X-space): nilai prediktor jauh dari pusat sebaran.
- Leverage point: observasi dengan nilai X ekstrem (potensi memengaruhi slope).
- Influential point: observasi yang jika dihapus, koefisien berubah nyata (diukur mis. Cook’s distance).
Secara umum, regresi OLS meminimalkan jumlah kuadrat residual
\[\hat{\beta} = \arg\min_{\beta} \sum_i (y_i - x_i^\top \beta)^2\]
sehingga residual besar akan dikuadratkan dan diberi bobot tidak proporsional — ini yang membuat OLS rentan terhadap outlier.
3.2 Identifikasi Outlier
3.2.1 Pemeriksaan Grafis
Scatterplot (untuk regresi sederhana) atau added-variable plot (untuk multivariat).
Residual plot untuk melihat pola non-acak dan titik residual ekstrem.
QQ-plot residual untuk melihat deviasi dari normalitas.
set.seed(123)
x <- 1:60
y <- 2 + 0.6*x + rnorm(length(x), 0, 3)
y[c(7, 35, 50)] <- y[c(7, 35, 50)] + c(18,-22, 25) # tiga outlier
ols <- lm(y~ x)
par(mfrow=c(1,3))
plot(x, y, pch=19, main="Scatterplot y vs x"); abline(ols, lwd=2)
plot(ols, which=1) # Residuals vs Fitted
qqnorm(rstandard(ols), main="QQ-plot Standardized Residual"); qqline(rstandard(ols))
3.2.2 Diagnostik Statistik
Notasi: \(\hat{e}_i\) residual OLS, \(\hat{\sigma}^2\) ragam residual, \(h_{ii}\) nilai diagonal hat matrix \(H = X(X^\top X)^{-1}X^\top\).
- Standardized/Studentized residual
\[r_i = \dfrac{\hat{e}_i}{\hat{\sigma}\,\sqrt{1 - h_{ii}}}, \quad t_i = \dfrac{\hat{e}_i}{\hat{\sigma}_{(i)}\,\sqrt{1 - h_{ii}}}.\]
Kriteria praktis: \(|r_i| > 2\) atau \(|t_i| > 3\) untuk kandidat outlier.
- Leverage
\[h_{ii} = x_i^\top (X^\top X)^{-1} x_i$, ambang populer: $h_{ii} > 2p\].
dengan 𝑝 jumlah parameter (termasuk intercept), 𝑛 ukuran sampel.
- Cook’s Distance
\[D_i = \dfrac{r_i^2}{p}\cdot \dfrac{h_{ii}}{1 - h_{ii}}.\] Aturan praktis: 𝐷𝑖 > 1 perlu ditinjau.
- DFFITS dan DFBETAS
\[\mathrm{DFFITS}_i \;=\; \dfrac{\hat{y}_i - \hat{y}_i^{(i)}}{\hat{\sigma}_{(i)}\sqrt{h_{ii}}}, \qquad \mathrm{DFBETAS}_{ij} \;=\; \dfrac{\hat{\beta}_j - \hat{\beta}_{j}^{(i)}}{SE_{(i)}(\hat{\beta}_j)}.\]
library(broom)
library(dplyr)
aug <- augment(ols) %>%
mutate(
h_ii = hatvalues(ols),
rstd = rstandard(ols),
cooks = cooks.distance(ols)
)
head(aug, 10)## # A tibble: 10 × 11
## y x .fitted .resid .hat .sigma .cooksd .std.resid h_ii rstd
## <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.919 1 3.05 -2.13 0.0650 5.59 5.49e-3 -0.397 0.0650 -0.397
## 2 2.51 2 3.65 -1.14 0.0618 5.60 1.49e-3 -0.213 0.0618 -0.213
## 3 8.48 3 4.26 4.22 0.0587 5.57 1.91e-2 0.784 0.0587 0.784
## 4 4.61 4 4.86 -0.249 0.0557 5.60 6.27e-5 -0.0461 0.0557 -0.0461
## 5 5.39 5 5.46 -0.0756 0.0528 5.60 5.46e-6 -0.0140 0.0528 -0.0140
## 6 10.7 6 6.07 4.68 0.0500 5.56 1.97e-2 0.865 0.0500 0.865
## 7 25.6 7 6.67 18.9 0.0474 4.97 3.03e-1 3.49 0.0474 3.49
## 8 3.00 8 7.27 -4.27 0.0448 5.57 1.45e-2 -0.787 0.0448 -0.787
## 9 5.34 9 7.88 -2.54 0.0424 5.59 4.83e-3 -0.467 0.0424 -0.467
## 10 6.66 10 8.48 -1.82 0.0400 5.59 2.33e-3 -0.334 0.0400 -0.334
## # ℹ 1 more variable: cooks <dbl>par(mfrow=c(1,2))
plot(aug$rstd, pch=19, main="Standardized Residual", ylab="r_i"); abline(h=c(-2,2), lty=2)
plot(aug$cooks, pch=19, main="Cook's Distance", ylab="D_i"); abline(h=1, lty=2)
Catatan: Gunakan beberapa metrik sekaligus untuk keputusan yang lebih andal (mis. kandidat pengaruh jika rstd besar dan \(h_{ii}\) tinggi dan \(D_i\) besar).
3.2.3 Strategi Penanganan Outlier
3.2.3.1 Verifikasi Kualitas Data
Koreksi data entry error atau coding error.
Konfirmasi konteks substantif (outlier bermakna versus noise).
3.2.3.2 Transformasi dan Rekayasa Fitur
Transformasi log, sqrt, Box–Cox/Yeo–Johnson untuk menstabilkan ragam dan mengurangi pengaruh ekstrem.
Binning atau winsorizing (mengganti nilai di luar persentil 1–99).
y_log <- log(pmax(y - min(y) + 1, 1)) # contoh log; pastikan argumen positif
par(mfrow = c(1, 2))
plot(x, y, pch = 19, main = "Asli")
plot(x, y_log, pch = 19, main = "Transformasi log(y)")
3.2.3.3 Model & Estimasi yang Tahan
- Weighted least squares (WLS): beri bobot kecil pada observasi dengan varians besar.
- Model robust (M-estimator, LAD, RLM-IRLS).
- Quantile regression: memusat pada median/kuantil lain, lebih stabil terhadap outlier pada mean.
- Model alternatif: GLM dengan link dan family sesuai sifat error (bila masalahnya bukan sekadar outlier).
3.2.3.4 Dokumentasi Keputusan
- Laporkan kriteria, titik yang terpengaruh, perubahan koefisien, serta implikasi substantif.
- Jangan sekadar menghapus data tanpa alasan yang dapat diaudit.
3.2.4 Ide & Metode Robust Regression
3.2.4.1 Mengapa Robust?
Fungsi kuadrat memberi pengaruh tak terbatas (unbounded influence) pada residual besar. Robust regression mengganti fungsi kerugian kuadrat dengan kerugian less than quadratic atau bounded, sehingga pengaruh outlier dibatasi.
3.2.4.2 Keluarga M‑Estimator
M‑estimasi meminimalkan
\[ \hat{\beta} \;=\; \arg\min_{\beta}\; \sum_{i=1}^{n} \rho\!\left(\frac{y_i - x_i^{\top}\beta}{\sigma}\right). \]
dengan fungsi loss 𝜌(⋅) tertentu dan fungsi influence 𝜓 = 𝜌′. Dua fungsi populer:
- Huber loss dengan ambang 𝑐:
$$ - Tukey’s bisquare (redescending):
\[ \rho(e) = \begin{cases} \dfrac{c^{2}}{6}\left[1 - \left(1 - \left(\dfrac{e}{c}\right)^{2}\right)^{3}\right], & |e|\le c,\\[8pt] \dfrac{c^{2}}{6}, & |e|>c, \end{cases} \qquad \psi(e) = \begin{cases} e\left(1 - \left(\dfrac{e}{c}\right)^{2}\right)^{2}, & |e|\le c,\\[8pt] 0, & |e|>c. \end{cases} \] Estimasi praktis sering menggunakan IRLS (Iteratively Reweighted Least Squares) dengan bobot
\[ w_i = \frac{\psi(e_i)}{e_i}. \] sehingga solusi pada iterasi \(t\) adalah pemecahan WLS dengan bobot \(w_i^{(t)}\).
3.2.4.3 Menentukan nilai \(c\)**
Nilai \(c\) adalah tuning constant yang menentukan kapan sebuah residual dianggap outlier.
Jika \(|e| \le c\), maka residual dianggap normal → diperlakukan seperti OLS.
Jika \(|e| > c\), maka residual dianggap outlier → diberi pengaruh terbatas.
Cara memilih \(c\)
- Huber Loss
Nilai \(c \approx 1.345\,\hat{\sigma}\)
sering digunakan karena menghasilkan estimator dengan efisiensi ~95%
dibanding OLS saat error normal.
- Tukey’s Bisquare
Nilai \(c \approx 4.685\,\hat{\sigma}\)
banyak dipakai agar estimator efisien mendekati OLS di distribusi
normal.
Artinya: \(c\) dikalibrasi berdasarkan skala error (\(\hat{\sigma}\)), supaya fungsi robust tetap efisien kalau data normal, tapi membatasi outlier.
Fungsi influence
Fungsi influence (\(\psi(e)\)) adalah
turunan dari fungsi loss \(\rho(e)\):
\[
\psi(e) = \frac{d}{de}\,\rho(e).
\] Interpretasi
- \(\psi(e)\) menunjukkan seberapa besar pengaruh sebuah residual terhadap estimasi koefisien.
- Pada OLS: \(\psi(e)=e\). Residual makin besar → pengaruh makin besar (tanpa batas).
- Pada Huber dan Tukey: \(\psi(e)\) dibatasi. Residual besar tidak bisa “mengendalikan” hasil.
Contoh - Huber: \[ \psi(e) = \begin{cases} e, & |e|\le c,\\[6pt] c\,\operatorname{sign}(e), & |e|>c \end{cases} \]
Efek outlier dibatasi pada \(\pm c\).
Tukey: \[ \psi(e) = \begin{cases} e\left(1 - \left(\dfrac{e}{c}\right)^2\right)^2, & |e|\le c,\\[8pt] 0, & |e|>c \end{cases} \]
outlier ekstrem bahkan diberi bobot nol.
Bobot dalam Iteratively Reweighted Least Squares (IRLS)
Dalam IRLS, setiap iterasi menghitung bobot \(w_i\) berdasarkan residual saat itu:
\[ w_i \;=\; \dfrac{\psi(e_i)}{e_i} \] Sehingga regresi robust menjadi Weighted Least Squares (WLS) dengan bobot berbeda di tiap observasi.
Contoh Perhitungan Bobot
Misal residual dari model OLS awal adalah: \(e=\{-0.5,\; 1.2,\; 5.0\}\), \(c=1.345\).
Huber loss:
- Untuk \(e = -0.5\): \(|e| < c \Rightarrow \psi(e) = -0.5 \Rightarrow w = 1.\)
- Untuk \(e = 1.2\): \(|e| < c \Rightarrow \psi(e) = 1.2 \Rightarrow w = 1.\)
- Untuk \(e = 5.0\): \(|e| > c \Rightarrow \psi(e) = 1.345 \Rightarrow w = 1.345/5 = 0.269.\)
Observasi ke-3 (outlier) hanya dapat bobot \(0.269\), jauh lebih kecil dibanding observasi lain.
Iterasi IRLS
- Fit OLS \(\rightarrow\) dapat residual awal.
- Hitung \(w_i\) untuk semua observasi.
- Fit ulang dengan WLS (berat \(w_i\)).
- Update residual \(\rightarrow\) hitung bobot baru \(\rightarrow\) ulangi.
- Berhenti saat perubahan koefisien kecil (konvergen).
Contoh R (Huber dengan IRLS)
library(MASS)
# Data dengan outlier
set.seed(123)
x <- 1:20
y <- 2*x + rnorm(20, 0, 2)
y[c(5, 15)] <- y[c(5, 15)] + 20 # tambahkan outlier
# Robust regression (Huber M-estimator)
fit_rob <- rlm(y~ x, psi = psi.huber)
# Lihat bobot yang diberikan
weights(fit_rob)## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1Output weights(fit_rob) akan menunjukkan:
- Hampir semua observasi bobot ~1
- Observasi outlier → bobot jauh lebih kecil (<1)
Estimasi Skala \(\hat{\sigma}\)
\(\hat{\sigma}\) pada robust regression
tidak sama dengan standar deviasi OLS.
Karena OLS sensitif terhadap outlier, dipakai estimasi skala robust.
Median Absolute Deviation (MAD)
\[
\hat{\sigma} = 1.4826 \cdot \operatorname{median}\bigl(\,|e_i -
\operatorname{median}(e)|\,\bigr)
\]
Konstanta 1.4826 membuatnya konsisten dengan simpangan baku normal.
Update Iteratif dalam IRLS
Setelah setiap iterasi, skala \(\sigma\) diupdate dengan formula berbasis
\(\psi(e)\).
Implementasi di R (fungsi rlm() dari paket
MASS):
- Menggunakan MAD residual OLS sebagai starting value.
- Update \(\hat{\sigma}\) dengan metode robust (Huber’s Proposal 2).
## [1] 2.487055##
## Call: rlm(formula = y ~ x, psi = psi.huber)
## Residuals:
## Min 1Q Median 3Q Max
## -4.6241 -1.6307 -0.0158 1.7847 19.3656
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 0.9707 1.2954 0.7493
## x 1.9845 0.1081 18.3511
##
## Residual standard error: 2.487 on 18 degrees of freedomCatatan:
- Nilai \(c\): ambang residual dianggap outlier (biasanya proporsional dengan \(\hat{\sigma}\)).
- Influence function: ukuran pengaruh residual, dibatasi dalam robust regression.
- Bobot IRLS: membatasi kontribusi outlier, membuat model stabil.
- \(\hat{\sigma}\): estimasi skala robust (MAD atau update IRLS), bukan simpangan baku OLS.
3.2.4.4 LAD (Least Absolute Deviation)
Meminimalkan: \[ \min_{\beta}\ \sum_{i=1}^{n} \big|\, y_i - x_i^{\top}\beta \,\big| \]
- Lebih tahan terhadap outlier pada \(Y\), namun kurang efisien jika error
benar-benar normal.
- Dapat diselesaikan dengan pemrograman linear atau algoritme khusus.
3.2.4.5 Quantile Regression
Menggeneralisasi LAD ke kuantil \(\tau \in (0, 1)\) dengan meminimalkan check loss \[ \sum_{i=1}^{n} \rho_\tau(\varepsilon_i), \qquad \rho_\tau(u) = u \cdot \big(\tau - \mathbf{1}\{u < 0\}\big). \] Memungkinkan analisis efek heterogen di berbagai titik distribusi respons (median, P90, dst.).
3.2.5 Praktik di R: OLS vs Robust
library(MASS) # rlm (Huber & bisquare)
library(quantreg) # quantile regression
library(dplyr); library(broom)3.2.5.1 Simulasi & OLS
set.seed(42)
n <- 120
x1 <- runif(n, 0, 10)
x2 <- rnorm(n)
y <- 3 + 2*x1- 1.5*x2 + rnorm(n, 0, 2)
# Tambah outlier respons & leverage
y[c(10, 55)] <- y[c(10, 55)] + c(25,-30)
x1[c(80, 95)] <- x1[c(80, 95)] + c(15, 18)
dat <- data.frame(y, x1, x2)
fit_ols <- lm(y~ x1 + x2, dat)
glance(fit_ols)## # A tibble: 1 × 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.476 0.467 5.53 53.1 3.96e-17 2 -374. 756. 767.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>Diagnosis pengaruh
h <- hatvalues(fit_ols)
cooks <- cooks.distance(fit_ols)
rstd <- rstandard(fit_ols)
cand_idx <- which( (abs(rstd) > 2) | (h > 2*length(coef(fit_ols))/n) | (cooks > 1) )
length(cand_idx); cand_idx## [1] 7## 10 44 55 58 68 80 95
## 10 44 55 58 68 80 953.2.5.2 Robust (Huber & Bisquare)
fit_hub <- rlm(y ~ x1 + x2, data = dat, psi = psi.huber) # Huber
fit_bis <- rlm(y ~ x1 + x2, data = dat, psi = psi.bisquare) # Tukey bisquare
bind_rows(
broom::tidy(fit_ols) %>% mutate(model = "OLS"),
broom::tidy(fit_hub) %>% mutate(model = "Huber"),
broom::tidy(fit_bis) %>% mutate(model = "Bisquare")
) %>%
arrange(term, model)## # A tibble: 9 × 6
## term estimate std.error statistic p.value model
## <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 (Intercept) 2.82 0.353 8.00 NA Bisquare
## 2 (Intercept) 3.22 0.360 8.92 NA Huber
## 3 (Intercept) 5.53 0.926 5.97 2.62e- 8 OLS
## 4 x1 2.04 0.0532 38.2 NA Bisquare
## 5 x1 1.95 0.0543 35.9 NA Huber
## 6 x1 1.44 0.140 10.3 4.24e-18 OLS
## 7 x2 -1.64 0.215 -7.63 NA Bisquare
## 8 x2 -1.51 0.220 -6.88 NA Huber
## 9 x2 -0.605 0.564 -1.07 2.86e- 1 OLSPerbandingan garis regresi (kasus sederhana y~x1):
fit_ols1 <- lm(y~ x1, dat)
fit_hub1 <- rlm(y~ x1, dat, psi=psi.huber)
fit_bis1 <- rlm(y~ x1, dat, psi=psi.bisquare)
plot(dat$x1, dat$y, pch=19, main="OLS vs Robust (y ~ x1)", xlab="x1", ylab="y")
abline(fit_ols1, lwd=2, lty=2)
abline(fit_hub1, lwd=2)
abline(fit_bis1, lwd=2, lty=3)
legend("topleft", c("OLS","Huber","Bisquare"), lty=c(2,1,3), lwd=2, bty="n")
3.2.5.3 Quantile Regression (Median)
## # A tibble: 3 × 5
## term estimate conf.low conf.high tau
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.22 2.61 4.16 0.5
## 2 x1 1.97 1.74 2.08 0.5
## 3 x2 -1.50 -1.68 -0.924 0.53.2.5.4 Ringkasan Diagnostik & Rekomendasi
dx <- tibble(
idx = 1:n,
rstd = rstd,
leverage = h,
cooks = cooks
)
flagged <- (abs(rstd) > 2) | (h > 2 * length(coef(fit_ols)) / n) | (cooks > 1)
dx %>%
filter(flagged) %>%
arrange(desc(abs(rstd)))## # A tibble: 7 × 4
## idx rstd leverage cooks
## <int> <dbl> <dbl> <dbl>
## 1 55 -5.82 0.0496 0.589
## 2 95 -5.15 0.315 4.07
## 3 80 -4.64 0.0656 0.504
## 4 10 4.12 0.0151 0.0869
## 5 58 -0.768 0.0999 0.0218
## 6 68 0.725 0.0571 0.0106
## 7 44 -0.157 0.0528 0.000456Tahapan Analisis
- Eksplorasi grafis: scatter/QQ/residual plots.
- Hitung metrik: \(r_i\), \(h_{ii}\), \(D_i\), DFFITS/DFBETAS.
- Verifikasi data: benahi kesalahan input.
- Modelisasi: bandingkan OLS vs robust
(Huber/Bisquare/LAD/Quantile).
- Sensitivitas: laporkan perubahan koefisien/prediksi
setelah menandai/menangani outlier.
- Dokumentasi: catat kriteria, ambang, dan alasan keputusan.
Template Pelaporan - Data &
konteks: sumber, pembersihan, asumsi.
- Identifikasi outlier: metrik yang digunakan dan titik
yang ditandai.
- Penanganan: transformasi, winsorizing, atau model
robust; alasan pemilihan.
- Hasil: koefisien, interval kepercayaan/kuantil,
akurasi/predictive performance.
- Uji sensitivitas: perbandingan OLS vs robust.
- Catatan keterbatasan: potensi bias, ukuran sampel,
generalisasi.
Catatan Praktis - Gunakan beberapa indikator
bersamaan; jangan andalkan satu metrik.
- Ambang rule-of-thumb (mis. \(|r_i|>2\), \(D_i>1\)) bersifat panduan, bukan aturan
keras.
- Outlier dapat menyimpan informasi penting; pahami konteks sebelum
menghapus.
- Robust \(\neq\) kebal terhadap semua
masalah; tetap cek spesifikasi model, heteroskedastisitas, dan
multikolinearitas.
Referensi - Rousseeuw, P. J., & Leroy, A. M.
(2003). Robust Regression and Outlier Detection. Wiley.
- Huber, P. J., & Ronchetti, E. M. (2009). Robust
Statistics. Wiley.
- Fox, J. (2016). Applied Regression Analysis and Generalized Linear
Models. Sage.
- Koenker, R. (2005). Quantile Regression. Cambridge University
Press.
3.3 Tugas
Petunjuk
- Bagian A dan B dapat dijawab di kertas/jurnal Anda.
- Bagian C memerlukan R. Semua kode disediakan dan dapat langsung dijalankan.
- Kunci jawaban ada di bagian akhir dokumen.
3.3.1 Bagian A — Konsep (Jawaban Singkat)
- Jelaskan apa itu outlier pada regresi linier. Bedakan response outlier, predictor outlier, leverage point, dan influential point.
- Mengapa OLS sensitif terhadap outlier? Jelaskan memakai fungsi kerugian OLS.
- Jelaskan apa yang dimaksud influence function 𝜓(𝑒) dalam robust regression dan bedanya pada OLS, Huber, dan Tukey’s bisquare.
- Bagaimana memilih nilai ambang 𝑐 pada Huber dan Tukey? Jelaskan peran ̂𝜎 (skala robust).
jawab
1 -Outlier: observasi yang
menyimpang jauh dari pola data umum.
- Response outlier: nilai residu sangat besar
(penyimpangan pada variabel dependen).
- Predictor outlier: nilai pada variabel independen
ekstrem dibanding observasi lain.
- Leverage point: observasi dengan nilai variabel
independen jauh dari rata-rata, sehingga berpotensi memberi pengaruh
besar pada garis regresi.
- Influential point: observasi yang jika dihilangkan,
akan mengubah hasil estimasi regresi secara signifikan (misalnya
terdeteksi dengan Cook’s Distance).
2 OLS meminimalkan Sum of Squared Errors (SSE):
\[ SSE = \sum_{i=1}^n (y_i - \hat{y}_i)^2 \]
Karena kuadrat digunakan, residu besar (outlier) akan memiliki pengaruh yang kuat terhadap estimasi parameter. Akibatnya, satu outlier saja dapat menggeser garis regresi secara drastis.
3 - OLS: influence function \(\psi(e) = e\). Linear, sehingga outlier
punya pengaruh besar.
- Huber: kombinasi linear dan konstan.
\[
\psi(e) = \begin{cases} e & |e| \leq c \\ c \cdot \text{sign}(e)
& |e| > c \end{cases}
\] → Outlier besar tidak terus-menerus memperbesar
pengaruhnya.
- Tukey’s bisquare: membatasi pengaruh lebih ketat.
Untuk \(|e| > c\), \(\psi(e)=0\). Outlier jauh sama sekali tidak
berpengaruh.
4 - Nilai \(c\)
dipilih berdasarkan keseimbangan antara efisiensi (untuk data normal)
dan ketahanan terhadap outlier.
- Umumnya, \(c = 1.345 \times
\hat{\sigma}\) (Huber) digunakan agar efisiensi ~95% pada data
normal.
- Pada Tukey, nilai \(c\) menentukan
batas residu yang dianggap sebagai outlier (misalnya \(c=4.685\)).
- Peran c: sebagai tuning constant yang
mengontrol sensitivitas metode terhadap outlier.
3.3.2 Bagian B — Analisis (Teori)
5) Bobot Huber dan identifikasi outlier
Diberikan residual OLS: \(e = \{-0.5, 1.2, 5.0\}\) dan \(c = 1.345\).
- Hitung bobot \(w_i\) untuk Huber loss.
- Observasi mana yang dianggap outlier oleh pendekatan ini?
jawab
- Hitung bobot w untuk Huber loss:
resid <- c(-0.5, 1.2, 5.0)
c_val <- 1.345
# Fungsi bobot Huber
w_huber <- function(e, c){
ifelse(abs(e) <= c, 1, c/abs(e))
}
weights <- w_huber(resid, c_val)
data.frame(residual = resid, huber_weight = weights)## residual huber_weight
## 1 -0.5 1.000
## 2 1.2 1.000
## 3 5.0 0.269- Interpretasi:
- Residual -0.5 → dalam batas, bobot ≈ 1.
- Residual 1.2 → dalam batas, bobot ≈ 1.
- Residual 5.0 → outlier, bobot < 1 → pengaruh dikurangi.
- Residual -0.5 → dalam batas, bobot ≈ 1.
Observasi mana yang dianggap outlier? → Residual 5.0.
6) Cook’s Distance Diberikan Cook’s Distance untuk tiga observasi: \(\{0.2, 0.8, 1.5\}\).
- Mana yang paling berpengaruh?
- Mengapa Cook’s Distance berguna dalam deteksi observasi berpengaruh?
jawab
- Aturan umum: observasi dengan \(D_i > 1\) dianggap berpengaruh (influential).
cooksD <- c(0.2, 0.3, 1.5)
influential <- cooksD > 1
data.frame(observation = 1:3, cooksD = cooksD, influential = influential)## observation cooksD influential
## 1 1 0.2 FALSE
## 2 2 0.3 FALSE
## 3 3 1.5 TRUE- Mana yang paling berpengaruh? → Observasi ke-3 (D =
1.5).
- Mengapa Cook’s Distance berguna?
Karena menggabungkan informasi dari leverage (X) dan residual (Y) sehingga mendeteksi observasi yang benar-benar mengubah hasil estimasi regresi.
3.3.3 Bagian C — Praktik dengan R
C.1 Simulasi data dengan outlier
set.seed(123)
x <- 1:20
y <- 2*x + rnorm(20, 0, 2)
y[c(5, 15)] <- y[c(5, 15)] + 20 # tambahkan outlier
dat <- data.frame(y, x)Estimasi OLS vs Robust (Huber)
library(MASS) # rlm
fit_ols <- lm(y~ x, dat)
fit_rob <- rlm(y~ x, dat, psi = psi.huber)
coef(fit_ols)## (Intercept) x
## 2.935974 1.937836## (Intercept) x
## 0.9707001 1.9844546C.3 Visualisasi garis OLS vs Robust
plot(dat$x, dat$y, pch=19, xlab="x", ylab="y", main="OLS vs Robust (Huber)")
abline(fit_ols, lwd=2, lty=2, col="gray40")
abline(fit_rob, lwd=2, col="black")
legend("topleft", c("OLS","Robust (Huber)"),
lty=c(2,1), lwd=2, bty="n")
C.4 Diagnostik standar (leverage, residual, Cook’s D)
h_ii <- hatvalues(fit_ols)
rstd <- rstandard(fit_ols)
cooks <- cooks.distance(fit_ols)
idx_flag <- which( (abs(rstd)>2) | (h_ii > 2*length(coef(fit_ols))/nrow(dat)) | (cooks > 1) )
data.frame(idx=1:nrow(dat), rstd=round(rstd,3), h_ii=round(h_ii,3), cooks=round(cooks,3),
flagged = 1:nrow(dat) %in% idx_flag)[idx_flag, ]## idx rstd h_ii cooks flagged
## 5 5 2.906 0.095 0.446 TRUE
## 15 15 2.760 0.080 0.333 TRUEC.5 Bobot IRLS dari model robust
w <- weights(fit_rob)
data.frame(idx=1:length(w), weight=round(w,3))[order(w), ][1:5, ] # tampilkan 5 bobot terkecil## idx weight
## 1 1 1
## 2 2 1
## 3 3 1
## 4 4 1
## 5 5 13.4 Quis
- Waktu pengerjaan: 15–20 menit.
- Tuliskan jawaban singkat dan jelas.
- Bagian Kunci Jawaban & Rubrik dapat disembunyikan saat dibagikan ke mahasiswa.
3.4.1 Bagian A — Pilihan Ganda (PG)
(Masing-masing 2 poin, total 10 poin)
A1. Mengapa OLS sangat sensitif terhadap outlier?
A. Menggunakan fungsi absolute residual.
B. Menggunakan fungsi kuadrat residual.
C. Menggunakan median residual.
D. Tidak menggunakan residual.
A2. Dalam regresi, leverage point adalah …
A. Observasi dengan nilai Y ekstrem.
B. Observasi dengan nilai X jauh dari pusat
distribusi.
C. Observasi dengan \(|r_i| >
2\).
D. Observasi dengan Cook’s \(D >
1\).
A3. Jika residual \(e_i = 5\) dan ambang Huber \(c = 1.345\), maka bobot \(w_i\) kira-kira …
A. 1
B. 0.269
C. 5
D. 1.345
A4. Pada Tukey’s bisquare, apa yang terjadi jika \(|e| > c\)?
A. Residual dihitung penuh.
B. Residual diabaikan (bobot nol).
C. Residual diganti dengan rata-rata.
D. Residual dibagi dua.
A5. Estimasi skala \(\hat{\sigma}\) yang tahan outlier biasanya menggunakan …
A. Rata-rata residual.
B. SD residual OLS.
C. Median Absolute Deviation (MAD).
D. Variansi OLS.
3.4.2 Bagian B — Isian Singkat
(Masing-masing 5 poin, total 15 poin)
B6. Tuliskan formula Cook’s Distance dan jelaskan arti jika \(D_i > 1\).
jawab
Formula (bentuk residual-based):
\[ D_i \;=\; \frac{( \hat{y} - \hat{y}_{(i)} )^\top ( \hat{y} - \hat{y}_{(i)} )}{p\,\hat{\sigma}^2} \]
di mana \(\) adalah vektor prediksi menggunakan semua observasi, \(_{(i)}\) adalah vektor prediksi ketika observasi ke-\(i\) dihapus, \(p\) adalah jumlah parameter (termasuk intercept), dan \(^2\) adalah estimasi varians residual (MSE).
Formula (ekivalen, menggunakan residual baku dan leverage):
\[ D_i \;=\; \frac{e_i^2}{p\,\hat{\sigma}^2}\cdot\frac{h_{ii}}{(1-h_{ii})^2} \]
dengan \(e_i\) residual terestimasi untuk observasi ke-\(i\) dan \(h_{ii}\) diagonal dari hat matrix \(H=X(X^\top X)^{-1}X^\top\) (leverage).
Interpretasi:
- Nilai \(D_i\) mengukur seberapa besar pengaruh observasi ke-\(i\)
terhadap hasil estimasi regresi (perubahan pada vektor prediksi jika
observasi dihilangkan).
- Aturan praktis sering digunakan: jika \(D_i > 1\) maka observasi
dianggap berpengaruh (influential). Namun aturan ini bersifat heuristik;
alternatif: bandingkan dengan kuantil F (misal \(D_i > \frac{4}{n}\)
atau gunakan cutoff lain bergantung konteks).
- Observasi dengan nilai Cook’s distance besar harus diperiksa: apakah
data salah entry, outlier, atau memang observasi sah yang menunjukkan
fenomena penting? Setelah diperiksa, pertimbangkan transformasi, model
robust, atau analisis sensitifitas (fit model tanpa observasi
tersebut).
B7. Sebutkan dua metode robust regression selain Huber loss.
jawab
Beberapa metode/regresi robust selain Huber: - Tukey’s
Biweight (bisquare) (M-estimator dengan fungsi psi bisquare) —
mengurangi bobot outlier kuat (memberi bobot nol untuk residu sangat
besar).
- Least Trimmed Squares (LTS) — meminimalkan jumlah
kuadrat untuk subset data (trimmed), tahan terhadap banyak
outlier.
- S-estimator, MM-estimator — estimasi robust dengan
breakdown point tinggi.
- RANSAC — algoritma iteratif yang menemukan model yang
cocok untuk subset data inliers.
Di bawah ini contoh implementasi di R menggunakan paket
MASS untuk Tukey’s biweight dan
robustbase untuk LTS.
# Jika belum terinstal, uncomment baris berikut:
# install.packages(c("MASS","robustbase"))
library(MASS) # rlm()
library(robustbase) # ltsReg()
# Simulasi data contoh
set.seed(1)
x <- 1:30
y <- 2 + 0.5*x + rnorm(30,0,1)
# tambahkan beberapa outlier
y[c(5,15,28)] <- y[c(5,15,28)] + c(8,-6,10)
data <- data.frame(x,y)
# Ordinary least squares
lm_fit <- lm(y ~ x, data=data)
# Robust: Tukey's biweight via rlm (psi = psi.bisquare)
rlm_bisq <- rlm(y ~ x, data=data, psi = psi.bisquare, maxit=100)
# Robust: Least Trimmed Squares (LTS)
lts_fit <- ltsReg(y ~ x, data=data)
summary(lm_fit)##
## Call:
## lm(formula = y ~ x, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3536 -0.9315 -0.1856 0.2953 7.9477
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.35959 0.95000 2.484 0.0193 *
## x 0.50793 0.05351 9.492 3.01e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.537 on 28 degrees of freedom
## Multiple R-squared: 0.7629, Adjusted R-squared: 0.7544
## F-statistic: 90.09 on 1 and 28 DF, p-value: 3.005e-10##
## Call: rlm(formula = y ~ x, data = data, psi = psi.bisquare, maxit = 100)
## Residuals:
## Min 1Q Median 3Q Max
## -5.0260 -0.6060 0.1477 0.6127 8.3501
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 2.1183 0.3510 6.0353
## x 0.5022 0.0198 25.4004
##
## Residual standard error: 0.9455 on 28 degrees of freedom## Intercept x
## 2.1714993 0.5063499Ringkasan perbedaan praktis:
- Huber menurunkan bobot outlier, tetapi masih memberikan kontribusi
linear di luar batas; Tukey (bisquare) memberikan bobot nol pada residu
yang sangat besar (lebih agresif dalam menolak outlier).
- LTS memiliki high breakdown point (tahan terhadap banyak
outlier) tetapi tidak semulus M-estimators dalam hal efisiensi jika data
sebenarnya normal.
B8. Diberikan residual \(\{-0.5, 1.2, 5.0\}\) dengan \(c = 1.345\). Tentukan residual mana yang dianggap outlier oleh metode Huber dan berikan alasannya.
jawab
Prinsip Huber:
Fungsi psi Huber dengan parameter cut-off \(c\) adalah \[
\psi_c(r) = \begin{cases}
r & \text{jika } |r| \le c,\\[4pt]
c\cdot \operatorname{sign}(r) & \text{jika } |r| > c.
\end{cases}
\]
Dengan kata lain, Huber memperlakukan residu kecil secara linear (seperti OLS) tetapi membatasi kontribusi residu besar ke konstanta \(c\) (mengurangi pengaruhnya). Seringkali observasi dengan \(|r| > c\) dipandang sebagai outlier karena kontribusi liniernya dipotong.
Untuk vektor residual yang diberikan, hitung absolutnya terhadap \(c=1.345\):
- \(|-0.5| = 0.5 \) → bukan outlier.
- \(|1.2| = 1.2 \) → bukan outlier.
- \(|5.0| = 5.0 > 1.345\) → dianggap outlier (pengaruhnya dibatasi oleh Huber).
Di bawah ini kode R untuk menghitung fungsi psi dan bobot Huber:
r <- c(-0.5, 1.2, 5.0)
c <- 1.345
psi_huber <- function(r, c) {
ifelse(abs(r) <= c, r, c * sign(r))
}
# bobot = psi(r)/r (untuk r != 0)
weights_huber <- ifelse(r==0, 1, psi_huber(r,c)/r)
data.frame(r = r, psi = psi_huber(r,c), weight = weights_huber)## r psi weight
## 1 -0.5 -0.500 1.000
## 2 1.2 1.200 1.000
## 3 5.0 1.345 0.269Penjelasan hasil: residu ke-3 (5.0) memiliki \(|r|>c\), sehingga psi di-truncate menjadi \(c\cdot(r)=1.345\), dan bobotnya menjadi kecil (psi/r ≈ 1.345/5 = 0.269) — artinya observasi tersebut diberikan bobot lebih kecil dalam estimasi, sehingga tidak mendominasi parameter model. Residual pertama dan kedua diberi bobot 1 (tidak dikurangi).