Statistical Computing: Midterm

1 Inferensi dalam Regresi Linier

Setelah membaca bagian ini, kita akan dapat:

1. Memahami distribusi perkiraan regresi.
2. Membuat perkiraan interval untuk parameter regresi, respon rata-rata, dan prediksi.
3. Dapat melakukan uji signifikansi regresi

Model regresi linear yang sederhana, biasa di modelkan sebagai berikut, \[Y_i = \beta_0 + \beta_1 x_i + \epsilon_i\]

dimana \(\epsilon_i \sim N(0, \sigma^2)\). Contohnya saja kita akan menggunakan observasi \((x_i,y_i)\) untuk \(i = 1,2,..., n\). Dalam menentukan nilai \(\beta_0\) dan \(\beta_1\) model yang akan diminimalkan adalah, \[f(\beta_0, \beta_1) = \sum_{i = 1}^{n}(y_i - (\beta_0 + \beta_1 x_i))^2\] Lalu, kita masukkan nilai-nilai \(\hat{\beta}_0\) dan \(\hat{\beta}_1\) yang sudah kita temukan, \[\begin{aligned} \hat{\beta}_1 &= \frac{S_{xy}}{S_{xx}} = \frac{\sum_{i = 1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sum_{i = 1}^{n}(x_i - \bar{x})^2}\\ \hat{\beta}_0 &= \bar{y} - \hat{\beta}_1 \bar{x}.\end{aligned}\] Kita juga akan memperkirakan \(\sigma ^2\) menggunakan \(s_e^2\). Dengan begitu, kita menemukan \(s_e^2\) adalah sama dengan kita memperkirakan \(\sigma\), dimana; \[s_e = \text{RSE} = \sqrt{\frac{1}{n - 2}\sum_{i = 1}^n e_i^2}\] biasa kita menyebutnya RSE (Residual Standar Error).

Ketika kita aplikasikan perhitungan ke data “cars”, kita memperoleh hasilnya sebagai berikut;

data(cars)
stop_dist_model = lm(dist ~ speed, data = cars)
summary(stop_dist_model)

## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

Bagian bawah dari hasil coding, akan membahas nilai-nilai dari Estimate, RSE (Residual Standar Error) dan Multiple R-squared.

Selanjutnya, kita akan membahas semua informasi di bawah Koefisien dan juga F-Statistics, dengan coding sebagai berikut;

plot(dist ~ speed, data = cars,
     xlab = "Speed (in Miles Per Hour)",
     ylab = "Stopping Distance (in Feet)",
     main = "Stopping Distance vs Speed",
     pch  = 20,
     cex  = 2,
     col  = "grey")
abline(stop_dist_model, lwd = 5, col = "darkorange")

Untuk memulai, kita akan mencatat bahwa ada ekspresi lain yang setara untuk \(S_{xy}\) yang tidak kita lihat di bagian terakhir; \[S_{xy}= \sum_{i = 1}^{n}(x_i - \bar{x})(y_i - \bar{y}) = \sum_{i = 1}^{n}(x_i - \bar{x}) y_i\] Mungkin ini kesetaraan yang tidak terduga. Namun, ini akan berguna untuk mengilustrasikan konsep di bagian ini.

Perhatikan bahwa, \(\hat{\beta}_1\) adalah contoh statistik jika dihitung dengan data observasi seperti yang tertulis di atas \(\hat{\beta}_0\).

Pada bagian ini akan lebih mudah untuk menggunakan keduanya \(\hat{\beta}_1\) dan \(\hat{\beta}_0\) sebagai variabel acak. Kita harus mengamati masing-masing nilai \(Y_i\). Maka, kita akan menggunakan notasi yang sedikit berbeda, menggantikan huruf kapital \(Y_i\) dengan huruf kecil \(y_i\).

\[\begin{aligned} \hat{\beta}_1 &= \frac{\sum_{i = 1}^{n}(x_i - \bar{x}) Y_i}{\sum_{i = 1}^{n}(x_i - \bar{x})^2} \\ \hat{\beta}_0 &= \bar{Y} - \hat{\beta}_1 \bar{x} \end{aligned}\]

Berhubungan dengan bagian terakhir, kita berpendapat bahwa perkiraan parameter model yang tidak diketahui ini \(\beta_0\) dan \(\beta_1\) bagus karena kami mendapatkannya dengan meminimalkan kesalahan. Sekarang kita akan membahas teorema Gauss-Markov yang membawa gagasan ini lebih jauh, yang menunjukkan bahwa perkiraan ini sebenarnya adalah perkiraan “terbaik”, dari sudut pandang tertentu.

1.1 Teorema Gauss - Markov

Teorema Gauss - Markov memberitahukan kita bahwa ketika kita ingin memperkirakan parameter sederhana linear model regresion \(\beta_0\) dan \(\beta_1\), maka \(\hat{\beta}_0\) dan \(\hat{\beta}_1\) yang kita peroleh adalah perkiraan tidak bias linier terbaik, atau di singkat BLUE (Kondisi aktual untuk teorema Gauss-Markov lebih santai daripada model SLR).

Sekarang kita akan membahas linier, tidak bias dan terbaik yang berkaitan dengan perkiraan ini.

1.1.1 Linear

Dalam pengaturan SLR bahwa \(x_i\) nilai-nilai dianggap jumlah tetap dan diketahui. Kemudian perkiraan linier adalah salah satu yang dapat ditulis sebagai kombinasi linier dari \(Y_i\). Dalam kasus \(\hat{\beta}_1\) kita bisa melihat bahwa; \[\hat{\beta}_1 = \frac{\sum_{i = 1}^{n}(x_i - \bar{x}) Y_i}{\sum_{i = 1}^{n}(x_i - \bar{x})^2} = \sum_{i = 1}^n k_i Y_i = k_1 Y_1 + k_2 Y_2 + \cdots k_n Y_n\] dimana, \(k_i = \displaystyle\frac{(x_i - \bar{x})}{\sum_{i = 1}^{n}(x_i - \bar{x})^2}\) Dengan cara yang sama, kita dapat menunjukkan bahwa \(\hat{\beta}_0\) dapat ditulis sebagai kombinasi linier dari \(Y_i\). Jadi keduanya \(\hat{\beta}_0\) dan \(\hat{\beta}_1\) adalah penduga linier.

1.1.2 Tidak Bias

Sekarang kita tahu perkiraan linier kita, lalu seberapa baik perkiraan kita? Salah satu ukuran “Kebaikan” dari suatu perkiraan adalah biasnya. Secara khusus, pasti kita lebih menyukai perkiraan yang tidak bias, artinya adalah nilai yang diharapkan adalah parameter yang sedang di perkirakan.

Dalam kasus estimasi regresi, kita memiliki; \[\begin{aligned} \text{E}[\hat{\beta}_0] &= \beta_0 \\ \text{E}[\hat{\beta}_1] &= \beta_1\end{aligned}\] Kasus diatas memberi tahu kita bahwa, ketika kondisi model SLR terpenuhi, maka rata-rata perkiraan kita akan benar. Namun, seperti yang kita lihat pada bagian terakhir saat melakukan simulasi dari model SLR, hal itu tidak berarti bahwa setiap perkiraan individu akan benar. Hanya saja, jika kita mengulangi proses tersebut beberapa kali, maka rata-rata perkiraannya akan benar.

1.1.3 Terbaik

Sekarang, jika kita membatasi diri pada linier dan estimasi tidak bias, bagaimana kita dapat mendefinisikan estimasi terbaik? Estimasi dengan varian minimum.

Pertama, sangat mudah untuk membuat perkiraan \(\beta_1\) yang memiliki varian yang sangat rendah, tetapi tidak bias. Contohnya saja, diberikan; \[\hat{\theta}_{BAD} = 5\] Lalu, kita tahu bahwa \(\hat{\theta}_{BAD}\) adalah nilai konstan, \[\text{Var}[\hat{\theta}_{BAD}] = 0\].

Namun sejak, \[\text{E}[\hat{\theta}_{BAD}] = 5\] kita bisa mengatakan \(\hat{\theta}_{BAD}\) adalah penaksir bias kecuali \(\beta_1=5\), yang tidak akan kita ketahui sebelumnya. Alasan tersebut adalah perkiraan yang buruk (kecuali secara kebetulan \(\beta_1=5\)) meskipun variannya bernilai sangat kecil. Ini adalah bagian dari alasan kita membatasi diri pada perkiraan yang tidak bias. Apa gunanya perkiraan, jika memperkirakan kuantitas yang salah?

Jadi sekarang, pertanyaan yang wajar adalah, apa variansnya \(\hat{\beta}_0\) dan \(\hat{\beta}_1\)?

Formula nya adalah sebagai berikut; \[\begin{aligned} \text{Var}[\hat{\beta}_0] &= \sigma^2 \left(\frac{1}{n} + \frac{\bar{x}^2}{S_{xx}}\right) \\ \text{Var}[\hat{\beta}_1] &= \frac{\sigma^2}{S_{xx}}\end{aligned}\]

Persamaan diatas mengukur variabilitas dari estimasi, karena peluang acak selama pengambilan sampel. Apakah ini “yang terbaik”? Apakah varians ini sekecil kemungkinan yang bisa kita dapatkan? Kalian hanya perlu menerima kata-kata kita bahwa mereka memang benar karena menunjukkan bahwa ini benar berada di luar cakupan kursus ini.

1.2 Distribusi Sampling

Kita telah “mendefinisikan ulang” perkiraan untuk \(\hat{\beta}_0\) dan \(\hat{\beta}_1\) sebagai variabel acak, kemudian kita dapat membahas distribusi samplingnya, yaitu distribusi ketika suatu statistik dianggap sebagai variabel acak.

Sejak keduanya \(\hat{\beta}_0\) dan \(\hat{\beta}_1\) adalah kombinasi linier dari \(Y_i\) dan masing-masing nilai \(Y_i\) berdistribusi normal, maka keduanya \(\hat{\beta}_0\) dan \(\hat{\beta}_1\) juga mengikuti berdistribusi normal.

Kemudian, dengan menggabungkan semua hal diatas, kita sampai pada distribusi \(\hat{\beta}_0\) dan \(\hat{\beta}_1\).

Dimana untuk \(\hat{\beta}_1\) dapat kita katakan bahwa; \[\hat{\beta}_1 = \frac{S_{xy}}{S_{xx}}= \frac{\sum_{i = 1}^{n}(x_i - \bar{x}) Y_i}{\sum_{i = 1}^{n}(x_i - \bar{x})^2} \sim N\left( \beta_1, \ \frac{\sigma^2}{\sum_{i = 1}^{n}(x_i - \bar{x})^2} \right)\] Atau singkatnya kita katakan, \[\hat{\beta}_1 \sim N\left(\beta_1, \frac{\sigma^2}{S_{xx}}\right)\]

Dan untuk \(\hat{\beta}_0\), sebagai berikut; \[\hat{\beta}_0 = \bar{Y} - \hat{\beta}_1 \bar{x} \sim N\left( \beta_0, \ \frac{\sigma^2 \sum_{i = 1}^{n}x_i^2}{n \sum_{i = 1}^{n}(x_i - \bar{x})^2} \right)\] Atau singkatnya kita katakan, \[\hat{\beta}_0 \sim N\left( \beta_0, \sigma^2 \left(\frac{1}{n} + \frac{\bar{x}^2}{S_{xx}}\right) \right)\]

Pada titik ini kita telah lalai untuk membuktikan sejumlah hasil. ALih-alih bekerja melalui turunan yang membosankan dari distribusi sampel ini, kita memilih akan membenarkan hasil ini untuk diri kita sendiri menggunakan simulasi.

Catatan untuk kalian para pembaca: Derivasi dan bukti ini dapat ditambahkan ke lampiran di lain waktu. Kalian juga dapat menemukan hasil ini di hampir semua buku teks regresi linier standar. Di UIUC, hasil ini kemungkinan besar akan disajikan di STAT 424 dan STAT 425. Namun, karena kalian tidak akan diminta untuk melakukan penurunan jenis ini dalam kursus ini, hasil tersebut untuk saat ini dihilangkan.

1.2.1 Simulasi Distribusi Sampel

Untuk memverifikasi hasil diatas, maka kita akan mensimulasikan sampel ukuran \(n\) = 100 dari model; \[Y_i = \beta_0 + \beta_1 x_i + \epsilon_i\] dimana \(\epsilon_i \sim N(0, \sigma^2)\).

Dalam hal ini, parameternya adalah; 1. \(\beta_0=3\) 2. \(\beta_1=6\) 3. \(\sigma^2=4\)

Kemudian, berdasarkan data diatas, kita harus menemukan, \[\hat{\beta}_1 \sim N\left(\beta_1, \frac{\sigma^2}{S_{xx}} \right)\] dan \[\hat{\beta}_0 \sim N\left( \beta_0, \sigma^2 \left(\frac{1}{n} + \frac{\bar{x}^2}{S_{xx}}\right) \right)\]

Pertama kita perlu memutuskan apa yang akan menjadi nilai \(x\) kita untuk simulasi ini. Hal tersebut dikarenakan nilai \(x\) dalam SLR juga dianggap kuantitas yang diketahui. Pilihan nilai-nilai \(x\) itu sewenang-wenang. Disini kita juga mengambil bilangan untuk mengacakan, dan menghitung \(S_{xx}\) yang kita perlukan kedepannya.

set.seed(42)
sample_size = 100 # this is n
x = seq(-1, 1, length = sample_size)
Sxx = sum((x - mean(x)) ^ 2)

Kita juga mengambil nilai parameter kita yang telah kita tentukan diatas.

beta_0= 3
beta_1= 6
sigma = 2

Dengan informasi yang ada, kita akan tahu distribusi sampling kita yang seharusnya;

(var_beta_1_hat = sigma ^ 2 / Sxx)

## [1] 0.1176238

(var_beta_0_hat = sigma ^ 2 * (1 / sample_size + mean(x) ^ 2 / Sxx))

## [1] 0.04

\[\hat{\beta}_1 \sim N( 6, 0.1176238)\] dan, \[\hat{\beta}_0 \sim N( 3, 0.04)\] Maksudnya adalah, \[\begin{aligned} \text{E}[\hat{\beta}_1] &= 6 \\ \text{Var}[\hat{\beta}_1] &= 0.1176238 \end{aligned}\] dan \[\begin{aligned} \text{E}[\hat{\beta}_0] &= 3 \\ \text{Var}[\hat{\beta}_0] &= 0.04\end{aligned}\]

Sekarang kita mensimulasikan data dari model ini sebanyak 10.000 kali. Perhatikan bahwa mungkin ini bukan cara yang baik dalam R untuk melakukan simulasi. Namun, kita melakukan simulasi dengan cara seperti ini agar kita mendapatkan kejelasan. Misalnya, kita bisa menggunakan fungsi sim_slr() dari bagian sebelumnya. Kita juga memilih untuk menyimpan variabel di lingkungan global ketimbang membuat bingkai data untuk setiap kumpulan data simulasi baru.

num_samples = 10000
beta_0_hats = rep(0, num_samples)
beta_1_hats = rep(0, num_samples)

for (i in 1:num_samples) {
  eps = rnorm(sample_size, mean = 0, sd = sigma)
  y   = beta_0 + beta_1 * x + eps
  
  sim_model = lm(y ~ x)
  
  beta_0_hats[i] = coef(sim_model)[1]
  beta_1_hats[i] = coef(sim_model)[2]
}

Setiap kita mesimulasikan data, kita akan memperoleh nilai koefisien yang diperkirakan. Variabel \(\hat{\beta}_0\) dan \(\hat{\beta}_1\) sekarang menyimpan 10.000 nilai simulasi \(\hat{\beta}_0\) dan \(\hat{\beta}_1\) masing-masing.

Kita pertama-tama harus memverifikasi distribusi dari \(\hat{\beta}_1\).

mean(beta_1_hats) # empirical mean

## [1] 6.001998

beta_1            # true mean

## [1] 6

var(beta_1_hats)  # empirical variance

## [1] 0.11899

var_beta_1_hat    # true variance

## [1] 0.1176238

Kita melihat bahwa cara empiris dan benar dan varians sangat mirip. Kita juga memverifikasi bahwa distribusi empiris adalah normal. Untuk melakukannya, kami merencanakan histogram dari \(\hat{\beta}_1\), dan menambahkan kurva untuk distribusi sebennarnya dari \(\hat{\beta}_1\). Kita akan menggunakan prob = TRUE untuk meletakkan histogram pada skala yang sama dengan kurva normal.

# note need to use prob = TRUE
hist(beta_1_hats, prob = TRUE, breaks = 20, 
     xlab = expression(hat(beta)[1]), main = "", border = "dodgerblue")
curve(dnorm(x, mean = beta_1, sd = sqrt(var_beta_1_hat)), 
      col = "darkorange", add = TRUE, lwd = 3)

Kemudian, kita mengulangi proses yang sama untuk \(\hat{\beta}_0\)

mean(beta_0_hats) # empirical mean

## [1] 3.001147

beta_0            # true mean

## [1] 3

var(beta_0_hats)  # empirical variance

## [1] 0.04017924

var_beta_0_hat    # true variance

## [1] 0.04

hist(beta_0_hats, prob = TRUE, breaks = 25, 
     xlab = expression(hat(beta)[0]), main = "", border = "dodgerblue")
curve(dnorm(x, mean = beta_0, sd = sqrt(var_beta_0_hat)),
      col = "darkorange", add = TRUE, lwd = 3)

Dalam simulasi ini, kita hanya mensimulasikan sejumlah sampel yang terbatas. Untuk benar-benar memverifikasi hasil distribusi, kita perlu mengamati jumlah sampel yang tidak terbatas. Namun, plot berikut harus menjelaskan bahwa jika kita terus melakukan simulasi, hasil empirisnya akan semakin mendekati apa yang diharapkan.

par(mar = c(5, 5, 1, 1)) # adjusted plot margins, otherwise the "hat" does not display
plot(cumsum(beta_1_hats) / (1:length(beta_1_hats)), type = "l", ylim = c(5.95, 6.05),
     xlab = "Number of Simulations",
     ylab = expression("Empirical Mean of " ~ hat(beta)[1]),
     col  = "dodgerblue")
abline(h = 6, col = "darkorange", lwd = 2)

par(mar = c(5, 5, 1, 1)) # adjusted plot margins, otherwise the "hat" does not display
plot(cumsum(beta_0_hats) / (1:length(beta_0_hats)), type = "l", ylim = c(2.95, 3.05),
     xlab = "Number of Simulations",
     ylab = expression("Empirical Mean of " ~ hat(beta)[0]),
     col  = "dodgerblue")
abline(h = 3, col = "darkorange", lwd = 2)

1.3 Standar Kesalahan (Standar Errors)

Jadi sekarang kita percaya pada dua hasil distribusi, yaitu; \[\begin{aligned} \hat{\beta}_0 &\sim N\left( \beta_0, \sigma^2 \left(\frac{1}{n} + \frac{\bar{x}^2}{S_{xx}}\right) \right) \\ \hat{\beta}_1 &\sim N\left( \beta_1, \frac{\sigma^2}{S_{xx}} \right) \end{aligned}\]

Kemudian dengan menstandarkan hasil ini, kita menemukan, \[\frac{\hat{\beta}_0 - \beta_0}{\text{SD}[\hat{\beta}_0]} \sim N(0, 1)\] dan \[\frac{\hat{\beta}_1 - \beta_1}{\text{SD}[\hat{\beta}_1]} \sim N(0, 1)\] dimana, \[\text{SD}[\hat{\beta}_0] = \sigma\sqrt{\frac{1}{n} + \frac{\bar{x}^2}{S_{xx}}}\] dan \[\text{SD}[\hat{\beta}_1] = \frac{\sigma}{\sqrt{S_{xx}}}\]

Karena kita tidak tahu \(\sigma\) dalam praktiknya, maka kita harus memperkirakannya menggunakan \(s_e\), yang akan kita masukkan ke ekspresi kita yang ada untuk stnadar deviasi perkiraan kita.

Kedua ekspresi baru ini disebut standar kesalahan (standar errors) yang merupakan perkiraan standar deviasi dari distribusi pengambilan sampel.

\[\text{SE}[\hat{\beta}_0] = s_e\sqrt{\frac{1}{n} + \frac{\bar{x}^2}{S_{xx}}}\]

\[\text{SE}[\hat{\beta}_1] = \frac{s_e}{\sqrt{S_{xx}}}\]

Sekarang jika kita membagi dengan standar kesalahan, alih-alih standar deviasi. Kita akan memperoleh hasil berikut yang memungkinkan kita untuk membuat interval kepercayaan dan melakukan pengujian hipotesis.

\[\frac{\hat{\beta}_0 - \beta_0}{\text{SE}[\hat{\beta}_0]} \sim t_{n-2}\]

\[\frac{\hat{\beta}_1 - \beta_1}{\text{SE}[\hat{\beta}_1]} \sim t_{n-2}\]

Untuk melihat ini, pertama kita harus memperhatikan bahwa,

\[\frac{\text{RSS}}{\sigma^2} = \frac{(n-2)s_e^2}{\sigma^2} \sim \chi_{n-2}^2\]

Juga perlu diingat bahwa variabel acak \(T\) didefinisikan sebagai,

\[T = \frac{Z}{\sqrt{\frac{\chi_{d}^2}{d}}}\]

mengikuti ditribusi \(t\) dengan derajat kebebasan \(d\), dimana \(\chi_{d}^2\) adalah variabel acak \(\chi^2\) dengan derajat kebebasan \(d\).

Kita menulis, \[T \sim t_d\] untuk mengatakan bahwa variabel acak \(T\) mengikuti distribusi \(t\) dengan derajat kebebasan \(d\).

Kemudian kita menggunakan trik klasik dengan “kalikan dengan1” dan menyusun ulang beberapa untuk sampai pada

\[\begin{aligned} \frac{\hat{\beta}_1 - \beta_1}{\text{SE}[\hat{\beta}_1]} &= \frac{\hat{\beta}_1 - \beta_1}{s_e / \sqrt{S_{xx}}} \\ &= \frac{\hat{\beta}_1 - \beta_1}{s_e / \sqrt{S_{xx}}} \cdot \frac{\sigma / \sqrt{S_{xx}}}{\sigma / \sqrt{S_{xx}}} \\ &= \frac{\hat{\beta}_1 - \beta_1}{\sigma / \sqrt{S_{xx}}} \cdot \frac{\sigma / \sqrt{S_{xx}}}{s_e / \sqrt{S_{xx}}} \\ &= \frac{\hat{\beta}_1 - \beta_1}{\sigma / \sqrt{S_{xx}}} \bigg/ \sqrt{\frac{s_e^2}{\sigma^2}} \\ &= \frac{\hat{\beta}_1 - \beta_1}{\text{SD}[\hat{\beta}_1]} \bigg/ \sqrt{\frac{\frac{(n - 2)s_e^2}{\sigma^2}}{n - 2}} \sim \frac{Z}{\sqrt{\frac{\chi_{n-2}^2}{n-2}}} \sim t_{n-2} \end{aligned}\]

dimana \(Z \sim N(0,1)\)

Ingatlah bahwa distribusinya \(t\) mirip dengan standar normal, tetapi dengan ekor yang lebih berat. Ketika derajat kebebasan meningkat, maka distribusi \(t\) menjadi lebih dan lebih seperti standar normal. Dibawah ini kita akan memplot distribusi normal yang standar serta dua contoh distribusi \(t\) dengan derajat kebebasan yang berbeda. Perhatikan bagaimana distribusi \(t\) dengan derajat kebebasan yang lebih besar akan lebih mirip dengan kurva standar normal.

# define grid of x values
x = seq(-4, 4, length = 100)

# plot curve for standard normal
plot(x, dnorm(x), type = "l", lty = 1, lwd = 2,
     xlab = "x", ylab = "Density", main = "Normal vs t Distributions")
# add curves for t distributions
lines(x, dt(x, df = 1), lty = 3, lwd = 2, col = "darkorange")
lines(x, dt(x, df = 10), lty = 2, lwd = 2, col = "dodgerblue")

# add legend
legend("topright", title = "Distributions",
       legend = c("t, df = 1", "t, df = 10", "Standard Normal"), 
       lwd = 2, lty = c(3, 2, 1), col = c("darkorange", "dodgerblue", "black"))

1.4 Interval Keyakinan untuk Kemiringan dan Perpotongan (Intercept)

Ingatlah bahwa interval kepercayaan sering kali berbentuk:

\[\text{EST} \pm \text{CRIT} \cdot \text{SE}\] atau \[\text{EST} \pm \text{MARGIN}\]

Dimana: \(\text{EST}\) adalah perkiraan untuk parameter yang diminati, \(\text{SE}\) adalah perkiraan standar kesalahan (standar errors), dan \(\text{MARGIN} = \text{CRIT} \cdot \text{SE}\)

Kemudian, untuk \(\beta_0\) dan \(\beta_1\) kita dapat membuat interval kepercayaan menggunakan;

\[\hat{\beta}_0 \pm t_{\alpha/2, n - 2} \cdot \text{SE}[\hat{\beta}_0] \quad \quad \quad \hat{\beta}_0 \pm t_{\alpha/2, n - 2} \cdot s_e\sqrt{\frac{1}{n}+\frac{\bar{x}^2}{S_{xx}}}\] dan \[\hat{\beta}_1 \pm t_{\alpha/2, n - 2} \cdot \text{SE}[\hat{\beta}_1] \quad \quad \quad \hat{\beta}_1 \pm t_{\alpha/2, n - 2} \cdot \frac{s_e}{\sqrt{S_{xx}}}\]

dimana \(t_{\alpha/2, n - 2}\) adalah nilai kritis sedemikian rupa \(P(t_{n-2} > t_{\alpha/2, n - 2}) = \alpha/2\)

1.5 Pengujian Hipotesis

“We may speak of this hypothesis as the ‘null hypothesis’, and it should be noted that the null hypothesis is never proved or established, but is possibly disproved, in the course of experimentation” - Ronald Aylmer Fisher

Ingat bahwa uji statistik (\(\text{TS}\)) untuk sebuah pengujian, sering berupa;

\[\text{TS} = \frac{\text{EST} - \text{HYP}}{\text{SE}}\]

dimana, \(\text{Est}\) adalah perkiraan untuk parameter yang diminati, \(\text{HYP}\) adalah nilai hipotesis dari parameter, dan \(\text{SE}\) adalah perkiraan standar kesalahan (standar errors).

Jadi, untuk mengujinya; \[H_0: \beta_0 = \beta_{00} \quad \text{vs} \quad H_1: \beta_0 \neq \beta_{00}\]

kita menggunakan uji statistik,

\[t = \frac{\hat{\beta}_0 - \beta_{00}}{\text{SE}[\hat{\beta}_0]} = \frac{\hat{\beta}_0-\beta_{00}}{s_e\sqrt{\frac{1}{n} + \frac{\bar{x}^2}{S_{xx}}}}\]

dimana, dibawah hipotesis nol, dan mengikuti distribusi \(t\) dengan derajat kebebasan \(n-2\). Kita gunakan \(\beta_{00}\) untuk menunjukkan nilai hipotesis \(\beta_0\).

Begitu pula untuk menguji \[H_0: \beta_1 = \beta_{10} \quad \text{vs} \quad H_1: \beta_1 \neq \beta_{10}\]

kita menggunakan uji statistik \[t = \frac{\hat{\beta}_1-\beta_{10}}{\text{SE}[\hat{\beta}_1]} = \frac{\hat{\beta}_1-\beta_{10}}{s_e / \sqrt{S_{xx}}}\]

yang lagi-lagi, dibawah hipotesis nol, mengikuti distribusi \(t\) dengan derajat kebebasan \(n-2\). Kita sekarang menggunakan \(\beta_{10}\) untuk menunjukkan nilai hipotesis \(\beta_1\).

1.6 Contoh cars

Sekarang kita kembali ke contoh data cars dari bagian terakhir untuk mengilustrasikan konsep ini. Kita pertama-tama menyesuaikan model menggunakan lm(), kemudian menggunakan summary() untuk melihat hasil secara lebih rinci.

stop_dist_model = lm(dist ~ speed, data = cars)
summary(stop_dist_model)

## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

1.6.1 Pengujian dalam R

Sekarang kita akan membahas hasil yang ditampilkan yang disebut coefficients. Pertama, ingatlah bahwa kita dapat mengekstrak informasi ini secara langsung.

names(summary(stop_dist_model))

##  [1] "call"          "terms"         "residuals"     "coefficients" 
##  [5] "aliased"       "sigma"         "df"            "r.squared"    
##  [9] "adj.r.squared" "fstatistic"    "cov.unscaled"

summary(stop_dist_model)$coefficients

##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) -17.579095  6.7584402 -2.601058 1.231882e-02
## speed         3.932409  0.4155128  9.463990 1.489836e-12

Fungsi names() memberitahu kita tentang informasi apa saja yang tersedia, dan kemudian kita menggunakan oeprator $ dan coefficients untuk mengekstrak informasi yang kita inginkan. Dua nilai disini harus segera akrab.

\[\hat{\beta}_0 = -17.5790949\] dan \[\hat{\beta}_1 = 3.9324088\]

yang merupakan perkiraan kita untuk parameter model \(\beta_0\) dan \(\beta_1\).

Sekarang mari kita fokus dengan keluaran baris kedua, yang relevan dengan \(\beta_1\).

summary(stop_dist_model)$coefficients[2,]

##     Estimate   Std. Error      t value     Pr(>|t|) 
## 3.932409e+00 4.155128e-01 9.463990e+00 1.489836e-12

Sekali lagi, nilai pertama Estimate adalah \[\hat{\beta}_1 = 3.9324088\]

Nilai kedua Std. Error adalah standar kesalahan (standar errors) dari \(\hat{\beta_1}\), \[\text{SE}[\hat{\beta}_1] = \frac{s_e}{\sqrt{S_{xx}}} = 0.4155128.\]

Nilai ketiga t value adalah nilai uji statistik untuk pengujian \(H_0: \beta_1 = 0\) vs \(H_1: \beta_1 \neq 0\), \[t = \frac{\hat{\beta}_1-0}{\text{SE}[\hat{\beta}_1]} = \frac{\hat{\beta}_1-0}{s_e / \sqrt{S_{xx}}} = 9.46399\]

Terakhir Pr(>|t|) memberi kita nilai p dari pengujian itu.

\[\text{p-value} = 1.4898365\times 10^{-12}\]

Perhatikan disini, kita secara khusus menguji apakah benar atau tidak bahwa \(\beta_1 = 0\).

Baris pertama keluaran melaporkan nilai yang sama, tetapi untuk \(\beta_0\)

summary(stop_dist_model)$coefficients[1,]

##     Estimate   Std. Error      t value     Pr(>|t|) 
## -17.57909489   6.75844017  -2.60105800   0.01231882

Singkatnya, coding berikut menyimpan informasi summary(stop_dist_model)$coefficients dalam variabel baru stop_dist_model_test_info, lalu mengekstrak setiap elemen menjadi variabel baru yang akan menjelaskan informasi yang dikandungnya.

stop_dist_model_test_info = summary(stop_dist_model)$coefficients

beta_0_hat      = stop_dist_model_test_info[1, 1] # Estimate
beta_0_hat_se   = stop_dist_model_test_info[1, 2] # Std. Error
beta_0_hat_t    = stop_dist_model_test_info[1, 3] # t value
beta_0_hat_pval = stop_dist_model_test_info[1, 4] # Pr(>|t|)

beta_1_hat      = stop_dist_model_test_info[2, 1] # Estimate
beta_1_hat_se   = stop_dist_model_test_info[2, 2] # Std. Error
beta_1_hat_t    = stop_dist_model_test_info[2, 3] # t value
beta_1_hat_pval = stop_dist_model_test_info[2, 4] # Pr(>|t|)

Kami kemudian dapat memverifikasi beberapa ekspresi yang setara: statistik uji \(t\) untuk \(\hat{\beta_1}\) dan nilai p dua sisi yang terkait dengan uji statistik tersebut.

(beta_1_hat - 0) / beta_1_hat_se

## [1] 9.46399

beta_1_hat_t

## [1] 9.46399

2 * pt(abs(beta_1_hat_t), df = length(resid(stop_dist_model)) - 2, lower.tail = FALSE)

## [1] 1.489836e-12

beta_1_hat_pval

## [1] 1.489836e-12

1.6.2 Signifikansi Regresi, Uji-t

Kita berhenti sejenak untuk membahas pentingnya uji regresi. Pertama, perhatikan bahwa berdasarkan hasil distribusi di atas, kita dapat menguji \(\beta_0\) dan \(\beta_1\) terhadap nilai tertentu, dan melakukan pengujian satu dan dua sisi.

Namun, satu tes yang sangat spesifik, \[H_0: \beta_1 = 0 \quad \text{vs} \quad H_1: \beta_1 \neq 0\] paling sering digunakan. Mari kita pikirkan pengujian ini dalam kaitannya dengan model regresi linier sederhana, \[Y_i = \beta_0 + \beta_1 x_i + \epsilon_i\]

Jika kita mengasumsikan hipotesis nol benar, maka \(\beta_1=0\) dan kita memiliki modelnya, \[Y_i = \beta_0 + \epsilon_i\]

Dalam model ini, respon tidak bergantung pada prediktor. Jadi kita bisa memikirkan tes ini dengan cara berikut,

1. Dibawah \(H_0\) tidak ada hubungan linier yang signifikan antara \(x\) dan \(y\).
2. Dibawah \(H_1\) ada hubungan linier yang signifikan antara \(x\) dan \(y\).

Untuk contoh cars,

1. Dibawah \(H_0\) tidak ada hubungan linier yang signifikan antara kecepatan dan jarak berhenti.
2. Dibawah \(H_1\) ada hubungan linier yang signifikan antara kecepatan dan jarak berhenti.

Sekali lagi, pengujian itu terlihat pada keluaran dari summary(), \[\text{p-value} = 1.4898365\times 10^{-12}\]

Dengan nilai p yang sangat rendah ini, kita akan menolak hipotesis nol dengan alasan apapun tingkat \(\alpha\), katakanlah misalnya \(\alpha = 0,01\). Jadi kita katakan ada hubungan linier yang signifikan antara kecepatan dan jarak berhenti. Perhatikan bahwa kita menekankan linier.

Dalam plot data simulasi ini, kita bisa melihat hubungan yang jelas antara \(x\) dan \(y\). Namun, ini bukan hubungan linier. Jika kita memasukkan garis ke data ini, itu sangat datar. Tes yang dihasilkan untuk \(H_0: \beta_1 = 0\) vs \(H_1: \beta_1 \neq 0\) memberikan nilai p yang besar, dalam kasus ini bernilai \(0.7564548\). Jadi, kita akan gagal untuk menolak dan mengatakan bahwa tidak ada hubungan linier yang signifikan antara \(x\) dan \(y\). Kita akan melihat nanti bagaimana menyesuaikan kurva ke data ini menggunakan model “linier” \(H_0: \beta_1 = 0\) vs \(H_1: \beta_1 \neq 0\) hanya dapat mendeteksi hubungan yang garis lurus (linier).

1.6.3 Interval Keyakinan Dalam R

Dengan menggunakan R kita dapat dengan mudah mendapatkan interva kepercayaan untk \(\beta_0\) dan \(\beta_1\).

confint(stop_dist_model, level = 0.99)

##                  0.5 %    99.5 %
## (Intercept) -35.706610 0.5484205
## speed         2.817919 5.0468988

Ini secara otomatis akan menghitung interval kepercayaan 99% untuk keduanya \(\beta_0\) dan \(\beta_1\), baris pertama untuk \(\beta_0\), dan baris kedua untuk \(\beta_1\).

cars sebagai contoh saat menginterpretasikan interval ini, kita katakan 99% yakin bahwa untuk peningkatan kecepatan 1 mil per jam, peningkatan rata-rata jarak berhenti adalah antara \(2.8179187\) dan \(5.0468988\) kaki, yang merupakan interval untuk \(\beta_1\).

Perhatikan bahwa interval kepercayaan 99% ini tidak mengandung nilai hipotesis 0. Karena tidak mengandung 0, ini setara dengan menolak pengujian \(H_0: \beta_1 = 0\) vs \(H_1: \beta_1 \neq 0\) di \(\alpha = 0,01\), yang telah kita lihat sebelumnya.

Kalian harus curiga terhadap interval kepercayaan untuk \(\beta_0\), karena mencakup nilai negatif, yang sesuai dengan jarak henti negatif. Secara teknis interpretasinya adalah 99% kita yakin bahwa jarak berhenti rata-rata sebuah mobil yang menempuh jarak 0 mil per jam adalah antara \(-35.7066103\) dan \(0.5484205\) kaki, tetapi kita tidak begitu percaya, karena kita benar-benar yakin bahwa itu akan menjadi non-negatif.

Catatan, kita dapat mengekstrak nilai tertentu dari keluaran ini dengan beberapa cara. Coding ini tidak dijalankan, dan sebagai gantinya, kalian harus memeriksa hubungannya dengan keluaran coding diatas.

confint(stop_dist_model, level = 0.99)[1,]

##       0.5 %      99.5 % 
## -35.7066103   0.5484205

confint(stop_dist_model, level = 0.99)[1, 1]

## [1] -35.70661

confint(stop_dist_model, level = 0.99)[1, 2]

## [1] 0.5484205

confint(stop_dist_model, parm = "(Intercept)", level = 0.99)

##                 0.5 %    99.5 %
## (Intercept) -35.70661 0.5484205

confint(stop_dist_model, level = 0.99)[2,]

##    0.5 %   99.5 % 
## 2.817919 5.046899

confint(stop_dist_model, level = 0.99)[2, 1]

## [1] 2.817919

confint(stop_dist_model, level = 0.99)[2, 2]

## [1] 5.046899

confint(stop_dist_model, parm = "speed", level = 0.99)

##          0.5 %   99.5 %
## speed 2.817919 5.046899

Kita juga dapat memverifikasi kalkulasi yang R lakukan untuk selang \(\beta_1\).

# store estimate
beta_1_hat = coef(stop_dist_model)[2]

# store standard error
beta_1_hat_se = summary(stop_dist_model)$coefficients[2, 2]

# calculate critical value for two-sided 99% CI
crit = qt(0.995, df = length(resid(stop_dist_model)) - 2)

# est - margin, est + margin
c(beta_1_hat - crit * beta_1_hat_se, beta_1_hat + crit * beta_1_hat_se)

##    speed    speed 
## 2.817919 5.046899

1.7 Interval Keyakinan untuk Respon Rata-Rata

Selain interval kepercayaan untuk \(\beta_0\) dan \(\beta_1\), ada dua perkiraan interval umum lainnya yang digunakan dengan regresi. Pertama disebut interval kepercayaan untuk respons rata-rata. Seringkali, kita menginginkan perkiraan interval untuk mean, \(E[Y \mid X = x]\) untuk nilai \(x\) tertentu.

Dalam situasi ini, kita menggunakan \(\hat{y}(x)\) sebagai perkiraan \(E[Y \mid X = x]\) kita. Kita sedikit memodifikasi notasi kita untuk memperjelas bahwa nilai prediksi adalah fungsi dari nilai \(x\). \[\hat{y}(x) = \hat{\beta}_0 + \hat{\beta}_1 x\] Ingat itu, \[\text{E}[Y \mid X = x] = \beta_0 + \beta_1 x\] Jadi, \(\hat{y}(x)\) adalah perkiraan yang bagus karena tidak bias; \[\text{E}[\hat{y}(x)] = \beta_0 + \beta_1 x\] Kita kemudian bisa mendapatkan, \[\text{Var}[\hat{y}(x)] = \sigma^2 \left(\frac{1}{n}+\frac{(x-\bar{x})^2}{S_{xx}}\right)\] Seperti perkiraan lain yang telah kita lihat, \(\hat{y}(x)\) juga mengikuti distribusi normal. Sejak \(\hat{\beta_0}\) dan \(\hat{\beta_1}\) adalah kombinasi linier dari variabel acak normal, begitu juga dengan \(\hat{y}(x)\). \[\hat{y}(x) \sim N \left(\beta_0 + \beta_1 x, \sigma^2 \left(\frac{1}{n}+\frac{(x-\bar{x})^2}{S_{xx}}\right) \right)\] Dan terakhir, kerena kita perlu memperkirakan varians, kita sampai pada standar kesalahan (standar errors) dari perkiraan kita. \[\text{SE}[\hat{y}(x)] = s_e \sqrt{\frac{1}{n}+\frac{(x-\bar{x})^2}{S_{xx}}}\] Kita kemudian, dapat menggunakan ini untuk menemukan interval kepercayaan untuk respons rata-rata, \[\hat{y}(x) \pm t_{\alpha/2, n - 2} \cdot s_e\sqrt{\frac{1}{n}+\frac{(x-\bar{x})^2}{S_{xx}}}\] Untuk menemukan interval kepercayaan untuk respons rata-rata dalam R, maka kita menggunakan fungsi predict(). Kita memberikan fungsi model pas kita serta data baru, lalu disimpan sebagai data frame (hal ini penting dilakukan agar R mengetahui nama variabel prediktor).

Disini, kita menemukan interval kepercayaan untuk jarak berhenti rata-rata saat mobil melaju 5 mil perjam dan saat mobil melaju 21 mil per jam.

new_speeds = data.frame(speed = c(5, 21))
predict(stop_dist_model, newdata = new_speeds, 
        interval = c("confidence"), level = 0.99)

##         fit       lwr      upr
## 1  2.082949 -10.89309 15.05898
## 2 65.001489  56.45836 73.54462

1.8 Prediksi Interval untuk Pengamatan Baru

Terkadang kita menginginkan perkiraan interval untuk pengamatan baru, \(Y\), untuk nilai \(x\) tertentu. Ini sangat mirip dengan interval untuk respons rata-rata, \(\text{E}[Y \mid X = x]\), tetapi berbeda dalam satu hal yang sangat penting.

Tebakan terbaik kita untuk observasi baru masih \(\hat{y}(x)\). Perkiraan rata-rata masih merupakan prediksi terbaik yang bisa kita buat. Perbedaannya terletak pada jumlah variabilitas. Kita tahu bahwa pengamatan akan bervariasi tentang garis regresi yang sebenarnya menurut distribusi \(N(0, \sigma^2)\). Karena ini, kita menambahkan faktor tambahan \(\sigma^2\) ke variabilitas perkiraan kita, untuk menjelaskan variabilitas pengamatan tentang garis regresi (linier regression). \[\begin{aligned} \text{Var}[\hat{y}(x) + \epsilon] &= \text{Var}[\hat{y}(x)] + \text{Var}[\epsilon] \\[2ex] &= \sigma^2 \left(\frac{1}{n}+\frac{(x-\bar{x})^2}{S_{xx}}\right) + \sigma^2 \\[2ex] &= \sigma^2 \left(1 + \frac{1}{n}+\frac{(x-\bar{x})^2}{S_{xx}}\right) \end{aligned}\] \[\hat{y}(x) + \epsilon \sim N \left(\beta_0 + \beta_1 x, \ \sigma^2 \left(1 + \frac{1}{n}+\frac{(x-\bar{x})^2}{S_{xx}}\right) \right)\] \[\text{SE}[\hat{y}(x) + \epsilon] = s_e \sqrt{1 + \frac{1}{n}+\frac{(x-\bar{x})^2}{S_{xx}}}\] Kita kemudian dapat menemukan interval prediksi menggunakan,

\[\hat{y}(x) \pm t_{\alpha/2, n - 2} \cdot s_e\sqrt{1 + \frac{1}{n}+\frac{(x-\bar{x})^2}{S_{xx}}}\]

Untuk menghitung ini, untuk satu set poin dalam R pemberitahuan hanya ada sedikit perubahan dalam sintaks dari menemukan interval kepercayaan untuk respons rata-rata.

predict(stop_dist_model, newdata = new_speeds, 
        interval = c("prediction"), level = 0.99)

##         fit       lwr       upr
## 1  2.082949 -41.16099  45.32689
## 2 65.001489  22.87494 107.12803

Perhatikan juga bahwa kedua interval ini lebih lebar daripada interval kepercayaan yang sesuai untuk repons rata-rata.

1.9 Keyakinan dan Pita Prediksi

Seringkali kita ingin memplot kedua interval keyakinan untuk respons rata-rata dan interval prediksi untuk semua kemungkinan nilai \(x\). Kita menyebutnya keyakinan dan pita prediksi.

speed_grid = seq(min(cars$speed), max(cars$speed), by = 0.01)
dist_ci_band = predict(stop_dist_model, 
                       newdata = data.frame(speed = speed_grid), 
                       interval = "confidence", level = 0.99)
dist_pi_band = predict(stop_dist_model, 
                       newdata = data.frame(speed = speed_grid), 
                       interval = "prediction", level = 0.99) 

plot(dist ~ speed, data = cars,
     xlab = "Speed (in Miles Per Hour)",
     ylab = "Stopping Distance (in Feet)",
     main = "Stopping Distance vs Speed",
     pch  = 20,
     cex  = 2,
     col  = "grey",
     ylim = c(min(dist_pi_band), max(dist_pi_band)))
abline(stop_dist_model, lwd = 5, col = "darkorange")

lines(speed_grid, dist_ci_band[,"lwr"], col = "dodgerblue", lwd = 3, lty = 2)
lines(speed_grid, dist_ci_band[,"upr"], col = "dodgerblue", lwd = 3, lty = 2)
lines(speed_grid, dist_pi_band[,"lwr"], col = "dodgerblue", lwd = 3, lty = 3)
lines(speed_grid, dist_pi_band[,"upr"], col = "dodgerblue", lwd = 3, lty = 3)
points(mean(cars$speed), mean(cars$dist), pch = "+", cex = 3)

Beberapa hal yang perlu diperhatikan:

1. Kita menggunakan argumen ylim untuk meregangkan sumbu-\(y\) dari plot, karena pita lebih panjang daripada titik.

2. Kita menambahkan titik pada titik (\(x,y\)) tersebut.

   • Ini adalah titik yang akan selalu dilewati garis regresi. (Pikirkan alasannya)
   • Ini adalah titi dimana keyakinan dan pita prediksi menjadi yang tersempit. Lihatlah kesalahan standar keduanya untuk memahami mengapa.

3. Pita prediksi (titik-titik biru) kurang melengkung dibandingkan dengan pita kepercayaan (biru putus-putus). Ini adalah hasil dari faktor \(\sigma^2\) tambahan yang ditambahkan ke varians pada nilai \(x\) berapapun.

1.10 Signifikansi Regresi, Uji-F

Dalam kasus regresi linier sederhana, uji signifikansi regresi \(t\) setara dengan uji lain, yaitu \(F\) menguji signifikansi regresi. Kesetaraan ini hanya akan berlaku benar untuk regresi linier sederhana, dan di bagian selanjutnya kita hanya akan menggunakan \(F\) untuk menguji signifikansi regresi.

Ingat dari bagian terakhir dekomposisi varians yang kita lihat sebelumnya menghihitung \(R^2\), \[\sum_{i=1}^{n}(y_i - \bar{y})^2 = \sum_{i=1}^{n}(y_i - \hat{y}_i)^2 + \sum_{i=1}^{n}(\hat{y}_i - \bar{y})^2\]

atau singkatnya, \[\text{SST} = \text{SSE} + \text{SSReg}.\]

Untuk mengembangkan \(F\) test, kita akan mengatur informasi ini dalam tabel ANOVA

ANOVA, atau Analisis Varians akan menjadi konsep yang sering kita bahas dalam bagian ini. Untuk saat ini, kita akan fokus pada hasil tabel, yaitu \(F\) statistik. \[F = \frac{\sum_{i=1}^{n}(\hat{y}_i - \bar{y})^2 / 1}{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2 / (n - 2)} \sim F_{1, n - 2}\]

yang mengikuti distribusi \(F\) dengan derajat kebebasan \(1\) dan \(n-2\) dibawah hipotesis nol. Sebuah distribusi \(F\) adalah distribusi kontinu yang hanya mengambil nilai positif dan memiliki dua parameter, yaitu dua derajat kebebasan.

Ingatlah, pentingnya uji regresi, \(Y\) tidak tergantung pada \(x\) dalam hipotesis nol. \[H_0: \beta_1 = 0 \quad \quad Y_i = \beta_0 + \epsilon_i\] Sedangkan pada hipotesis alternatif \(Y\) mungkin bergantung pada \(x\) \[H_1: \beta_1 \neq 0 \quad \quad Y_i = \beta_0 + \beta_1 x_i + \epsilon_i\] Kita bisa menggunakan \(F\) statistik untuk melakukan tes ini. \[F = \frac{\sum_{i=1}^{n}(\hat{y}_i - \bar{y})^2 / 1}{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2 / (n - 2)}\] Secara khusus, kita akan menolak null saat \(F\) statistiknya benar, yaitu, ketika ada kemungkinan rendah bahwa pengamatan bisa datang dari model nol secara kebetulan. Kita akan menghitung nilai p untuk kita dengan R.

Untuk melakukan \(F\) tes di R kalian dapat melihat baris terakhir dari output yang dalam fungsi summary() disebut F-statistics yang memberikan nilai uji statistik, derajat kebebasan yang relevan, serta \(p-value\) dari pengujian

summary(stop_dist_model)

## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

Selain itu, kalian dapat menggunakan fungsi anova() dalam menampilkan informasidalam tabel ANOVA.

anova(stop_dist_model)

## Analysis of Variance Table
## 
## Response: dist
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## speed      1  21186 21185.5  89.567 1.49e-12 ***
## Residuals 48  11354   236.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ini juga memberikan nilai p (\(p-value\)) untuk pengujian. Kalian harus memperhatikan bahwa nilai p (\(p-value\)) dari distribusi \(t\) adalah sama. Kalian mungkin juga memperhatikan bahwa nilai statistik pengujian untuk uji \(t\), \(9.46399\), dapat dikuadratkan untuk mendapatkan nilai \(F\) statistik, \(89.5671065\).

Perhatikan bahwa ada cara lain yang setara untuk melakukan ini dalam R, dimana kita akan sering sekali untuk membandingkan dua model.

anova(lm(dist ~ 1, data = cars), lm(dist ~ speed, data = cars))

## Analysis of Variance Table
## 
## Model 1: dist ~ 1
## Model 2: dist ~ speed
##   Res.Df   RSS Df Sum of Sq      F   Pr(>F)    
## 1     49 32539                                 
## 2     48 11354  1     21186 89.567 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pernyataan model lm(dist ~ 1, data = cars) menerapkan model \(Y_i = \beta_0 + \epsilon_i\) ke data mobil (data cars). Catat \(\hat{y} = \bar{y}\) itu ketika \(Y_i = \beta_0 + \epsilon_i\)

Pernyataan model lm(dist ~ speed, data = cars) menerapkan model \(Y_i = \beta_0 + \beta_1 x_i + \epsilon_i\). Kita kemudian dapat menganggap penggunaan anova() ini sebagai perbandingan langsung kedua model (Perhatikan bahwa kita mendapatkan nilai p (\(p-value\)) yang sama lagi).