1.1 Konsep Analisis Regresi dan Korelasi

Data

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

Korelasi

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

Regresi

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

Visualisasi

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

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

1.9.1 Persiapan Data dan Eksplorasi

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

1.9.2 Estimasi OLS & Ringkasan

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

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

diag_df <- tibble(
fitted = fitted(fit1), resid = resid(fit1),
stdres = rstandard(fit1), hat = hatvalues(fit1), cook = cooks.distance(fit1)
)

p1 <- ggplot(diag_df, aes(fitted, resid)) +
geom_point(alpha = 0.7) +
geom_hline(yintercept = 0, linetype = 2) +
labs(x = "Fitted", y = "Residuals")

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.9.4 Multikolinearitas (VIF)

car::vif(fit1)
##       wt       hp     disp 
## 4.844618 2.736633 7.324517

1.9.5 Robust SE (White HC3) & Interpretasi

coeftest(fit1) # SE OLS
## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) 37.10550527  2.11081525 17.5788 < 2.2e-16 ***
## wt          -3.80089058  1.06619064 -3.5649  0.001331 ** 
## hp          -0.03115655  0.01143579 -2.7245  0.010971 *  
## disp        -0.00093701  0.01034974 -0.0905  0.928507    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(fit1, vcov = vcovHC(fit1, type = "HC3")) # SE robust (HC3)
## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) 37.10550527  2.71026981 13.6907 6.257e-14 ***
## wt          -3.80089058  1.08086207 -3.5165   0.00151 ** 
## hp          -0.03115655  0.01437457 -2.1675   0.03886 *  
## disp        -0.00093701  0.01012138 -0.0926   0.92690    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.9.6 Nonlinearitas: Polinomial vs Spline pada wt

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

1.9.7 Validasi Silang (K-fold) untuk Perbandingan Model

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

1.10 Studi Kasus 2: Simulasi Heteroskedastisitas (WLS vs OLS vs Robust)

n <- 400
x1 <- runif(n, 0, 10)
x2 <- rnorm(n, 5, 2)
sigma_i <- 0.6 + 0.2 * x1 # var naik dengan x1
err <- rnorm(n, 0, sigma_i)
y <- 2 + 1.5*x1- 1.2*x2 + err
sim <- tibble(y, x1, x2, w = 1/sigma_i^2)

fit_ols <- lm(y~ x1 + x2, data = sim)
fit_wls <- lm(y~ x1 + x2, data = sim, weights = w)

# Tabel perbandingan SE
comp_se <- tibble(
term = names(coef(fit_ols)),
OLS_SE = coef(summary(fit_ols))[, "Std. Error"],
Robust_SE = sqrt(diag(vcovHC(fit_ols, type = "HC3"))),
WLS_SE = coef(summary(fit_wls))[, "Std. Error"]
)
comp_se
# Uji Breusch-Pagan untuk heteroskedastisitas
lmtest::bptest(fit_ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  fit_ols
## BP = 29.415, df = 2, p-value = 4.098e-07
ggplot(data.frame(fitted=fitted(fit_ols), resid=resid(fit_ols)), aes(fitted, resid)) +
geom_point(alpha=.6) + geom_smooth(se=FALSE, method="loess") +
geom_hline(yintercept = 0, linetype=2) + labs(x="Fitted", y="Residuals")

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

1.12.1 Simulasi data dengan galat AR(1)

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

dar1 <- tibble::tibble(t = 1:T, x, y)
head(dar1)

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

library(nlme)
# OLS biasa
fit_ols_ar1 <- lm(y~ x, data = dar1)

# GLS dengan korelasi AR(1) di residual
fit_gls_ar1 <- gls(y~ x, data = dar1,
correlation = corAR1(form =~ t))
summary(fit_ols_ar1)
## 
## Call:
## lm(formula = y ~ x, data = dar1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3360 -0.8231  0.0026  0.8211  3.9013 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.95441    0.06043   15.79   <2e-16 ***
## x            2.06045    0.05162   39.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.209 on 398 degrees of freedom
## Multiple R-squared:  0.8002, Adjusted R-squared:  0.7997 
## F-statistic:  1594 on 1 and 398 DF,  p-value: < 2.2e-16
summary(fit_gls_ar1)
## Generalized least squares fit by REML
##   Model: y ~ x 
##   Data: dar1 
##        AIC      BIC    logLik
##   1139.541 1155.487 -565.7706
## 
## Correlation Structure: AR(1)
##  Formula: ~t 
##  Parameter estimate(s):
##       Phi 
## 0.5772332 
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 0.9552097 0.11651869  8.19791       0
## x           2.0449090 0.03743358 54.62766       0
## 
##  Correlation: 
##   (Intr)
## x 0.004 
## 
## Standardized residuals:
##          Min           Q1          Med           Q3          Max 
## -2.752477543 -0.674750151 -0.007333444  0.679546648  3.225627947 
## 
## Residual standard error: 1.210605 
## Degrees of freedom: 400 total; 398 residual

1.12.3 Diagnostik serial correlation

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

par(mfrow = c(1,1))
lmtest::dwtest(fit_ols_ar1) # indikasi korelasi serial pada OLS
## 
##  Durbin-Watson test
## 
## data:  fit_ols_ar1
## DW = 0.84993, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

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

# Estimasi rho dari Regresi resid_t ~ resid_{t-1}
res <- resid(fit_ols_ar1)
rho_hat <- coef(lm(res[-1]~ 0 + res[-length(res)]))[1]

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

1.13 Latihan

set.seed(123)

# Paket yang Diperlukan
pkgs <- c("dplyr","ggplot2","car","lmtest","sandwich","MASS","boot","broom")
new <- pkgs[!(pkgs %in% installed.packages()[,"Package"])]
if(length(new)) install.packages(new, dependencies = TRUE)
lapply(pkgs, library, character.only = TRUE)
## [[1]]
##  [1] "nlme"      "splines"   "modelr"    "car"       "carData"   "lmtest"   
##  [7] "zoo"       "sandwich"  "broom"     "lubridate" "forcats"   "stringr"  
## [13] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [19] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [25] "methods"   "base"     
## 
## [[2]]
##  [1] "nlme"      "splines"   "modelr"    "car"       "carData"   "lmtest"   
##  [7] "zoo"       "sandwich"  "broom"     "lubridate" "forcats"   "stringr"  
## [13] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [19] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [25] "methods"   "base"     
## 
## [[3]]
##  [1] "nlme"      "splines"   "modelr"    "car"       "carData"   "lmtest"   
##  [7] "zoo"       "sandwich"  "broom"     "lubridate" "forcats"   "stringr"  
## [13] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [19] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [25] "methods"   "base"     
## 
## [[4]]
##  [1] "nlme"      "splines"   "modelr"    "car"       "carData"   "lmtest"   
##  [7] "zoo"       "sandwich"  "broom"     "lubridate" "forcats"   "stringr"  
## [13] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [19] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [25] "methods"   "base"     
## 
## [[5]]
##  [1] "nlme"      "splines"   "modelr"    "car"       "carData"   "lmtest"   
##  [7] "zoo"       "sandwich"  "broom"     "lubridate" "forcats"   "stringr"  
## [13] "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"    "ggplot2"  
## [19] "tidyverse" "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [25] "methods"   "base"     
## 
## [[6]]
##  [1] "MASS"      "nlme"      "splines"   "modelr"    "car"       "carData"  
##  [7] "lmtest"    "zoo"       "sandwich"  "broom"     "lubridate" "forcats"  
## [13] "stringr"   "dplyr"     "purrr"     "readr"     "tidyr"     "tibble"   
## [19] "ggplot2"   "tidyverse" "stats"     "graphics"  "grDevices" "utils"    
## [25] "datasets"  "methods"   "base"     
## 
## [[7]]
##  [1] "boot"      "MASS"      "nlme"      "splines"   "modelr"    "car"      
##  [7] "carData"   "lmtest"    "zoo"       "sandwich"  "broom"     "lubridate"
## [13] "forcats"   "stringr"   "dplyr"     "purrr"     "readr"     "tidyr"    
## [19] "tibble"    "ggplot2"   "tidyverse" "stats"     "graphics"  "grDevices"
## [25] "utils"     "datasets"  "methods"   "base"     
## 
## [[8]]
##  [1] "boot"      "MASS"      "nlme"      "splines"   "modelr"    "car"      
##  [7] "carData"   "lmtest"    "zoo"       "sandwich"  "broom"     "lubridate"
## [13] "forcats"   "stringr"   "dplyr"     "purrr"     "readr"     "tidyr"    
## [19] "tibble"    "ggplot2"   "tidyverse" "stats"     "graphics"  "grDevices"
## [25] "utils"     "datasets"  "methods"   "base"

1. Simulasi Data

n <- 200
X1 <- rnorm(n, 0, 1)
# X2 berkorelasi ~0.7 dengan X1
rho <- 0.7
X2 <- rho*X1 + sqrt(1-rho^2)*rnorm(n)
X3 <- rnorm(n)

# Heteroskedastisitas: sd error tergantung |X1|
eps <- rnorm(n, sd = 0.5 + 0.5*abs(X1))

# Model benar (ground truth): y = 2 + 1.5*X1 - 0.8*X2 + 0.5*X3 + eps
y <- 2 + 1.5*X1- 0.8*X2 + 0.5*X3 + eps

# Tambah beberapa outlier pada respon
idx_out <- sample(1:n, 3)
y[idx_out] <- y[idx_out] + rnorm(3, mean = 6, sd = 1)
dat <- data.frame(y, X1, X2, X3)
summary(dat)
##        y                X1                 X2                 X3          
##  Min.   :-3.194   Min.   :-2.30917   Min.   :-2.36565   Min.   :-2.80977  
##  1st Qu.: 1.157   1st Qu.:-0.62576   1st Qu.:-0.68926   1st Qu.:-0.55753  
##  Median : 2.040   Median :-0.05874   Median : 0.03268   Median : 0.07583  
##  Mean   : 2.066   Mean   :-0.00857   Mean   : 0.02408   Mean   : 0.03178  
##  3rd Qu.: 2.895   3rd Qu.: 0.56840   3rd Qu.: 0.69036   3rd Qu.: 0.68098  
##  Max.   : 9.427   Max.   : 3.24104   Max.   : 3.00049   Max.   : 2.43023

2. Estimasi Parameter (OLS)

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

3. Uji Hipotesis

3.3 Uji Hipotesis Gabungan

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

3.4 Uji dengan SE robust (menghadapi heteroskedastisitas)

# Koefisien + SE robust (HC3)
coeftest(mod0, vcov. = sandwich::vcovHC(mod0, type = "HC3"))
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  2.077419   0.087811 23.6578 < 2.2e-16 ***
## X1           1.315357   0.163042  8.0676 7.027e-14 ***
## X2          -0.758334   0.123626 -6.1341 4.644e-09 ***
## X3           0.577370   0.107728  5.3595 2.330e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. Uji Asumsi

4.1 Multikolinearitas

# Variance Inflation Factor (VIF)
car::vif(mod0)
##       X1       X2       X3 
## 1.811923 1.816606 1.003660
# Kondisi Numerik (opsional): kappa > 30 indikasi masalah
kappa(model.matrix(mod0))
## [1] 2.263679

4.2 Heteroskedastisitas

# Breusch-Pagan
lmtest::bptest(mod0)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod0
## BP = 1.0283, df = 3, p-value = 0.7944
# White Test (via bptest dengan kuadrat fitted)
lmtest::bptest(mod0,~ fitted(mod0) + I(fitted(mod0)^2))
## 
##  studentized Breusch-Pagan test
## 
## data:  mod0
## BP = 5.3586, df = 2, p-value = 0.06861

4.3 Normalitas residual

res <- rstandard(mod0)
shapiro.test(res) # sensitif untuk n besar; gunakan juga QQ-plot
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.87436, p-value = 7.729e-12
qqnorm(res); qqline(res)

4.4 Autokorelasi residual (opsional)

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

4.5 Uji spesifikasi (linearitas) – Ramsey RESET

lmtest::resettest(mod0)
## 
##  RESET test
## 
## data:  mod0
## RESET = 4.2735, df1 = 2, df2 = 194, p-value = 0.01527

4.6 Plot Diagnostik Standar

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

par(mfrow = c(1, 1))

5. Seleksi Model

dat2 <- dat %>%
mutate(X1_sq = X1^2, X2_sq = X2^2, X3_sq = X3^2,
X1X2 = X1*X2, X1X3 = X1*X3, X2X3 = X2*X3)

mod_full <- lm(y~ X1 + X2 + X3 + X1_sq + X2_sq + X3_sq + X1X2 + X1X3 + X2X3, data = dat2)
AIC(mod0); AIC(mod_full)
## [1] 657.0574
## [1] 660.7653
# Stepwise berbasis AIC (dua arah)
mod_step <- MASS::stepAIC(mod_full, direction = "both", trace = FALSE)
summary(mod_step)
## 
## Call:
## lm(formula = y ~ X1 + X2 + X3 + X2_sq + X1X2, data = dat2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2286 -0.6853 -0.0298  0.4825  7.0242 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.19439    0.10487  20.925  < 2e-16 ***
## X1           1.31032    0.12513  10.472  < 2e-16 ***
## X2          -0.72799    0.12278  -5.929 1.37e-08 ***
## X3           0.57331    0.08996   6.373 1.32e-09 ***
## X2_sq       -0.28501    0.11895  -2.396   0.0175 *  
## X1X2         0.23671    0.13546   1.748   0.0821 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.22 on 194 degrees of freedom
## Multiple R-squared:  0.4571, Adjusted R-squared:  0.4431 
## F-statistic: 32.67 on 5 and 194 DF,  p-value: < 2.2e-16
broom::glance(mod_step)[,c("r.squared","adj.r.squared","AIC","BIC")]
# Validasi silang (10-fold CV) untuk Membandingkan RMSE:
# Fungsi CV RMSE Menggunakan boot::cv.glm
cv_rmse <- function(model, data, K=10){
  # cv.glm mengembalikan delta: [1] raw CV, [2] adjusted CV
  set.seed(123)
  res <- boot::cv.glm(data = data, glmfit = model, K = K)
  sqrt(res$delta[1])
}
rmse_mod0 <- cv_rmse(mod0, dat)
rmse_modstep<- cv_rmse(mod_step, dat2)
c(RMSE_mod0 = rmse_mod0, RMSE_modStep = rmse_modstep)
##    RMSE_mod0 RMSE_modStep 
##          NaN          NaN

6. Evaluasi Outlier & Titik Berpengaruh

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

# Ambang Sederhana
thr_res <- 3 # |studentized residual| > 3
thr_lev <- 2*length(coef(mod0))/nrow(dat) # leverage tinggi
thr_cook <- 4/nrow(dat) # Cook's distance besar

flag <- data.frame(
  id = 1:nrow(dat),
  studres = studres,
  leverage = lev,
  cookd = cookd
) %>%
mutate(
  flag_res = abs(studres) > thr_res,
  flag_lev = leverage > thr_lev,
  flag_cook = cookd > thr_cook,
  any_flag = flag_res | flag_lev | flag_cook
) %>%
arrange(desc(any_flag), desc(abs(studres)))

head(flag, 10)
car::outlierTest(mod0) # menguji outlier pada residual studentized dengan koreksi Bonferroni
##     rstudent unadjusted p-value Bonferroni p
## 162 6.450812         8.5823e-10   1.7165e-07
## 90  4.302679         2.6674e-05   5.3348e-03
## 18  4.188058         4.2583e-05   8.5167e-03
# Plot Cepat untuk Cook’s distance:
plot(cookd, type="h", main="Cook's Distance", ylab="D", xlab="Observasi")
abline(h = thr_cook, lty=2)

# 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)
broom::tidy(mod_refit)

1.14.2 Soal Latihan

1. Esai: OLS dan Heteroskedastisitas

2.2 Motivasi

2.2.1 Pertumbuhan Bakteri Mengikuti Kurva Logistik

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()

2.2.2 Respon Obat Mengikuti Fungsi Michaelis-Menten

# 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()

2.2.3 Ekonomi: Diminishing Return Mengikuti Bentuk Kurva Non-linear

# 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()

2.4 Estimasi & Inferensi

2.4.1 Estimasi Parameter Menggunakan Metode Iteratif (Gauss-Newton, Levenberg-Marquardt)

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

1.14.2 Soal Latihan

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

Jawaban:

OLS tetap tidak bias meskipun ada heteroskedastisitas karena syarat eksogenitas masih terpenuhi, artinya error tidak berkorelasi dengan variabel X sehingga rata-rata hasil estimasi tetap benar. Namun, OLS menjadi tidak efisien karena ada metode lain seperti GLS atau Weighted LS yang bisa menghasilkan estimasi dengan ragam lebih kecil ketika error memiliki variasi yang tidak seragam. Masalah terbesar heteroskedastisitas sebenarnya bukan pada nilai koefisiennya, tetapi pada standar error yang salah perhitungan sehingga uji-t, uji-F, maupun confidence interval bisa menyesatkan. Di sinilah robust standard errors berperan, yaitu dengan menghitung standar error menggunakan rumus yang konsisten terhadap heteroskedastisitas, sehingga inferensi atau kesimpulan statistik tetap valid meskipun data tidak memenuhi asumsi homoskedastisitas.

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

Jawaban:

Pada model regresi \(y_i = \beta_0 + \beta_1 x_i + u_i\) dengan ragam error \(Var(u_i \mid x_i) = \sigma^2 |x_i|\), metode Weighted Least Squares (WLS) digunakan untuk mengatasi heteroskedastisitas. Bobot yang digunakan adalah \(w_i = 1/|x_i|\), sehingga observasi dengan nilai \(x_i\) besar (ragam error besar) akan memiliki bobot kecil, sedangkan observasi dengan nilai \(x_i\) kecil akan memiliki bobot besar. Estimator WLS dalam bentuk matriks dituliskan sebagai \(\hat{\beta}_{WLS} = (X'WX)^{-1}X'Wy\), dengan \(W\) merupakan matriks diagonal berisi bobot \(w_i\). Estimator ini juga bisa diperoleh dengan cara mentransformasi data, yaitu membagi \(y_i\), \(x_i\), dan konstanta dengan \(\sqrt{|x_i|}\), kemudian menjalankan OLS pada data hasil transformasi tersebut. Untuk kasus regresi sederhana, koefisien slope dapat ditulis dalam bentuk rata-rata berbobot, yaitu \(\hat{\beta}_1 = \frac{\sum w_i (x_i - \bar{x}_w)(y_i - \bar{y}_w)}{\sum w_i (x_i - \bar{x}_w)^2}\) dan \(\hat{\beta}_0 = \bar{y}_w - \hat{\beta}_1 \bar{x}_w\), dengan \(\bar{x}_w = \frac{\sum w_i x_i}{\sum w_i}\) dan \(\bar{y}_w = \frac{\sum w_i y_i}{\sum w_i}\)

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

Jawaban:

Ridge Regression dan Lasso sama-sama memodifikasi regresi OLS yang semula meminimalkan \(\sum_{i=1}^n (y_i - x_i'\beta)^2\). Pada Ridge, fungsi yang diminimalkan ditambahkan penalti kuadrat koefisien sehingga menjadi \(\sum_{i=1}^n (y_i - x_i'\beta)^2 + \lambda \sum_{j=1}^p \beta_j^2\). Penalti ini membuat koefisien menyusut mendekati nol, tetapi jarang benar-benar nol. Akibatnya bias meningkat, namun variance menurun sehingga model lebih stabil, terutama ketika ada multikolinearitas atau jumlah variabel banyak. Sementara itu, Lasso menggunakan penalti absolut, yaitu \(\sum_{i=1}^n (y_i - x_i'\beta)^2 + \lambda \sum_{j=1}^p |\beta_j|\), yang tidak hanya menyusutkan koefisien, tetapi juga dapat memaksa sebagian koefisien menjadi nol. Hal ini menyebabkan bias biasanya lebih besar dibanding Ridge, tetapi variance bisa ditekan lebih jauh. Dari sisi seleksi variabel, Ridge tidak melakukan seleksi karena semua variabel tetap ada dalam model meskipun koefisiennya kecil, sedangkan Lasso mampu melakukan seleksi otomatis dengan meniadakan variabel yang dianggap tidak relevan. Dengan demikian, Ridge lebih cocok digunakan ketika semua variabel dianggap penting untuk dipertahankan, sedangkan Lasso lebih sesuai ketika diinginkan model yang lebih sederhana dengan variabel terpilih

4. Diberikan matriks 𝑋 dan vektor 𝑦, hitunĝ 𝛽 OLS serta 𝐻; identifikasi leverage maksimum

Jawaban:

Pada regresi linear dengan matriks desain \(X\) (ukuran \(n \times p\)) dan vektor respon \(y\) (ukuran \(n \times 1\)), estimator Ordinary Least Squares (OLS) diperoleh dari

\[ \hat{\beta} = (X^\top X)^{-1} X^\top y \]

Matriks proyeksi atau hat matrix didefinisikan sebagai

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

Leverage setiap observasi \(i\) adalah elemen diagonal \(H\) yang ditulis sebagai

\[ h_{ii} = x_i^\top (X^\top X)^{-1} x_i \]

dengan \(x_i\) adalah vektor baris ke-\(i\) dari \(X\). Nilai leverage maksimum ditentukan oleh

\[ \max_i h_{ii} \]

Secara umum, leverage memiliki sifat \(0 \leq h_{ii} \leq 1\) dan \(\sum_{i=1}^n h_{ii} = p\), sehingga rata-rata leverage adalah \(p/n\). Sebuah observasi dianggap memiliki leverage tinggi jika \(h_{ii}\) jauh lebih besar dari rata-rata, misalnya lebih dari \(2p/n\) atau \(3p/n\). Observasi dengan leverage maksimum merupakan pengamatan yang paling berpengaruh terhadap nilai prediksi model

5.A. VIF tinggi menunjukkan?

Jawaban:

Variance Inflation Factor (VIF) digunakan untuk mendeteksi adanya multikolinearitas dalam model regresi. Secara matematis, VIF didefinisikan sebagai

\[ \text{VIF}_j = \frac{1}{1 - R_j^2} \]

dengan \(R_j^2\) adalah koefisien determinasi dari regresi variabel \(X_j\) terhadap semua variabel penjelas lainnya. Nilai VIF yang tinggi menunjukkan bahwa variabel \(X_j\) sangat berkorelasi dengan variabel lain sehingga varians koefisien \(\hat{\beta}_j\) membesar. Secara praktis, VIF mendekati 1 menandakan hampir tidak ada multikolinearitas, sedangkan nilai VIF lebih dari 10 biasanya dianggap indikasi adanya masalah serius dalam model.

5.B. Cook’s D digunakan untuk?

Jawaban:

Cook’s Distance (Cook’s D) digunakan untuk mengukur seberapa besar pengaruh suatu observasi terhadap hasil estimasi koefisien regresi. Rumusnya adalah

\[ D_i = \frac{( \hat{\beta} - \hat{\beta}_{(i)} )^\top (X^\top X)( \hat{\beta} - \hat{\beta}_{(i)} )}{p \, \hat{\sigma}^2} \]

dengan \(\hat{\beta}\) adalah estimator OLS menggunakan semua data, \(\hat{\beta}_{(i)}\) adalah estimator tanpa observasi ke-\(i\), \(p\) adalah jumlah parameter, dan \(\hat{\sigma}^2\) adalah estimator varians error. Nilai Cook’s D yang besar menunjukkan bahwa penghapusan observasi \(i\) akan menyebabkan perubahan signifikan pada koefisien regresi, sehingga observasi tersebut dianggap sebagai influential point. Secara praktis, titik dengan \(D_i > \tfrac{4}{n}\) (dengan \(n\) adalah jumlah observasi) sering diidentifikasi sebagai pengamatan yang perlu diperiksa lebih lanjut.

5.C. Dalam AR(1), kovarians antara 𝜀𝑡 dan 𝜀𝑡−𝑘 adalah?

Jawaban:

Dalam model AR(1) yang dituliskan sebagai

\[ y_t = \phi y_{t-1} + \varepsilon_t \]

dengan \(\varepsilon_t\) berdistribusi white noise yaitu identik dan saling bebas dengan rataan nol dan ragam \(\sigma^2\), kovarians antar error adalah nol untuk jarak waktu berbeda. Artinya, \(\text{Cov}(\varepsilon_t, \varepsilon_{t-k}) = 0\) untuk setiap \(k \neq 0\), dan hanya bernilai \(\sigma^2\) ketika \(k = 0\). Hal ini terjadi karena error dalam model AR(1) diasumsikan tidak berkorelasi antarwaktu. Jika yang dimaksud bukan error melainkan proses \(y_t\), maka autokovariansnya bergantung pada parameter \(\phi\), yaitu

\[ \gamma_k = \phi^k \gamma_0, \quad \text{dengan } \gamma_0 = \frac{\sigma^2}{1-\phi^2} \]