Epidemiologi

1 Contoh

Ilustrasi dengan R: OLS vs Polinomial Tinggi (Bias–Variance secara numerik)

1.0.1 1. Simulasi data (DGP)

Kita buat fungsi sebenarnya:

\[ f(x) = x^2 \]

\[ Y = f(x) + \varepsilon \]

dengan

\[ \varepsilon \sim N(0,\sigma^2) \]

set.seed(123)

f_true <- function(x) x^2
sigma <- 0.8
n <- 40

# grid evaluasi untuk melihat bias/variance di berbagai x
x_grid <- seq(-2, 2, length.out = 120)
f_grid <- f_true(x_grid)

# contoh satu dataset
x <- runif(n, -2, 2)
y <- f_true(x) + rnorm(n, 0, sigma)

plot(x, y, pch=16, main="Satu sampel training: Y = x^2 + noise")
lines(x_grid, f_grid, lwd=2)
legend("topleft", legend=c("data", "f(x)=x^2"), pch=c(16, NA), lty=c(NA,1), bty="n")

1.0.2 Fit dua model: linear vs polinomial derajat tinggi

# model sederhana (bias cenderung lebih besar jika f non-linear)
m_lm <- lm(y ~ x)

# model fleksibel (rawan variance tinggi)
deg <- 15
m_poly <- lm(y ~ poly(x, deg, raw = TRUE))

# plot fit pada satu dataset
plot(x, y, pch=16, main="Fit pada satu dataset")
lines(x_grid, f_grid, lwd=2)
lines(x_grid, predict(m_lm, newdata=data.frame(x=x_grid)), lwd=2, lty=2)
lines(x_grid, predict(m_poly, newdata=data.frame(x=x_grid)), lwd=2, lty=3)
legend("topleft",
       legend=c("f(x)=x^2 (true)", "Linear (OLS)", paste0("Polynomial deg=",deg)),
       lty=c(1,2,3), bty="n")

1.1 Mengestimasi Bias^2 dan Variance lewat resampling (simulasi banyak training set)

B <- 300

pred_lm <- matrix(NA_real_, nrow=B, ncol=length(x_grid))
pred_poly <- matrix(NA_real_, nrow=B, ncol=length(x_grid))

set.seed(2026)
for (b in 1:B) {
  x_b <- runif(n, -2, 2)
  y_b <- f_true(x_b) + rnorm(n, 0, sigma)
  
  fit_lm <- lm(y_b ~ x_b)
  fit_poly <- lm(y_b ~ poly(x_b, deg, raw = TRUE))
  
  pred_lm[b,] <- predict(fit_lm,  newdata=data.frame(x_b=x_grid))
  pred_poly[b,] <- predict(fit_poly, newdata=data.frame(x_b=x_grid))
}

# fungsi ringkasan bias-variance
bv_summary <- function(pred_mat, f_grid, sigma2) {
  mu_hat <- colMeans(pred_mat)
  bias2 <- (mu_hat - f_grid)^2
  varhat <- apply(pred_mat, 2, var)
  msehat <- bias2 + varhat + sigma2
  list(mu_hat=mu_hat, bias2=bias2, varhat=varhat, msehat=msehat)
}

s_lm <- bv_summary(pred_lm, f_grid, sigma^2)
s_poly <- bv_summary(pred_poly, f_grid, sigma^2)

# plot bias^2 dan variance (per x)
par(mfrow=c(2,2), mar=c(4,4,2,1))

plot(x_grid, s_lm$bias2, type="l", lwd=2, main="Bias^2 vs x (Linear)", ylab="Bias^2", xlab="x")
plot(x_grid, s_lm$varhat, type="l", lwd=2, main="Variance vs x (Linear)", ylab="Variance", xlab="x")

plot(x_grid, s_poly$bias2, type="l", lwd=2, main=paste0("Bias^2 vs x (Poly ",deg,")"), ylab="Bias^2", xlab="x")
plot(x_grid, s_poly$varhat, type="l", lwd=2, main=paste0("Variance vs x (Poly ",deg,")"), ylab="Variance", xlab="x")

par(mfrow=c(1,1))

1.2 Ringkasan numerik (rata-rata di seluruh grid)

avg_metrics <- function(s) {
  c(
    Bias2 = mean(s$bias2),
    Variance = mean(s$varhat),
    Irreducible = sigma^2,
    MSE = mean(s$msehat)
  )
}

rbind(
  Linear_OLS = avg_metrics(s_lm),
  Poly_high_degree = avg_metrics(s_poly)
)

##                       Bias2     Variance Irreducible          MSE
## Linear_OLS         1.471226 1.382084e-01        0.64 2.249434e+00
## Poly_high_degree 521.285428 1.861103e+05        0.64 1.866322e+05

library(glmnet)

## Warning: package 'glmnet' was built under R version 4.4.3

## Loading required package: Matrix

## Loaded glmnet 4.1-10

# fungsi: fit ridge pada data (dengan fitur polinomial)
ridge_fit_predict <- function(x_train, y_train, x_grid, deg=15, lambda=1.0) {
  Xtr <- sapply(0:deg, function(k) x_train^k)
  Xg  <- sapply(0:deg, function(k) x_grid^k)
  
  fit <- glmnet(Xtr, y_train, alpha=0, lambda=lambda, standardize=TRUE)
  as.numeric(predict(fit, newx=Xg, s=lambda))
}

# ulangi ridge di banyak training set untuk estimasi bias-variance
lambda <- 3
pred_ridge <- matrix(NA_real_, nrow=B, ncol=length(x_grid))

set.seed(2026)
for (b in 1:B) {
  x_b <- runif(n, -2, 2)
  y_b <- f_true(x_b) + rnorm(n, 0, sigma)
  pred_ridge[b,] <- ridge_fit_predict(x_b, y_b, x_grid, deg=deg, lambda=lambda)
}

s_ridge <- bv_summary(pred_ridge, f_grid, sigma^2)

# bandingkan ringkasan
rbind(
  Linear_OLS = avg_metrics(s_lm),
  Poly_high_degree = avg_metrics(s_poly),
  Ridge_poly_features = avg_metrics(s_ridge)
)

##                           Bias2     Variance Irreducible          MSE
## Linear_OLS            1.4712258 1.382084e-01        0.64 2.249434e+00
## Poly_high_degree    521.2854284 1.861103e+05        0.64 1.866322e+05
## Ridge_poly_features   0.2869474 1.070658e-01        0.64 1.034013e+00

# plot mean prediction (rata-rata atas banyak training set)
plot(x_grid, f_grid, type="l", lwd=2, main="Rata-rata prediksi E[f^(x)] dan f(x)",
     ylab="y", xlab="x")
lines(x_grid, s_lm$mu_hat, lwd=2, lty=2)
lines(x_grid, s_poly$mu_hat, lwd=2, lty=3)
lines(x_grid, s_ridge$mu_hat, lwd=2, lty=4)
legend("topleft",
       legend=c("f(x) true", "E[f^(x)] Linear", paste0("E[f^(x)] Poly ",deg),
                paste0("E[f^(x)] Ridge (lambda=",lambda,")")),
       lty=c(1,2,3,4), bty="n")

2 Latihan 2.11.

2.0.1 Ubah fungsi sebenarnya menjadi f(x)=sin(2πx) dan ulangi analisis. Bagaimana berubah bias–variance?

set.seed(123)

f_true <- function(x) sin(2*pi*x)
sigma <- 0.8
n <- 40

x_grid <- seq(-2, 2, length.out = 120)
f_grid <- f_true(x_grid)

x <- runif(n, -2, 2)
y <- f_true(x) + rnorm(n, 0, sigma)

plot(x, y, pch=16, main="Satu sampel training: Y = sin(2πx) + noise")
lines(x_grid, f_grid, lwd=2)
legend("topleft", legend=c("data", "f(x)=sin(2πx)"), pch=c(16, NA), lty=c(NA,1), bty="n")

m_lm <- lm(y ~ x)

deg <- 15
m_poly <- lm(y ~ poly(x, deg, raw = TRUE))

plot(x, y, pch=16, main="Fit pada satu dataset")
lines(x_grid, f_grid, lwd=2)
lines(x_grid, predict(m_lm, newdata=data.frame(x=x_grid)), lwd=2, lty=2)
lines(x_grid, predict(m_poly, newdata=data.frame(x=x_grid)), lwd=2, lty=3)

legend("topleft",
       legend=c("f(x)=sin(2πx)", "Linear (OLS)", paste0("Polynomial deg=",deg)),
       lty=c(1,2,3), bty="n")

1. Linear

• Bias tinggi karena model linear hanya membuat garis lurus sehingga jauh dari fungsi sebenarnya yang berbentuk sinus (bergelombang)
• Variance kecil karena model sederhana sehingga prediksi relatif stabil antar sampel

2. Polinomial

• Bias kecil, mendekati bentuk sinus
• Variance besar karena model sangat sensitif terhadap noise dan perubahan data, sehingga berpotensi mengalami overfitting..

B <- 300

pred_lm <- matrix(NA_real_, nrow=B, ncol=length(x_grid))
pred_poly <- matrix(NA_real_, nrow=B, ncol=length(x_grid))

set.seed(2026)

for (b in 1:B) {
  
  x_b <- runif(n, -2, 2)
  y_b <- f_true(x_b) + rnorm(n, 0, sigma)
  
  fit_lm <- lm(y_b ~ x_b)
  fit_poly <- lm(y_b ~ poly(x_b, deg, raw = TRUE))
  
  pred_lm[b,] <- predict(fit_lm, newdata=data.frame(x_b=x_grid))
  pred_poly[b,] <- predict(fit_poly, newdata=data.frame(x_b=x_grid))
}

bv_summary <- function(pred_mat, f_grid, sigma2) {
  
  mu_hat <- colMeans(pred_mat)
  bias2 <- (mu_hat - f_grid)^2
  varhat <- apply(pred_mat, 2, var)
  msehat <- bias2 + varhat + sigma2
  
  list(mu_hat=mu_hat, bias2=bias2, varhat=varhat, msehat=msehat)
}

s_lm <- bv_summary(pred_lm, f_grid, sigma^2)
s_poly <- bv_summary(pred_poly, f_grid, sigma^2)

par(mfrow=c(2,2), mar=c(4,4,2,1))

plot(x_grid, s_lm$bias2, type="l", lwd=2,
     main="Bias^2 vs x (Linear)", ylab="Bias^2", xlab="x")

plot(x_grid, s_lm$varhat, type="l", lwd=2,
     main="Variance vs x (Linear)", ylab="Variance", xlab="x")

plot(x_grid, s_poly$bias2, type="l", lwd=2,
     main=paste0("Bias^2 vs x (Poly ",deg,")"), ylab="Bias^2", xlab="x")

plot(x_grid, s_poly$varhat, type="l", lwd=2,
     main=paste0("Variance vs x (Poly ",deg,")"), ylab="Variance", xlab="x")

avg_metrics <- function(s) {
  c(
    Bias2 = mean(s$bias2),
    Variance = mean(s$varhat),
    Irreducible = sigma^2,
    MSE = mean(s$msehat)
  )
}

rbind(
  Linear_OLS = avg_metrics(s_lm),
  Poly_high_degree = avg_metrics(s_poly)
)

##                        Bias2     Variance Irreducible          MSE
## Linear_OLS         0.4781825 5.881990e-02        0.64 1.177002e+00
## Poly_high_degree 578.6735752 1.912306e+05        0.64 1.918099e+05

1. Linear: Bias² vs x (kiri atas) Bias kuadrat berfluktuasi di sepanjang domain. Hal ini terjadi karena fungsi sebenarnya berbentuk gelombang sinus, sedangkan model linear hanya dapat merepresentasikan garis lurus, sehingga model tidak mampu menangkap pola periodik dan mengalami underfitting.

2. Linear: Variance vs x (kanan atas) Variance relatif kecil dan membentuk pola seperti huruf U: lebih kecil di tengah domain dan sedikit meningkat di tepi. Ini menunjukkan bahwa model linear cukup stabil terhadap perubahan data training.

3. Poly 15: Bias² vs x (kiri bawah) Bias kuadrat sangat kecil di sebagian besar domain, terutama di area tengah. Hal ini karena polinomial derajat tinggi cukup fleksibel untuk mendekati bentuk gelombang sinus. Namun di tepi domain bias meningkat sangat tajam, menunjukkan model tidak stabil pada boundary.

4. Poly 15: Variance vs x (kanan bawah) Variance hampir nol di tengah domain tetapi melonjak sangat besar di ujung domain. Ini merupakan tanda overfitting dan instabilitas numerik pada polinomial derajat tinggi, di mana prediksi menjadi sangat sensitif terhadap noise data.

library(glmnet)
ridge_fit_predict <- function(x_train, y_train, x_grid, deg=15, lambda=1.0) {
  
  # membuat fitur polinomial
  Xtr <- sapply(0:deg, function(k) x_train^k)
  Xg  <- sapply(0:deg, function(k) x_grid^k)
  
  # fit ridge regression
  fit <- glmnet(Xtr, y_train, alpha=0, lambda=lambda, standardize=TRUE)
  
  # prediksi pada grid x
  as.numeric(predict(fit, newx=Xg, s=lambda))
}

lambda <- 3

pred_ridge <- matrix(NA_real_, nrow=B, ncol=length(x_grid))

set.seed(2026)

for (b in 1:B) {
  
  x_b <- runif(n, -2, 2)
  y_b <- f_true(x_b) + rnorm(n, 0, sigma)
  
  pred_ridge[b,] <- ridge_fit_predict(
    x_train = x_b,
    y_train = y_b,
    x_grid = x_grid,
    deg = deg,
    lambda = lambda
  )
  
  s_ridge <- bv_summary(pred_ridge, f_grid, sigma^2)
  
  rbind(
  Linear_OLS = avg_metrics(s_lm),
  Poly_high_degree = avg_metrics(s_poly),
  Ridge_poly_features = avg_metrics(s_ridge)
)
}

plot(x_grid, f_grid,
     type="l", lwd=2,
     main="Rata-rata prediksi E[f^(x)] dan fungsi sebenarnya",
     xlab="x", ylab="y")

lines(x_grid, s_lm$mu_hat, lwd=2, lty=2)
lines(x_grid, s_poly$mu_hat, lwd=2, lty=3)
lines(x_grid, s_ridge$mu_hat, lwd=2, lty=4)

legend("topleft",
       legend=c("f(x) true",
                "E[f^(x)] Linear",
                paste0("E[f^(x)] Poly ",deg),
                paste0("E[f^(x)] Ridge (lambda=",lambda,")")),
       lty=c(1,2,3,4),
       bty="n")

1. \(E[\hat f(x)]\) Linear

Garis putus-putus untuk model linear terlihat hampir seperti garis lurus dan jauh dari bentuk gelombang sinus. Hal ini menunjukkan bahwa model linear tidak cukup fleksibel untuk menangkap pola periodik pada fungsi sebenarnya sehingga menghasilkan bias yang tinggi (underfitting). Model ini memiliki variance yang relatif kecil.

2. \(E[\hat f(x)]\)} Polinomial derajat 15

Model polinomial derajat tinggi mampu mengikuti pola fungsi sebenarnya di beberapa bagian domain, terutama di tengah, sehingga bias relatif kecil pada area tersebut. Namun pada bagian tepi domain kurva menjadi tidak stabil dan menyimpang dari fungsi sebenarnya. Hal ini menunjukkan bahwa model memiliki variance besar, sehingga prediksi sangat sensitif terhadap perubahan kecil pada data training dan rentan terhadap overfitting.

3. \(E[\hat f(x)]\)} Ridge, Lambda 3

Ridge masih mampu menangkap sebagian pola non-linear, tetapi bentuknya lebih halus dan stabil dibandingkan polinomial tinggi. Hal ini menunjukkan bahwa regularization (\(\lambda = 3\)) bekerja dengan menaikkan bias sedikit namun secara signifikan menurunkan variance

2.0.2 Ubah ukuran sampel n dari 40 menjadi 200. Apa yang terjadi pada variance model fleksibel?

set.seed(123)

f_true <- function(x) sin(2*pi*x)
sigma <- 0.8
n <- 200

x_grid <- seq(-2, 2, length.out = 120)
f_grid <- f_true(x_grid)

x <- runif(n, -2, 2)
y <- f_true(x) + rnorm(n, 0, sigma)

plot(x, y, pch=16, main="Satu sampel training: Y = sin(2πx) + noise")
lines(x_grid, f_grid, lwd=2)
legend("topleft", legend=c("data", "f(x)=sin(2πx)"), pch=c(16, NA), lty=c(NA,1), bty="n")

m_lm <- lm(y ~ x)

deg <- 15
m_poly <- lm(y ~ poly(x, deg, raw = TRUE))

plot(x, y, pch=16, main="Fit pada satu dataset")
lines(x_grid, f_grid, lwd=2)
lines(x_grid, predict(m_lm, newdata=data.frame(x=x_grid)), lwd=2, lty=2)
lines(x_grid, predict(m_poly, newdata=data.frame(x=x_grid)), lwd=2, lty=3)

legend("topleft",
       legend=c("f(x)=sin(2πx)", "Linear (OLS)", paste0("Polynomial deg=",deg)),
       lty=c(1,2,3), bty="n")

Berdasarkan gambar output, ketika ukuran sampel ditingkatkan menjadi n=200, jumlah titik data menjadi lebih banyak dan tersebar lebih rapat di sepanjang domain x. Pola fungsi sebenarnya menjadi lebih jelas terlihat dibandingkan ketika jumlah sampel lebih sedikit. Dengan informasi data yang lebih banyak, model dapat menangkap pola fungsi dengan lebih stabil sehingga variasi prediksi model cenderung lebih kecil.

B <- 300

pred_lm <- matrix(NA_real_, nrow=B, ncol=length(x_grid))
pred_poly <- matrix(NA_real_, nrow=B, ncol=length(x_grid))

set.seed(2026)

for (b in 1:B) {
  
  x_b <- runif(n, -2, 2)
  y_b <- f_true(x_b) + rnorm(n, 0, sigma)
  
  fit_lm <- lm(y_b ~ x_b)
  fit_poly <- lm(y_b ~ poly(x_b, deg, raw = TRUE))
  
  pred_lm[b,] <- predict(fit_lm, newdata=data.frame(x_b=x_grid))
  pred_poly[b,] <- predict(fit_poly, newdata=data.frame(x_b=x_grid))
}

bv_summary <- function(pred_mat, f_grid, sigma2) {
  
  mu_hat <- colMeans(pred_mat)
  bias2 <- (mu_hat - f_grid)^2
  varhat <- apply(pred_mat, 2, var)
  msehat <- bias2 + varhat + sigma2
  
  list(mu_hat=mu_hat, bias2=bias2, varhat=varhat, msehat=msehat)
}

s_lm <- bv_summary(pred_lm, f_grid, sigma^2)
s_poly <- bv_summary(pred_poly, f_grid, sigma^2)

par(mfrow=c(2,2), mar=c(4,4,2,1))

plot(x_grid, s_lm$bias2, type="l", lwd=2,
     main="Bias^2 vs x (Linear)", ylab="Bias^2", xlab="x")

plot(x_grid, s_lm$varhat, type="l", lwd=2,
     main="Variance vs x (Linear)", ylab="Variance", xlab="x")

plot(x_grid, s_poly$bias2, type="l", lwd=2,
     main=paste0("Bias^2 vs x (Poly ",deg,")"), ylab="Bias^2", xlab="x")

plot(x_grid, s_poly$varhat, type="l", lwd=2,
     main=paste0("Variance vs x (Poly ",deg,")"), ylab="Variance", xlab="x")

avg_metrics <- function(s) {
  c(
    Bias2 = mean(s$bias2),
    Variance = mean(s$varhat),
    Irreducible = sigma^2,
    MSE = mean(s$msehat)
  )
}

rbind(
  Linear_OLS = avg_metrics(s_lm),
  Poly_high_degree = avg_metrics(s_poly)
)

##                        Bias2   Variance Irreducible       MSE
## Linear_OLS       0.477622715 0.01100145        0.64 1.1286242
## Poly_high_degree 0.001518939 0.14745757        0.64 0.7889765

1. Linear: Bias² vs x (kiri atas) Bias kuadrat berfluktuasi di sepanjang domain. Hal ini terjadi karena fungsi sebenarnya berbentuk gelombang sinus, sedangkan model linear hanya dapat merepresentasikan garis lurus, sehingga model tidak mampu menangkap pola periodik dan mengalami underfitting.

2. Linear: Variance vs x (kanan atas) Variance relatif kecil dan membentuk pola seperti huruf U: lebih kecil di tengah domain dan meningkat di tepi. Dengan ukuran sampel yang lebih besar, variance menjadi lebih kecil dibandingkan sebelumnya.

3. Poly 15: Bias² vs x (kiri bawah) Bias kuadrat sangat kecil di hampir seluruh domain karena polinomial derajat tinggi cukup fleksibel untuk mendekati bentuk gelombang sinus. Peningkatan ukuran sampel membuat estimasi model lebih stabil sehingga bias di sebagian besar area rendah.

4. Poly 15: Variance vs x (kanan bawah) Variance sangat kecil di sebagian besar domain dan hanya meningkat di bagian tepi. Dibandingkan sampel kecil, lonjakan variance jauh lebih rendah, menunjukkan bahwa dengan data lebih banyak model fleksibel menjadi lebih stabil meskipun masih sensitif di boundary.

library(glmnet)

ridge_fit_predict <- function(x_train, y_train, x_grid, deg=15, lambda=1.0) {
  
  # membuat fitur polinomial
  Xtr <- sapply(0:deg, function(k) x_train^k)
  Xg  <- sapply(0:deg, function(k) x_grid^k)
  
  # fit ridge regression
  fit <- glmnet(Xtr, y_train, alpha=0, lambda=lambda, standardize=TRUE)
  
  # prediksi
  as.numeric(predict(fit, newx=Xg, s=lambda))
}

lambda <- 3

pred_ridge <- matrix(NA_real_, nrow=B, ncol=length(x_grid))

set.seed(2026)

for (b in 1:B) {
  
  x_b <- runif(n, -2, 2)
  y_b <- f_true(x_b) + rnorm(n, 0, sigma)
  
  pred_ridge[b,] <- ridge_fit_predict(
    x_train = x_b,
    y_train = y_b,
    x_grid = x_grid,
    deg = deg,
    lambda = lambda
  )
}

s_ridge <- bv_summary(pred_ridge, f_grid, sigma^2)

rbind(
  Linear_OLS = avg_metrics(s_lm),
  Poly_high_degree = avg_metrics(s_poly),
  Ridge_poly_features = avg_metrics(s_ridge)
)

##                           Bias2   Variance Irreducible       MSE
## Linear_OLS          0.477622715 0.01100145        0.64 1.1286242
## Poly_high_degree    0.001518939 0.14745757        0.64 0.7889765
## Ridge_poly_features 0.455969500 0.01160909        0.64 1.1075786

plot(x_grid, f_grid,
     type="l", lwd=2,
     main="Rata-rata prediksi E[f^(x)] dan fungsi sebenarnya",
     xlab="x", ylab="y")

lines(x_grid, s_lm$mu_hat, lwd=2, lty=2)
lines(x_grid, s_poly$mu_hat, lwd=2, lty=3)
lines(x_grid, s_ridge$mu_hat, lwd=2, lty=4)

legend("topleft",
       legend=c("f(x) true",
                "E[f^(x)] Linear",
                paste0("E[f^(x)] Poly ",deg),
                paste0("E[f^(x)] Ridge (lambda=",lambda,")")),
       lty=c(1,2,3,4),
       bty="n")

1. \(E[\hat f(x)]\) Linear

Garis putus-putus untuk model linear terlihat hampir seperti garis lurus dan jauh dari bentuk gelombang sinus. Hal ini menunjukkan bahwa meskipun jumlah sampel lebih banyak, model linear tetap tidak cukup fleksibel untuk menangkap pola periodik pada fungsi sebenarnya sehingga bias tetap tinggi (underfitting).

2. \(E[\hat f(x)]\)} Polinomial derajat 15

Model polinomial derajat tinggi mampu mengikuti pola fungsi sinus dengan sangat baik di hampir seluruh domain. Dengan ukuran sampel yang lebih besar, kurva prediksi menjadi lebih stabil dan lebih mendekati fungsi sebenarnya dibandingkan ketika jumlah sampel kecil.

3. \(E[\hat f(x)]\)} Ridge, Lambda 3

Kurva Ridge masih berada di antara model linear dan polinomial tinggi. Model ini tidak sepenuhnya mengikuti gelombang sinus, tetapi prediksinya lebih halus dan stabil tanpa penyimpangan besar di tepi domain, menunjukkan bahwa regularization membantu menjaga keseimbangan antara bias dan variance ketika data lebih banyak.

2.0.3 Coba beberapa nilai lambda pada ridge. Kapan model mulai underfit?

# install.packages("glmnet")
library(glmnet)

ridge_fit_predict <- function(x_train, y_train, x_grid, deg=15, lambda=1.0) {
  
  Xtr <- sapply(0:deg, function(k) x_train^k)
  Xg  <- sapply(0:deg, function(k) x_grid^k)
  
  fit <- glmnet(Xtr, y_train, alpha=0, lambda=lambda, standardize=TRUE)
  
  as.numeric(predict(fit, newx=Xg, s=lambda))
}

lambda_values <- c(0.01, 0.1, 1, 3, 10, 50)
ridge_results <- matrix(NA, nrow=length(lambda_values), ncol=4)

for(i in 1:length(lambda_values)){
  
  lambda <- lambda_values[i]
  
  pred_ridge <- matrix(NA_real_, nrow=B, ncol=length(x_grid))
  
  set.seed(2026)
  
  for (b in 1:B) {
    
    x_b <- runif(n, -2, 2)
    y_b <- f_true(x_b) + rnorm(n, 0, sigma)
    
    pred_ridge[b,] <- ridge_fit_predict(x_b, y_b, x_grid, deg=deg, lambda=lambda)
  }
  
  s_ridge <- bv_summary(pred_ridge, f_grid, sigma^2)
  
  ridge_results[i,] <- avg_metrics(s_ridge)
}

rownames(ridge_results) <- paste0("lambda=", lambda_values)

ridge_results

##                  [,1]        [,2] [,3]     [,4]
## lambda=0.01 0.3710712 0.034679027 0.64 1.045750
## lambda=0.1  0.4229368 0.028566376 0.64 1.091503
## lambda=1    0.4511554 0.016919557 0.64 1.108075
## lambda=3    0.4559695 0.011609088 0.64 1.107579
## lambda=10   0.4659028 0.007555439 0.64 1.113458
## lambda=50   0.4846272 0.005882295 0.64 1.130510

cols <- rainbow(length(lambda_values))

plot(x_grid, f_grid, type="l", lwd=3, col="black",
     main="Ridge untuk beberapa nilai lambda",
     ylab="y", xlab="x")

for(i in 1:length(lambda_values)){
  
  lambda <- lambda_values[i]
  
  pred_ridge <- matrix(NA_real_, nrow=B, ncol=length(x_grid))
  
  for (b in 1:B) {
    
    x_b <- runif(n, -2, 2)
    y_b <- f_true(x_b) + rnorm(n, 0, sigma)
    
    pred_ridge[b,] <- ridge_fit_predict(x_b, y_b, x_grid, deg=deg, lambda=lambda)
  }
  
  s_ridge <- bv_summary(pred_ridge, f_grid, sigma^2)
  
  lines(x_grid, s_ridge$mu_hat, lwd=2, col=cols[i])
}

legend("topright",
       legend=c("True function", paste0("lambda=", lambda_values)),
       col=c("black", cols),
       lty=1, bty="n")

• \(\lambda\) kecil (0.01 – 0.1): Model masih sangat fleksibel sehingga dapat mengikuti pola data dengan baik. Pada kondisi ini bias kecil, tetapi variance relatif lebih besar karena model lebih sensitif terhadap variasi data training.

• \(\lambda\) sedang (sekitar 1 – 3): Regularisasi mulai menstabilkan model. Variance menurun dan bias masih relatif kecil, sehingga menghasilkan keseimbangan bias–variance yang lebih baik.

• \(\lambda\) besar (\(\geq 10\)): Penalti terhadap koefisien menjadi sangat kuat sehingga koefisien menyusut mendekati nol. Model menjadi terlalu sederhana sehingga bias meningkat dan variance sangat kecil, sehingga model mulai underfit.

2.0.4 Coba deg=5 vs deg=15. Jelaskan perubahan kompleksitas.

lambda <- 3

deg5 <- 5
deg15 <- 15

m_poly5 <- lm(y ~ poly(x, deg5, raw = TRUE))
m_poly15 <- lm(y ~ poly(x, deg15, raw = TRUE))

plot(x, y, pch=16, main="Perbandingan Polynomial deg=5 vs deg=15")
lines(x_grid, f_grid, lwd=3)

lines(x_grid, predict(m_poly5, newdata=data.frame(x=x_grid)), lwd=2, lty=2)
lines(x_grid, predict(m_poly15, newdata=data.frame(x=x_grid)), lwd=2, lty=3)

legend("topleft",
       legend=c("True function",
                "Polynomial deg=5",
                "Polynomial deg=15"),
       lty=c(1,2,3),
       bty="n")

Polynomial derajat 5 Bias lebih kecil daripada model linear karena model sudah bisa menangkap sebagian pola sinus. Variance relatif kecil atau sedang, sehingga prediksi masih cukup stabil.

Polynomial derajat 15 Bias sangat kecil karena model sangat fleksibel dan dapat mengikuti bentuk fungsi sebenarnya. Namun variance besar karena model menjadi sangat sensitif terhadap perubahan data dan noise.

pred_poly5 <- matrix(NA_real_, nrow=B, ncol=length(x_grid))
pred_poly15 <- matrix(NA_real_, nrow=B, ncol=length(x_grid))

for (b in 1:B) {
  
  x_b <- runif(n, -2, 2)
  y_b <- f_true(x_b) + rnorm(n, 0, sigma)
  
  fit_poly5 <- lm(y_b ~ poly(x_b, deg5, raw = TRUE))
  fit_poly15 <- lm(y_b ~ poly(x_b, deg15, raw = TRUE))
  
  pred_poly5[b,] <- predict(fit_poly5, newdata=data.frame(x_b=x_grid))
  pred_poly15[b,] <- predict(fit_poly15, newdata=data.frame(x_b=x_grid))
}

s_poly5 <- bv_summary(pred_poly5, f_grid, sigma^2)
s_poly15 <- bv_summary(pred_poly15, f_grid, sigma^2)

par(mfrow=c(2,2), mar=c(4,4,2,1))

plot(x_grid, s_poly5$bias2, type="l", lwd=2,
     main="Bias^2 vs x (Poly 5)", ylab="Bias^2", xlab="x")

plot(x_grid, s_poly5$varhat, type="l", lwd=2,
     main="Variance vs x (Poly 5)", ylab="Variance", xlab="x")

plot(x_grid, s_poly15$bias2, type="l", lwd=2,
     main="Bias^2 vs x (Poly 15)", ylab="Bias^2", xlab="x")

plot(x_grid, s_poly15$varhat, type="l", lwd=2,
     main="Variance vs x (Poly 15)", ylab="Variance", xlab="x")

1. Poly 5: Bias² vs x (kiri atas) Bias kuadrat masih cukup terlihat dan berfluktuasi di sepanjang domain. Polinomial derajat 5 belum cukup fleksibel untuk sepenuhnya mengikuti pola gelombang sinus, sehingga masih terdapat perbedaan antara rata-rata prediksi model dan fungsi sebenarnya.

2. Poly 5: Variance vs x (kanan atas) Variance relatif kecil di sebagian besar domain dan sedikit meningkat di bagian tepi. Ini menunjukkan bahwa model polinomial derajat 5 cukup stabil terhadap perubahan data training karena kompleksitasnya masih moderat.

3. Poly 15: Bias² vs x (kiri bawah) Bias kuadrat sangat kecil di hampir seluruh domain karena polinomial derajat 15 sangat fleksibel sehingga mampu mendekati bentuk fungsi sinus dengan baik. Hal ini menunjukkan bahwa model memiliki kemampuan representasi yang tinggi terhadap pola data.

4. Poly 15: Variance vs x (kanan bawah) Variance sangat besar di bagian tepi domain dan mendekati nol di tengah domain. Hal ini menunjukkan bahwa model polinomial derajat tinggi sangat sensitif terhadap perubahan data training, terutama pada boundary, sehingga berpotensi mengalami overfitting.

rbind(
  Poly_deg5 = avg_metrics(s_poly5),
  Poly_deg15 = avg_metrics(s_poly15)
)

##                   Bias2   Variance Irreducible       MSE
## Poly_deg5  0.4433923667 0.03372826        0.64 1.1171206
## Poly_deg15 0.0006239756 0.14304473        0.64 0.7836687

# install.packages("glmnet")
library(glmnet)

# fungsi fit ridge dengan fitur polinomial
ridge_fit_predict <- function(x_train, y_train, x_grid, deg=15, lambda=1.0) {
  
  Xtr <- sapply(0:deg, function(k) x_train^k)
  Xg  <- sapply(0:deg, function(k) x_grid^k)
  
  fit <- glmnet(Xtr, y_train, alpha=0, lambda=lambda, standardize=TRUE)
  
  as.numeric(predict(fit, newx=Xg, s=lambda))
}

# install.packages("glmnet")
library(glmnet)

# fungsi fit ridge dengan fitur polinomial
ridge_fit_predict <- function(x_train, y_train, x_grid, deg=15, lambda=1.0) {
  
  Xtr <- sapply(0:deg, function(k) x_train^k)
  Xg  <- sapply(0:deg, function(k) x_grid^k)
  
  fit <- glmnet(Xtr, y_train, alpha=0, lambda=lambda, standardize=TRUE)
  
  as.numeric(predict(fit, newx=Xg, s=lambda))
}

s_ridge <- bv_summary(pred_ridge, f_grid, sigma^2)

rbind(
  Linear_OLS = avg_metrics(s_lm),
  Poly_deg5  = avg_metrics(s_poly5),
  Poly_deg15 = avg_metrics(s_poly15),
  Ridge_poly15 = avg_metrics(s_ridge)
)

##                     Bias2    Variance Irreducible       MSE
## Linear_OLS   0.4776227150 0.011001453        0.64 1.1286242
## Poly_deg5    0.4433923667 0.033728257        0.64 1.1171206
## Poly_deg15   0.0006239756 0.143044734        0.64 0.7836687
## Ridge_poly15 0.4846974799 0.005374844        0.64 1.1300723

plot(x_grid, f_grid,
     type="l", lwd=2,
     main="Rata-rata prediksi E[f^(x)] dan fungsi sebenarnya",
     ylab="y", xlab="x")

lines(x_grid, s_lm$mu_hat, lwd=2, lty=2)
lines(x_grid, s_poly5$mu_hat, lwd=2, lty=3)
lines(x_grid, s_poly15$mu_hat, lwd=2, lty=4)
lines(x_grid, s_ridge$mu_hat, lwd=2, lty=5)

legend("topleft",
       legend=c("f(x) true",
                "E[f^(x)] Linear",
                "E[f^(x)] Poly 5",
                "E[f^(x)] Poly 15",
                paste0("E[f^(x)] Ridge (λ=",lambda,")")),
       lty=c(1,2,3,4,5),
       bty="n")

1. Linear Rata-rata prediksi model linear hampir membentuk garis lurus dan jauh dari bentuk fungsi sinus. Hal ini menunjukkan bahwa model terlalu sederhana sehingga tidak mampu menangkap pola non-linear pada data dan menghasilkan bias yang tinggi (underfitting).

2. Polynomial derajat 5 Model polinomial derajat 5 mulai mampu mengikuti sebagian pola gelombang sinus, meskipun masih terdapat perbedaan dengan fungsi sebenarnya. Hal ini menunjukkan bias lebih kecil dibanding model linear dengan kompleksitas yang masih moderat.

3. Polynomial derajat 15 Model polinomial derajat 15 sangat fleksibel sehingga mampu mengikuti bentuk fungsi sebenarnya di beberapa bagian domain. Namun di bagian tepi kurva menjadi tidak stabil dan menyimpang dari fungsi sebenarnya, yang menunjukkan variance tinggi dan potensi overfitting.

4. Ridge (λ = 3) Ridge menghasilkan kurva yang lebih halus dan stabil dibandingkan polinomial derajat tinggi. Regularization membantu mengurangi kompleksitas model sehingga variance menurun, meskipun dengan sedikit peningkatan bias.