Responsi 9 - MSR

Intermezo

Materi Kuliah

Disini kita akan coba membuktikan pendugaan parameter dari distribusi eksponensial dari dua rumus berikut.

Problem Statement

Diketahui dua penduga untuk mean populasi eksponensial \(\theta\) yaitu:

• \(\theta_1 = \frac{1}{n} \sum_{i=1}^{n} X_i\)
• \(\theta_2 = \frac{1}{n-1} \sum_{i=1}^{n} X_i\)

Apakah \(\theta_1\) dan \(\theta_2\) tak bias dan/atau konsisten? Apakah \(\theta_1\) lebih efisien dari \(\theta_2\)? Tunjukkan dengan simulasi statistika!

Planning: Aims

Tujuan Simulasi:

1. Menunjukkan apakah \(\theta_1\) dan \(\theta_2\) merupakan penduga tak bias bagi mean \(\theta\) untuk populasi sebaran eksponensial.
2. Menunjukkan apakah \(\theta_1\) lebih efisien dibandingkan \(\theta_2\).
3. Menunjukkan apakah \(\theta_1\) dan \(\theta_2\) konsisten.

Planning: Data Generating Mechanism

Data dibangkitkan berdasarkan sebaran Eksponensial dengan rata-rata \(\theta = 10\) dan ukuran contoh \(n = 5, 15, 30, 60, 100, 10000\).

Planning: Estimand

Estimand yang digunakan adalah \(\theta = 10\).

Planning: Methods

1. Data dari populasi Eksponensial(\(\theta = 10\)) dengan \(n = 5, 15, 30, 60, 100, 10000\) dibangkitkan sebanyak \(nsim = 1000\) ulangan.
2. Untuk setiap ulangan, hitung \(\theta_1\) dan \(\theta_2\).
3. Hitung performa estimasi untuk masing-masing ukuran contoh.
4. Buat grafik perbandingan bias dan variansi.

Planning: Performance

Ukuran performa yang digunakan:

• Bias: \(\text{Bias}(\hat{\theta}) = \frac{1}{nsim} \sum_{i=1}^{nsim} (\hat{\theta}_i - \theta)\)
• Relative Bias: \(\text{Relative Bias}(\hat{\theta}) = \frac{\text{Bias}(\hat{\theta})}{\theta}\)
• Variansi Empirik: \(\text{Var}(\hat{\theta}) = \frac{1}{nsim-1} \sum_{i=1}^{nsim} (\hat{\theta}_i - \bar{\hat{\theta}})^2\)

Coding and Execution

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.1
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(tibble)

theta <- 2        # true mean
lambda <- 1/theta # rate parameter

n_obs_list <- c(5, 15, 30, 60, 100, 10000)
n_sim <- 1000

results <- list()

for (n_obs in n_obs_list) {
  W1_vals <- numeric(n_sim)
  W2_vals <- numeric(n_sim)
  
  for (i in 1:n_sim) {
    samp <- rexp(n_obs, rate = lambda)
    W1_vals[i] <- mean(samp)
    W2_vals[i] <- (n_obs / (n_obs - 1)) * mean(samp)
  }
  
  df <- tibble(
    estimator = c(rep("W1", n_sim), rep("W2", n_sim)),
    estimate = c(W1_vals, W2_vals),
    n_obs = n_obs
  )
  
  results[[as.character(n_obs)]] <- df
}

final_data <- bind_rows(results)

summary_table <- final_data %>%
  group_by(estimator, n_obs) %>%
  summarise(
    bias = mean(estimate - theta),
    rel_bias = mean((estimate - theta) / theta),
    emp_var = var(estimate),
    .groups = "drop"
  )

summary_table

## # A tibble: 12 × 5
##    estimator n_obs      bias  rel_bias  emp_var
##    <chr>     <dbl>     <dbl>     <dbl>    <dbl>
##  1 W1            5 -0.0381   -0.0191   0.725   
##  2 W1           15  0.00433   0.00217  0.241   
##  3 W1           30 -0.0167   -0.00834  0.137   
##  4 W1           60 -0.00534  -0.00267  0.0662  
##  5 W1          100 -0.0142   -0.00711  0.0394  
##  6 W1        10000  0.000481  0.000241 0.000402
##  7 W2            5  0.452     0.226    1.13    
##  8 W2           15  0.147     0.0737   0.277   
##  9 W2           30  0.0517    0.0259   0.147   
## 10 W2           60  0.0285    0.0142   0.0684  
## 11 W2          100  0.00585   0.00292  0.0402  
## 12 W2        10000  0.000681  0.000341 0.000402

Simpangan (Deviation)

1. Teori Singkat:

• Definisi Simpangan:

Simpangan (atau deviasi) adalah selisih antara nilai data dengan nilai rata-rata atau nilai pusat dari data tersebut. Simpangan memberikan gambaran tentang seberapa jauh data tersebar dari pusatnya. Jika simpangan kecil, maka data lebih terkonsentrasi di sekitar rata-rata. Sebaliknya, jika simpangan besar, data lebih tersebar jauh dari rata-rata.

\[ \text{Simpangan} = x_i - \bar{x} \] di mana:

• \(x_i\) adalah nilai data ke-i

• \(\bar{x}\) adalah rata-rata dari data tersebut.

• Contoh dalam Kehidupan Nyata:

Misalkan kita mengukur tinggi badan sekelompok siswa. Jika sebagian besar siswa memiliki tinggi badan yang dekat dengan rata-rata, maka simpangan tinggi badan siswa tersebut kecil. Namun, jika ada beberapa siswa yang sangat tinggi atau sangat pendek dibandingkan rata-rata, maka simpangan tinggi badan tersebut akan besar.

2. Penerapan Simpangan di R:

a. Simpangan dari Suatu Data Sebaran:

Kita mulai dengan menghitung simpangan dari data yang dihasilkan secara acak berdasarkan distribusi normal. Di sini, kita membangkitkan data dengan distribusi normal, lalu menghitung simpangan masing-masing nilai terhadap rata-rata data tersebut.

# Membangkitkan data normal
population <- rnorm(1000, mean=50, sd=10)  # Data populasi normal, mean=50, sd=10

# Mengambil sampel dari populasi
sample <- sample(population, 30)  # Sampel acak sebanyak 30 data

# Menghitung simpangan
deviations <- sample - mean(sample)  # Simpangan dihitung sebagai selisih antara data dan rata-rata

# Menampilkan histogram dari distribusi simpangan
hist(deviations, main="Distribusi Simpangan", xlab="Simpangan", col="lightblue", border="black")

Penjelasan Kode: - set.seed(123): Mengatur seed agar hasil yang dihasilkan dapat direproduksi (menghasilkan hasil yang sama setiap kali dijalankan).

• rnorm(1000, mean=50, sd=10): Membuat 1000 data acak dari distribusi normal dengan rata-rata 50 dan deviasi standar 10.

• sample(population, 30): Mengambil sampel sebanyak 30 data dari populasi yang telah dibuat.

• deviations <- sample - mean(sample): Menghitung simpangan dengan cara mengurangi setiap nilai sampel dengan rata-rata sampel.

• hist(deviations, ...): Membuat histogram untuk memvisualisasikan distribusi simpangan.

b. Simpangan pada Data dengan Variasi yang Berbeda:

Sekarang, kita akan membandingkan dua dataset yang memiliki deviasi standar (SD) berbeda. Ini akan menunjukkan bagaimana variasi dalam data mempengaruhi simpangan.

# Membangkitkan dua set data dengan standar deviasi yang berbeda
set.seed(123)
data1 <- rnorm(1000, mean=50, sd=10)  # SD = 10
data2 <- rnorm(1000, mean=50, sd=30)  # SD = 30

# Menghitung simpangan dari masing-masing data
deviations1 <- data1 - mean(data1)
deviations2 <- data2 - mean(data2)

# Membuat histogram dari kedua simpangan
par(mfrow=c(1, 2))  # Membagi layar menjadi dua bagian untuk dua histogram
hist(deviations1, main="Simpangan Data 1 (SD=10)", xlab="Simpangan", col="lightblue", border="black")
hist(deviations2, main="Simpangan Data 2 (SD=30)", xlab="Simpangan", col="lightcoral", border="black")

c. Simpangan pada Data yang Mengandung Outlier:

Sekarang, mari kita lihat bagaimana simpangan dihitung pada data yang mengandung nilai ekstrim (outlier). Outlier dapat memperbesar simpangan, sehingga kita akan melihat pengaruhnya terhadap distribusi simpangan.

# Data tanpa outlier
data_without_outlier <- c(50, 52, 54, 53, 51, 51, 49, 48, 47)  # Data normal tanpa outlier

# Data dengan outlier
data_with_outlier <- c(50, 52, 54, 53, 51, 100, 51, 49, 48, 47)  # Nilai 100 adalah outlier

# Menghitung simpangan untuk data tanpa outlier
deviations_without_outlier <- data_without_outlier - mean(data_without_outlier)

# Menghitung simpangan untuk data dengan outlier
deviations_with_outlier <- data_with_outlier - mean(data_with_outlier)

# Menampilkan simpangan untuk kedua data
deviations_without_outlier

## [1] -0.5555556  1.4444444  3.4444444  2.4444444  0.4444444  0.4444444 -1.5555556
## [8] -2.5555556 -3.5555556

deviations_with_outlier

##  [1] -5.5 -3.5 -1.5 -2.5 -4.5 44.5 -4.5 -6.5 -7.5 -8.5

# Membuat histogram perbandingan antara simpangan dengan dan tanpa outlier
par(mfrow=c(1, 2))  # Membagi layar menjadi dua bagian untuk dua histogram
hist(deviations_without_outlier, main="Simpangan Tanpa Outlier", xlab="Simpangan", col="lightblue", border="black")
hist(deviations_with_outlier, main="Simpangan Dengan Outlier", xlab="Simpangan", col="lightcoral", border="black")

Galat (Error)

1. Teori Singkat:

1. Definisi Formal:

Dalam teori estimasi, error (atau galat) didefinisikan sebagai selisih antara penduga parameter dan parameter yang sebenarnya: \[ \text{Error} = \hat{\theta} - \theta \]

di mana:

• \(\hat{\theta}\) adalah estimator atau penduga parameter (misalnya, \(\bar{x}\), rata-rata sampel)

• \(\theta\) adalah parameter yang sebenarnya (misalnya, \(\mu\), rata-rata populasi)

Error ini menunjukkan seberapa jauh estimasi kita dari nilai yang sebenarnya, dan kita biasanya tidak dapat mengukurnya secara langsung karena parameter yang sebenarnya tidak diketahui.

2. Konteks Penting yang Perlu Dipahami:

Konsep	Contoh	Sifat	Dapat Diobservasi?
Error (Galat)	\(\bar{x} - \mu\)	Teoretis	Tidak (karena \(\mu\) tidak diketahui)
Residual (Sisaan)	\(y_i - \hat{y}_i\)	Empiris	Ya
Bias	\(E[\hat{\theta}] - \theta\)	Ekspektasi matematis	Dapat Diestimasi

Penjelasan:

• Error (Galat) adalah teoretis karena kita tidak bisa mengetahui nilai parameter yang sebenarnya (seperti \(\mu\)) secara pasti. Error ini hanya bisa dihitung jika kita tahu parameter yang sesungguhnya.

• Residual (Sisaan) adalah empiris, dihitung sebagai perbedaan antara nilai yang diamati dan nilai yang diprediksi oleh model. Residual bisa diobservasi karena kita memiliki data yang dapat digunakan untuk menghitungnya.

• Bias adalah perbedaan antara ekspektasi dari penduga dan parameter yang sebenarnya. Bias menggambarkan seberapa besar kesalahan sistematik dalam estimasi kita.

3. Nuansa Kritis:

• Dalam Model Regresi:

  • Error (\(\epsilon\)) = \(y - E(y|x)\) (komponen stokastik yang tidak terobservasi)

  • Residual (\(e\)) = \(y - \hat{y}\) (versi estimasi error berdasarkan model)

• Dalam Sampling Theory:

  • Sampling Error: Variabilitas alami karena proses pengambilan sampel. Ini menggambarkan seberapa banyak estimasi kita bisa bervariasi dari sampel ke sampel.

  • Non-Sampling Error: Kesalahan yang terjadi karena masalah pengukuran, ketidaksesuaian sampel, atau faktor lain yang bukan berasal dari pengambilan sampel.

2. Penerapan Error di R:

a. Model Linier dengan Error Acak:

Pada bagian ini, kita akan membuat model linier sederhana dan menambahkan error acak untuk melihat bagaimana error mempengaruhi data.

# Membuat data linier sederhana
set.seed(123)
x <- 1:100
true_relationship <- 2 * x + 10  # Hubungan linear yang sebenarnya: y = 2x + 10

# Menambahkan error acak
error <- rnorm(100, mean=0, sd=15)  # Error acak dengan SD = 15
y <- true_relationship + error  # Nilai y dengan error acak

# Membuat plot untuk visualisasi
plot(x, y, main="Data dengan Error Acak", xlab="X", ylab="Y", pch=16, col="blue")
abline(a=10, b=2, col="red")  # Menambahkan garis model linier tanpa error

Penjelasan Kode:

• x adalah variabel independen (misalnya waktu atau jumlah unit yang diproduksi).

• true_relationship adalah nilai y yang seharusnya berdasarkan hubungan linier tanpa error.

• error adalah error acak yang ditambahkan ke nilai y, dengan distribusi normal yang memiliki rata-rata 0 dan standar deviasi 15.

• y adalah nilai yang dihasilkan dengan menambahkan error acak ke hubungan linier.

• Plot menampilkan data yang dihasilkan, dengan garis merah menunjukkan hubungan linier yang sebenarnya.

b. Menghitung Error Estimasi (Penduga Parameter):

Sekarang, kita akan menghitung error sebagai selisih antara penduga parameter (rata-rata sampel) dan parameter yang sebenarnya (rata-rata populasi).

# Menetapkan parameter sebenarnya
set.seed(123)
mu <- 100          # Parameter sebenarnya
samp <- rnorm(30, mean = mu, sd = 15)  # Sampel acak dengan rata-rata 100 dan SD 15
x_bar <- mean(samp)  # Rata-rata sampel

# Menghitung error sebagai selisih antara rata-rata sampel dan parameter yang sebenarnya
error <- x_bar - mu
error

## [1] -0.7065563

Penjelasan Kode:

• mu <- 100: Menetapkan nilai parameter yang sebenarnya (rata-rata populasi).

• samp <- rnorm(30, mean = mu, sd = 15): Membuat sampel acak sebanyak 30 data dengan rata-rata 100 dan deviasi standar 15.

• x_bar <- mean(samp): Menghitung rata-rata sampel.

• error <- x_bar - mu: Menghitung error sebagai selisih antara rata-rata sampel (penduga parameter) dan nilai parameter yang sebenarnya.

Bias Estimator

Pengertian Bias

Bias dari sebuah estimator \(( \hat{\theta} )\) terhadap parameter \(( \theta )\) didefinisikan sebagai:

\(\text{Bias}(\hat{\theta}) = E[\hat{\theta}] - \theta\)

Sebuah estimator dikatakan tak bias jika nilai harapannya sama dengan parameter yang diestimasi.

Ilustrasi Eksplorasi Bias pada Estimator Rata-rata Eksponensial

kita akan mengeksplorasi apakah estimator rata-rata dari distribusi eksponensial juga merupakan penduga tak bias. Kita menggunakan distribusi eksponensial dengan parameter laju \(( \lambda = 1 )\), sehingga nilai harapan dari sampel eksponensial adalah \(( 1/\lambda = 1 )\).

Perhitungan Teoritis (Rata-rata Eksponensial)

• Untuk distribusi eksponensial \(( X \sim \text{Exp}(\lambda))\), nilai harapan teoritisnya adalah:

  \([E[X] = \frac{1}{\lambda}]\)

  Dengan \(( \lambda = 1 )\), maka \(( E[X] = 1 )\).

• Estimator rata-rata sampel:

  \([\bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i]\)

  adalah tak bias karena:

  \([E[\bar{X}] = \frac{1}{n} \cdot n \cdot \frac{1}{\lambda} = \frac{1}{\lambda}]\)

Simulasi R

1. Menentukan Distribusi Eksponensial: Di sini kita menyebutkan bahwa kita menggunakan distribusi eksponensial dengan parameter \(\lambda = 1\), sehingga nilai harapan teoritis adalah 1.

2. Membuat Simulasi: Langkah ini menjelaskan kita akan melakukan 1000 simulasi dengan ukuran sampel \(n = 10\).

3. Estimator Rata-rata Sampel: Di sini kita menyatakan rumus untuk menghitung estimator rata-rata sampel.

set.seed(321)
sim_expo <- replicate(1000, mean(rexp(10, rate = 1)))
mean(sim_expo)

## [1] 0.9902134

Interpretasi

Estimator rata-rata dari distribusi eksponensial adalah tak bias karena nilai harapan dari \(( \bar{X} )\) adalah \(( 1/\lambda = 1 )\). Hasil rata-rata dari simulasi seharusnya mendekati 1. Jika berbeda, bisa jadi karena ukuran sampel kecil (n = 10), namun secara teori tetap tak bias.

Residual

Pengertian Residual

Residual adalah selisih antara nilai observasi dengan nilai prediksi dari model:

\([ e_i = y_i - \hat{y}_i ]\)

Residual biasa digunakan untuk mengevaluasi kesesuaian model regresi.

Asumsi Residual dalam Model Linier

1. Rata-rata residual = 0
2. Varians residual konstan (homoskedastisitas)
3. Residual berdistribusi normal
4. Tidak ada autokorelasi antar residual

Simulasi Regresi Sederhana

set.seed(42)
n <- 100
x <- rnorm(n, mean = 10, sd = 2)
y <- 3 + 2 * x + rnorm(n, sd = 2)

model <- lm(y ~ x)
residuals <- resid(model)

# Plot residual terhadap x
plot(x, residuals, main = "Plot Residual", xlab = "x", ylab = "Residual")
abline(h = 0, col = "red")

# Plot QQ-plot residual
qqnorm(residuals)
qqline(residuals, col = "blue")

Penjelasan Kode

• set.seed(42): Mengatur random seed agar hasil yang diperoleh bisa direproduksi.

• n <- 100: Menetapkan jumlah observasi.

• x <- rnorm(n, mean = 10, sd = 2): Membuat prediktor dari distribusi normal.

• y <- 3 + 2 * x + rnorm(n, sd = 2): Membuat respons y dengan menambahkan error acak.

• model <- lm(y ~ x): Membuat model regresi linear sederhana.

• residuals <- resid(model): Mengambil nilai residual dari model.

• plot(x, residuals, …): Mengecek asumsi homoskedastisitas dan linearitas.

• qqnorm(residuals): Mengecek apakah residual berdistribusi normal.

• qqline(residuals, col = “blue”): Garis referensi distribusi normal.

Interpretasi

Plot residual terhadap x menunjukkan sebaran acak di sekitar garis nol, yang mengindikasikan bahwa model cukup baik. QQ-plot menunjukkan residual mengikuti distribusi normal secara umum, walau bisa terdapat sedikit penyimpangan di ekor.

Simulasi Heteroskedastisitas

# Data dengan varians error tergantung x
x <- runif(n, 1, 10)
y <- 5 + 3 * x + rnorm(n, sd = x)  # varians meningkat seiring x

model_het <- lm(y ~ x)
resid_het <- resid(model_het)

# Plot data dan garis regresi
plot(x, y, main = "Scatterplot x vs y",
     xlab = "x", ylab = "y", pch = 16, col = "blue")
abline(model_het, col = "red", lwd = 2)  # Garis regresi
legend("topleft", legend = c("Data", "Regresi Linear"),
       col = c("blue", "red"), pch = c(16, NA), lty = c(NA, 1), lwd = c(NA, 2))

plot(x, resid_het, main = "Residual pada Data Heteroskedastik",
     xlab = "x", ylab = "Residual")
abline(h = 0, col = "red")

Diagnostik Model

Gunakan fungsi plot(model) untuk mengevaluasi asumsi residual secara otomatis:

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

Penjelasan Kode

• set.seed(123): Mengatur random seed agar hasil yang diperoleh dapat direproduksi, sehingga kode menghasilkan output yang konsisten setiap kali dijalankan.

• n <- 100: Menetapkan jumlah sampel atau observasi, di mana n berjumlah 100. Ini menentukan banyaknya data yang digunakan dalam simulasi.

• x <- runif(n, 1, 10): Membuat variabel prediktor x dari distribusi uniform antara 1 dan 10. Distribusi uniform digunakan di sini untuk mensimulasikan variabel bebas yang memiliki rentang nilai yang konsisten.

• y <- 5 + 3 * x + rnorm(n, sd = x): Membuat variabel respons y sebagai hasil dari model linear y = 5 + 3x + e, di mana e adalah error acak yang berdistribusi normal dengan standar deviasi yang bergantung pada nilai x. Ini menciptakan kondisi heteroskedastisitas, di mana varians error meningkat seiring bertambahnya nilai x.

• model_het <- lm(y ~ x): Membuat model regresi linear sederhana antara y dan x. Fungsi lm() digunakan untuk mengestimasi koefisien regresi linear pada data yang telah disimulasikan.

• resid_het <- resid(model_het): Mengambil residual dari model regresi yang telah dibangun. Residual ini adalah selisih antara nilai yang diprediksi oleh model dan nilai observasi yang sebenarnya.

• plot(x, y, main = "Scatterplot x vs y", xlab = "x", ylab = "y", pch = 16, col = "blue"): Membuat plot scatter untuk visualisasi hubungan antara x dan y, di mana titik-titik data diberi warna biru dan bentuk bulat (pch = 16).

• plot(x, resid_het, main = "Residual pada Data Heteroskedastik", xlab = "x", ylab = "Residual"): Membuat plot residual terhadap x untuk memeriksa pola residual. Dalam kasus heteroskedastisitas, residual seharusnya tidak tersebar secara merata, dan ini akan terlihat pada grafik.

• par(mfrow = c(2, 2)): Mengatur tampilan grafik menjadi 2x2, sehingga empat plot diagnostik dapat ditampilkan dalam satu kanvas.

• plot(model_het): Menampilkan empat plot diagnostik untuk memeriksa asumsi regresi linear. Plot ini termasuk residual vs fitted, normal Q-Q, scale-location, dan residuals vs leverage, yang masing-masing membantu memeriksa masalah seperti heteroskedastisitas, normalitas residual, dan pengaruh data ekstrim.

Interpretasi

Dari plot residual, tampak terjadinya heteroskedastisitas, yaitu varians residual meningkat seiring nilai x. Hal ini melanggar asumsi homoskedastisitas dalam regresi linier.