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 awal

GG1 <- ggplot(dat1, aes(wt, mpg)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = FALSE) +
labs(x = "Weight (1000 lbs)", y = "Miles per Gallon (mpg)",
title = "Hubungan mpg vs wt pada mtcars")
GG1

1.9.2 Estimasi OLS & Ringkasan

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

Koefisien dengan CI 95%

tidy(fit1, conf.int = TRUE)

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"
# 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
# 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
# Estimasi Parameter (OLS)
broom::glance(mod0) # ringkas: R^2, Adj R^2, AIC, dll.
# 3.3 Uji hipotesis gabungan (contoh:
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
# 4.1 Multikolinearitas
# 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
# 4.2 Heteroskedastisitas
# 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
# 4.3 Normalitas residual
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)

# 4.6 Plot Diagnostik Standar
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
# 5. Seleksi Model
# 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
# 5. Seleksi Model
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)

Tugas

Chapter 2

Regresi Nonlinier: Model, Estimasi, dan Aplikasi

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 nonlinear.

# 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

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

Plot

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

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

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

Satu langkah iterasi Gauss–Newton (GN) diperoleh dengan menyelesaikan persamaan normal: (𝐽^⊤𝐽)Δ = 𝐽^⊤.r

gn_step <- function(x, y, theta){
r <- resid_vec(x, y, theta)
J <- J_exp(x, theta)

# Normal equations
lhs <- crossprod(J, J) # J^T J
rhs <- crossprod(J, r) # J^T r
Delta <- solve(lhs, rhs) # langkah GN
list(Delta = as.numeric(Delta),
theta_new = theta + as.numeric(Delta),
r = r, J = J, lhs = lhs, rhs = rhs,
SSE = sum(r^2))
}
theta0 <- c(beta0 = 1.0, beta1 = 0.0) # tebakan awal sengaja "kasar"
step1 <- gn_step(x, y, theta0)
step1
## $Delta
## [1] 0.2542352 1.4943019
## 
## $theta_new
##    beta0    beta1 
## 1.254235 1.494302 
## 
## $r
## [1] 0.7758097 1.6076466 3.2677209 3.9474096 5.6919489 8.6494041
## 
## $J
##      d_beta0 d_beta1
## [1,]       1       0
## [2,]       1       1
## [3,]       1       2
## [4,]       1       3
## [5,]       1       4
## [6,]       1       5
## 
## $lhs
##         d_beta0 d_beta1
## d_beta0       6      15
## d_beta1      15      55
## 
## $rhs
##             [,1]
## d_beta0 23.93994
## d_beta1 86.00013
## 
## $SSE
## [1] 136.6569

Iterasi Beberapa Kali Hingga Konvergen

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

Latihan 1.14.2

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

Ordinary Least Squares (OLS) merupakan metode estimasi yang didasarkan pada asumsi Gauss-Markov. Salah satu asumsi penting adalah homoskedastisitas, yaitu varian error yang konstan untuk semua observasi. Jika asumsi ini dilanggar dan terjadi heteroskedastisitas, maka konsekuensinya berbeda untuk sifat unbiasedness dan efficiency dari estimator OLS.

OLS tetap unbiased di bawah heteroskedastisitas. Hal ini karena ketidakbiasan hanya mensyaratkan asumsi eksogenitas, yaitu nilai harapan error bersyarat pada variabel independen sama dengan nol, E(u|X)=0. Dengan kata lain, selama error tidak berkorelasi dengan variabel penjelas (prediktor), koefisien OLS masih mengestimasi nilai parameter populasi secara rata-rata dengan benar. Oleh sebab itu, adanya heteroskedastisitas tidak menggeser nilai tengah distribusi estimator OLS.

Namun, OLS menjadi tidak efisien dalam arti tidak lagi memiliki varians terkecil dibanding estimator linear tak bias lain (BLUE). Dalam konteks homoskedastisitas, teorema Gauss-Markov menjamin OLS sebagai BLUE. Tetapi ketika varians error berbeda-beda antar observasi, OLS tetap linear dan unbiased, hanya saja varians estimator-nya lebih besar dari estimator alternatif (misalnya Generalized Least Squares/GLS). Dengan demikian, efisiensi hilang karena OLS tidak lagi memanfaatkan informasi mengenai struktur varian yang tidak konstan tersebut.

Masalah yang lebih serius muncul dalam inferensi. Di bawah heteroskedastisitas, formula standard error (SE) OLS konvensional tidak valid, karena didasarkan pada asumsi varian error konstan. Jika standard error yang salah digunakan, maka uji t dan uji F menjadi tidak reliabel—biasanya menghasilkan tingkat signifikansi yang menyesatkan (size distortion).

Untuk mengatasi hal ini, digunakan heteroskedasticity-robust standard errors (White’s robust SE). Robust SE menghitung ulang matriks varian kovarian dari estimator OLS dengan memperhitungkan heteroskedastisitas, tanpa perlu mengetahui bentuk spesifiknya. Secara intuitif, robust SE membobotkan kontribusi masing-masing error sesuai dengan besar kecilnya varians yang dimiliki, sehingga estimasi varian parameter menjadi lebih akurat. Dengan demikian, meskipun OLS tetap digunakan untuk estimasi koefisien, hasil inferensi seperti uji hipotesis dan interval kepercayaan tetap valid, karena standar error yang dipakai tidak lagi bias akibat adanya heteroskedastisitas.

Singkatnya, OLS tetap konsisten dan tidak bias meskipun terjadi heteroskedastisitas, tetapi kehilangan efisiensi. Robust SE berperan penting untuk memperbaiki masalah inferensi, sehingga peneliti tetap dapat melakukan pengujian hipotesis secara valid.

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

Berikut merupakan persamaan model regresi linear klasik.
\[ y = X\beta + u \]

Dengan asumsi eksogenitas:
\[ E(u|X) = 0 \]

dan varian residual diketahui proporsional terhadap \(x_i\), yaitu:
\[ Var(u_i | X) = \sigma^2 h_i, \quad h_i = |x_i| \]

Artinya, semakin besar nilai absolut \(x_i\), maka semakin besar pula penyebaran error. Dalam situasi seperti ini OLS masih bisa dipakai, tetapi tidak lagi efisien karena semua observasi dianggap sama penting, padahal sebagian observasi jauh lebih “berisik”. Untuk memperbaikinya, digunakan Weighted Least Squares (WLS) yang memberi bobot pada setiap observasi berdasarkan besar kecilnya varian error.

Estimator OLS umum ketika \(\Omega\) diketahui adalah:
\[ \hat{\beta}_{OLS} = (X'\Omega^{-1}X)^{-1}X'\Omega^{-1}y \]

Dengan \(\Omega\) adalah matriks kovarian error. Karena \(\Omega = \sigma^2 H\), dengan
\[ H = \text{diag}(x_1, \dots, x_n), \]

maka
\[ \Omega^{-1} = \frac{1}{\sigma^2} H^{-1} \]

Faktor \(\sigma^2\) dapat dieliminasi, sehingga bentuk akhirnya sebagai berikut:
\[ \hat{\beta}_{WLS} = (X'H^{-1}X)^{-1}X'H^{-1}y \]

Interpretasi lebih lanjut dapat diperoleh melalui transformasi bobot. Definisikan matriks diagonal:
\[ W = \text{diag}\left(\frac{1}{x_1}, \dots, \frac{1}{x_n}\right) \]

Dengan mengalikan persamaan regresi dengan \(W\), diperoleh model transformasi:
\[ Wy = WX\beta + Wu \]

Di mana vektor error baru memiliki varian konstan \(\sigma^2 I\). Estimasi OLS pada model transformasi tersebut identik dengan hasil WLS pada model awal. Secara implikatif, bobot yang digunakan dalam WLS adalah kebalikan dari besarnya varians residual. Karena dalam kasus ini varians sebanding dengan \(x_i\), maka setiap observasi dibobotkan dengan faktor \(1/x_i\). Estimator WLS tersebut menghasilkan sifat efisiensi yang tidak dimiliki OLS dalam kondisi heteroskedastisitas.

Untuk memahami perannya, kita bisa lihat dari sisi bobot. Jika varian error sebanding dengan \(x_i\), maka bobot yang tepat adalah kebalikannya, yaitu \(1/x_i\). Observasi dengan \(x_i\) besar (error lebih besar) mendapat bobot kecil, sedangkan observasi dengan \(x_i\) kecil (error lebih kecil) mendapat bobot besar. Dengan cara ini, WLS membuat model lebih efisien karena informasi yang lebih “bersih” punya pengaruh lebih kuat dibanding data yang berisik.

Kesimpulan: ketika varian residual proporsional terhadap \(x_i\), estimator WLS diperoleh dengan memberi bobot \(1/x_i\) pada setiap observasi. Estimator ini mengatasi masalah heteroskedastisitas dengan cara menyeimbangkan kontribusi masing-masing data, sehingga hasil estimasi menjadi lebih efisien daripada OLS.

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

Ridge dan Lasso merupakan dua metode regularisasi yang dikembangkan untuk mengatasi masalah multikolinearitas dan overfitting dalam regresi linear. Keduanya menambahkan penalti terhadap besarnya koefisien, namun menggunakan bentuk penalti yang berbeda. Ridge menggunakan penalti kuadrat (L2), sedangkan Lasso menggunakan penalti absolut (L1). Perbedaan ini membawa implikasi penting pada bias–variance trade-off dan kemampuan seleksi variabel.

Dari sudut pandang bias–variance, baik Ridge maupun Lasso sama-sama menambah bias karena koefisien ditarik mendekati nol akibat penalti. Namun penambahan bias ini justru bermanfaat karena dapat menurunkan varians model. Pada data yang memiliki multikolinearitas tinggi atau jumlah variabel yang besar, OLS cenderung menghasilkan varians estimator yang sangat besar sehingga prediksi menjadi tidak stabil. Ridge dan Lasso mengorbankan sedikit ketepatan dengan menambah bias, tetapi sebagai gantinya memberikan penurunan varians yang signifikan. Hasil akhirnya seringkali justru meningkatkan akurasi prediksi pada data uji.

Perbedaan yang menonjol antara Ridge dan Lasso muncul pada seleksi variabel. Ridge hanya mengecilkan koefisien mendekati nol, tetapi tidak pernah membuatnya benar-benar nol. Artinya, Ridge tidak melakukan seleksi variabel secara eksplisit karena semua variabel tetap berada dalam model meskipun pengaruhnya diperkecil. Sebaliknya, Lasso dengan penalti absolut dapat benar-benar memaksa beberapa koefisien menjadi nol. Hal ini membuat Lasso mampu melakukan seleksi variabel secara otomatis, sehingga lebih cocok digunakan ketika terdapat dugaan bahwa hanya sebagian variabel yang relevan.

Dengan demikian, dapat disimpulkan bahwa baik Ridge maupun Lasso sama-sama mengatur trade-off antara bias dan varians agar model menjadi lebih stabil. Ridge lebih efektif ketika semua variabel dianggap relevan tetapi terdapat masalah multikolinearitas, sedangkan Lasso lebih bermanfaat ketika diperlukan model yang lebih sederhana dengan seleksi variabel secara otomatis.

Diberikan matriks \(X\) dan vektor \(y\), hitung OLS serta \(H\); identifikasi leverage maksimum.

Model regresi linear dalam bentuk matriks dituliskan sebagai berikut:
\[ y = X\beta + u \]

Dengan solusi OLS untuk parameter adalah:
\[ \hat{\beta}_{OLS} = (X'X)^{-1}X'y \]

Selanjutnya didefinisikan matriks proyeksi atau hat matrix, yaitu:
\[ H = X(X'X)^{-1}X' \]

Yang memetakan vektor observasi \(y\) ke nilai prediksi:
\[ \hat{y} = Hy \]

Matriks \(H\) memiliki sifat simetris dan idempoten, sehingga semua informasi leverage dapat diperoleh dari diagonal matriks ini.

Nilai leverage untuk observasi ke-\(i\) didefinisikan sebagai elemen diagonal \(h_{ii}\) dari matriks \(H\). Nilai leverage menunjukkan seberapa besar pengaruh observasi tersebut terhadap nilai fitted \(\hat{y}_i\). Observasi dengan nilai leverage tinggi cenderung memiliki pengaruh besar terhadap garis regresi.

Untuk mengidentifikasi leverage maksimum, langkahnya adalah menghitung semua diagonal \(h_{ii}\) dari matriks \(H\), lalu memilih nilai terbesar. Observasi yang memiliki nilai leverage maksimum inilah yang dianggap paling berpengaruh dalam menentukan hasil regresi.

Dengan demikian, penyelesaian soal ini melibatkan tiga komponen utama, yaitu:
1. Menghitung \(\hat{\beta}_{OLS}\) dengan formula \((X'X)^{-1}X'y\),
2. Membentuk matriks proyeksi \(H = X(X'X)^{-1}X'\),
3. Menentukan leverage maksimum dari elemen diagonal \(H\).

Pilihan Ganda:

  1. VIF tinggi menunjukkan? Variance Inflation Factor (VIF) yang tinggi menunjukkan adanya masalah multikolinearitas antar variabel independen. Jika VIF melebihi nilai ambang tertentu (misalnya 10), maka variabel tersebut memiliki korelasi linear yang kuat dengan variabel lain dalam model, sehingga varians estimasi koefisien regresi membesar dan membuat hasil estimasi kurang stabil.

  2. Cook’s D digunakan untuk? Cook’s Distance digunakan untuk mendeteksi observasi yang berpengaruh besar (influential points) terhadap hasil regresi. Observasi dengan nilai Cook’s D yang relatif tinggi mengindikasikan bahwa jika data tersebut dihilangkan, maka estimasi koefisien regresi akan berubah secara signifikan. Dengan demikian, Cook’s D menjadi alat diagnostik penting untuk menilai kestabilan model terhadap data outlier atau data yang berpengaruh kuat.

  3. Dalam AR(1), kovarians antara \(y_t\) dan \(y_{t-k}\) adalah? Untuk proses AR(1): \[ y_t = \phi y_{t-1} + u_t, \quad u_t \sim iid(0,\sigma^2) \]

Kovarians antara \(y_t\) dan \(y_{t-k}\): \[ Cov(y_t , y_{t-k}) = \phi^k \frac{\sigma^2}{1 - \phi^2}, \quad k \geq 0 \]

Artinya, kovarians menurun secara eksponensial seiring bertambahnya lag \(k\), dengan laju penurunan ditentukan oleh parameter \(\phi\).