Setelah membaca bagian ini, anda akan dapat:
Memahami distribusi dari estimasi regresi.
Membuat estimasi interval untuk parameter regresi, rata-rata respons, dan prediksi.
Melakukan uji signifikansi regresi.
Di bab terakhir, kami mendefinisikan model regresi linear sederhana,
\[Y_i = β_0 + β_1 x_i + ϵ_i\]
dimana \(ϵ_i∼N(0, σ^2)\). Lalu, kita memakai observasi \((x_i, y_i)\) untuk \(i = 1,2,...n\) untuk menemukan nilai dari \(β_0\) dan \(β_1\) yang meminimalkan \[ f(β_0, β_1) = ∑^n _{i=1} (y_i - (β_0 + β_1x_i))^2 \]
Kami menyebut nilai-nilainya \({\hat β_0}\) dan \({\hat β_1}\), dimana kita temukan,
\[ {\hat β_1} = \frac {S_{xy}} {S_{xx}} = \frac {∑^n _{i=1} (x_i-\bar x) ( y_i - \bar y)} {∑^n _{i=1} (x_i-\bar x)^2}\] \[ {\hat β_0} = \bar y - {\hat β_1} \bar x \]
Kami juga mengestimasi \(σ^2\) menggunakan \(s^2_e\). Dengan kata lain kami menemukan bahwa \(s_e\) adalah perkiraaan dari \(σ\)
\[ s_e = RSE = \sqrt {\frac {1} {n-2} ∑^n _{i=1}e^2_i }\]
dimana kamu menyebutnya sebagai RSE yaitu " Residual Standard Error " atau " Kesalahan Standar Relatif"
Saat kita mengaplikasikannya ke data cars, kami memperoleh hasil sebagai berikut:
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
Output dari summary() memberikan informasi yang banyak dan berguna. Bagian call meruju pada rumus yang digunakan. Sedangkan bagian Residual atau Sisa yang menunjukkan perbedaan antara nilai aktual dan nilai yang dipredisi dengan menggunakan rumus linear. Dalam persamaan linear, diharapkan residu sekecil mungkin dan simetris.
Bagian Coefficients menunjukkan koefisien aktual. Jadi dari hasil yang ditunjukkan oleh koefisien, kita dapat memembuat model linear seperti \(JarakBerhenti\) \(=3.9324*\) \(kecepatan-17.5791\). Std.Error adalah kesalahan standar deviasi dari perkiraaan koefisiensi dimana semakin kecil kesalahan standar deviasinya, semakin akuran estimasi atau perkirannya. Bagian terakhir, Pr(>|t|) merupakan sebuah parameter, yaitu, probabilitas bahwa kita kebetulan memiliki nilai intersep yang bukan nol murni secara kebetulan, tanpa korelasi apapun antara variabel. Probabilitas diharapkan sangat rendah sehingga parameter bisa lebih dipercaya.
Dari hasil di atas, koefisien terlihat sangat relevan karena probabilitasnya hampir nol, tetapi intersep kurang relevan karena ada di sekitar 0.01.
Di bab terakhir, kami hanya berdiskusi mengenai nilai dari Estimate (Estimasi), Residual standard Error (Kesalahan Standar Relatif), dan Multiple R-squared(R Majemuk). Di bab ini, kami akan membahas semua informasi mengenai Coefficients (Koefisien) serta F-statistic (F-statistik).
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")
Gambar di atas menunjukkan adanya garis regresi dan titik error.
Untuk memulai,kami memerhatikan bahwa ada persamaan lain yang setara untuk \(S_{xy}\) yang tidak kami lihat di bab terakhir
\[S_{xy} = ∑^n _{i=1} (x_i-\bar x) ( y_i - \bar y) = ∑^n _{i=1} (x_i-\bar x) y_i \]
Ini mungkin merupakan penyetaraan yang mengejutkan (Memungkinkan untuk dibuktikan). Namun, ini akan berguna untuk mengilustrasikan konsep di bab ini.
Perhatikan bahwa, \({\hat β_1}\), adalah sampel statistik jika dihitung dengan data observasi seperti yang tertulis di atas, sebagaimana adanya \({\hat β_0}\).
Namun, di bab ini, sering kali akan lebih mudah untuk menggunakan kedua \({\hat β_1}\) dan \({\hat β_0}\) sebagai variabel acak, yakni, kami belum mengamati nilai masing-masing dari \(Y_i\). Jika demikian, kami akan menggunakan notasi yang sedikit berbeda, menggantikan huruf kapital \(Y_i\) dengan huruf kecil \(y_i\).
\[ {\hat β_1} = \frac {∑^n _{i=1} (x_i- \bar x) Y_i } {∑^n _{i=1} (x_i-\bar x)^2}\]
\[ {\hat β_0} = \bar Y - {\hat β_1} \bar x \]
Di bab terakhir, kami berpendapat bahwa estimasi parameter model yang tidak diketahui yaitu \(β_0\) dan \(β_1\) adalah bagus karena kami mendapatkannya dengan meminimalkan kesalahan. Sekarang kita akan membahas teorema Gauss-Markov yang membawa gagasan ini lebih jauh, menunjukkan bahwa perkiraan atau estimasi ini sebenarnya adalah perkiraan “terbaik”, dari sudut pandang tertentu.
Teorema Gauss-Markov memberitahu kita bahwa ketika mengestimasi parameter model regresi linier sederhana \(β_0\) dan \(β_1\), \({\hat β_0}\) dan \({\hat β_1}\) yang kami peroleh adalah perkiraan tidak bias linier terbaik, atau disingkat BLUE. (Kondisi aktual untuk teorema Gauss-Markov lebih santai daripada model SLR.)
Sekarang, kami akan membahas linear, unbiased (tidak bias), dan best (terbaik) yang berkaitan dengan estimasi ini
Linear
Mengingat kembali, dalam pengaturan SLR bahwa nilai dari \(x_i\) dianggap sebagai besaran tetap dan diketahui. Kemudian, perkitaan linear atau estimasi linear adalah salah satu yang dapat ditulis sebagai kombinasi linear dari \(Y_i\). Dalam kasus \({\hat β_1}\) kami melihat
\[ \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, kamu dapat menunjukkan bahwa \(\hat{\beta}_0\) bisa ditulis sebagai kombinasi linear dari \(Y_i\). Jadi, keduanya \(\hat{\beta}_0\) dan \(\hat{\beta}_1\) adalah estimasi linear
Unbiased atau tidak linear
Sekarang kita tahu bahwa estimasi kita adalah linier, seberapa bagus estimasi-estimasi ini? Salah satu ukuran “goodness” dari suatu perkiraan atau estimasi adalah biasnya. Secara khusus, kami lebih menyukai perkiraan yang tidak bias, yang berarti nilai yang diharapkan adalah parameter yang sedang diperkirakan.
Dalam kasus estimasi regresi, kami memiliki,
\[ \begin{aligned} \text{E}[\hat{\beta}_0] &= \beta_0 \\ \text{E}[\hat{\beta}_1] &= \beta_1. \end{aligned} \]
Ini memberi tahu kita bahwa, ketika kondisi model SLR terpenuhi, rata-rata perkiraan kami akan benar. Namun, seperti yang kita lihat pada bab terakhir saat melakukan simulasi dari model SLR, itu tidak berarti bahwa setiap perkiraan individu akan benar. Hanya saja, jika kita mengulangi proses tersebut berkali-kali, rata-rata perkiraannya akan benar.
Best atau Terbaik
Sekarang, jika kita membatasi pada estimasi linier dan tidak bias, bagaimana kita mendefinisikan estimasi terbaik? Estimasi dengan varians minimum.
Pertama, sangat mudah untuk membuat perkiraan \(β_1\) yang memiliki varian yang sangat rendah, tetapi tidak bias. Misalnya, tentukan:
\[ \hat{\theta}_{BAD} = 5. \] Lalu, karene \(\hat{\theta}_{BAD}\) adalah nilai konstan, \[\text{Var}[\hat{\theta}_{BAD}] = 0.\]
Namun sejak, \[\text{E}[\hat{\theta}_{BAD}] = 5\]
Kami mengatakan \(\hat{θ}_{BAD}\) adalah penaksir bias kecuali \(\beta_1 = 5\), yang tidak akan kami ketahui sebelumnya. Untuk alasan ini, ini adalah perkiraan yang buruk (kecuali secara kebetulan \(\beta_1 = 5\) meskipun variannya sekecil mungkin. Ini adalah bagian dari alasan kami membatasi pada perkiraan yang tidak bias. Apa gunanya perkiraan, jika memperkirakan kuantitas yang salah?
Jadi sekarang, pertanyaan adalah, apa variansnya \(\hat{\beta}_0\) dan \(\hat{\beta}_1\)? Mereka adalah
\[ \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} \]
Ini mengukur variabilitas dari perkiraan karena peluang acak selama pengambilan sampel. Apakah ini “yang terbaik”? Apakah memungkinan untuk mendapatkan variansi kecil? Anda hanya perlu mempercayai apa yang kami katakan karena menunjukkan bahwa ini benar berada di luar cakupan kursus ini.
## 8.2 Distribusi Sampling
Sekarang kita telah “mendefinisikan ulang” perkiraan untuk \(\hat{\beta}_0\) dan \(\hat{\beta}_1\) sebagai variabel acak, kita dapat membahas distribusi samplingnya, yaitu distribusi ketika suatu statistik dianggap sebagai variabel acak
Karena \(\hat {\beta}_0\) dan \(\hat {\beta}_1\) adalah kombinasi linear dari \(Y_i\) dan setiap \(Y_i\) yang berdistribusi normal, maka \(\hat {\beta}_0\) dan \(\hat {\beta}_1\) juga berdistribusi normal.
Kemudian, gabungkan semua hal di atas, sampai pada distribusi \(\hat {\beta}_0\) dan \(\hat {\beta}_1.\)
Untuk \(\hat {\beta}_1\) kami katakan,
\[ \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). \]
Untuk \(\hat {\beta}_1\) kami katakan,
\[ \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 lebih singkatnya,
\[ \hat{\beta}_1 \sim N\left( \beta_1, \frac{\sigma^2}{S_{xx}} \right).\]
Dan untuk \(\hat{\beta}_0\)
\[ \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 lebih singkatnya,
\[ \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 kami telah lalai untuk membuktikan sejumlah hasil ini. Alih-alih bekerja melalui turunan yang membosankan dari distribusi sampling ini, kami akan membenarkan hasil ini untuk diri kami sendiri dengan menggunakan simulasi.
Catatan untuk pembaca saat ini: Derivasi dan bukti ini dapat ditambahkan ke lampiran di lain waktu. Anda 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 anda tidak akan diminta untuk melakukan penurunan atau derivasijenis ini dalam mata pelajaran ini, untuk saat ini dihilangkan.
### 8.2.1 Mensimulasikan Distribusi Sampling
Untuk memverifikasi hasil di atas, kami akan mensimulasikan sampel berukuran \(n=100\) dari model
\[ Y_i = \beta_0 + \beta_1 x_i + \epsilon_i\]
di mana \(\epsilon_i \sim N (0,\sigma^2).\) Dalam hal ini, parameternya dikenal sebagai:
\(\beta_0 = 3\)
\(\beta_1 = 6\)
\(\sigma^2 = 4\)
Kemudian, berdasarkan hal di atas, 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 sebelumnya berapa nilai \(x\) untuk simulasi ini, karena nilai \(x\) dalam SLR juga dianggap kuantitas yang diketahui. Pilihan nilai \(x\) bisa berubah-ubah. Di sini kami juga menetapkan pengacakan, dan menghitung \(S_ {xx}\) yang akan kami perlukan untuk selanjutnya.
set.seed(42)
sample_size = 100 # this is n
x = seq(-1, 1, length = sample_size)
Sxx = sum((x - mean(x)) ^ 2)
Kami juga memperbaiki nilai parameter kami
beta_0 = 3
beta_1 = 6
sigma = 2
Dengan informasi ini, kami mengetahui bahwa distribusi 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).\]
Itu 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}\]
Kami sekarang mensimulasikan data dari model ini 10.000 kali. Perhatikan bahwa ini bukan cara R dalam melakukan simulasi. Kami melakukan simulasi dengan cara ini dalam upaya untuk mendapatkan kejelasan. Misalnya, kita bisa menggunakan fungsi sim_slr() dari bab sebelumnya. Kami juga hanya menyimpan variabel di lingkungan global daripada 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 kali kami mensimulasikan data, kami memperoleh nilai koefisien yang diperkirakan. Variabel beta_0_hats dan beta_1_hats sekarang menyimpan 10.000 nilai simulasi dari \(\hat {\beta}_0\) dan \(\hat {\beta}_1\) masing-masing.
Pertama-tama, kami memverifikasi distribusi \(\hat {\beta}_1\).
mean(beta_1_hats) # rata-rata empiris
## [1] 6.001998
beta_1 # rata-rata aktual
## [1] 6
var(beta_1_hats) # variansi empiris
## [1] 0.11899
var_beta_1_hat # variansi aktual
## [1] 0.1176238
Kami melihat bahwa rata-rata dan varians dari empiris sangat mirip. Kami juga memverifikasi bahwa distribusi empiris normal. Untuk melakukannya, kami memplot histogram dari beta_1_hats, dan menambahkan kurva untuk distribusi \(\hat {\beta}_1\) yang sebenarnya. Kami menggunakan prob = TRUE untuk meletakkan histogram pada skala yang sama dengan kurva normal.
# Catatan untuk menggunakan 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)
Terlihat data berdistribusi normal.
Kami mengulangi proses untuk \(\hat {\beta}_0\)
mean(beta_0_hats) # rata-rata empiris
## [1] 3.001147
beta_0 # rata-rata aktual
## [1] 3
var(beta_0_hats) # variansi empiris
## [1] 0.04017924
var_beta_0_hat # Variansi aktual
## [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 studi simulasi ini, kami hanya mensimulasikan sejumlah sampel yang terbatas. Untuk benar-benar memverifikasi hasil distribusi, kita perlu mengamati jumlah sampel yang tak 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)
## 8.3 Kesalahan Standar atau Standard Error
Jadi sekarang kami mempercayai hasil dari dua distribusi
\[\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 membakukan atau menstandarisasikan hasil ini kami 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 mengetahui \(\sigma\) dalam praktiknya, kita harus memperkirakannya menggunakan \(s_e\), yang kita masukkan ke ekspresi yang ada untuk deviasi standar perkiraan kita.
Kedua ekspresi baru ini disebut kesalahan standar yang merupakan perkiraan deviasi standar 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 kesalahan standar, daripada dengan deviasi standar, kita memperoleh hasil berikut yang memungkinkan kita 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 perhatikan bahwa,
\[ \frac{\text{RSS}}{\sigma^2} = \frac{(n-2)s_e^2}{\sigma^2} \sim \chi_{n-2}^2.\]
Juga ingat bahwa variabel acak \(T\) didefinisikan sebagai,
\[T = \frac{Z}{\sqrt{\frac{\chi_{d}^2}{d}}}\]
mengikuti distribusi \(t\) dengan \(d\) derajat kebebasan, di mana \(\chi_{d}^2\) adalah \(\chi^2\) variabel acak dengan \(d\) derajat kebebasan.
Kamu menulis,
\[T \sim t_d\]
untuk mengatakan bahwa variabel acak \(T\) mengikuti distribusi \(t\) dengan derajat kebebasan \(d\).
Kemudian kami menggunakan trik klasik yaitu “kalikan dengan 1” dan menyusun kembali 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 distribusi \(t\) mirip dengan normal standar, tetapi dengan ekor (tails) yang lebih berat. Ketika derajat kebebasan meningkat, distribusi \(t\) menjadi semakin seperti normal standar. Di bawah ini kami memplot distribusi normal standar serta dua contoh distribusi a \(t\) dengan derajat kebebasan berbeda. Perhatikan bagaimana distribusi \(t\) dengan derajat kebebasan yang lebih besar lebih mirip dengan kurva normal standar.
# Tentukan kisi nilai x
x = seq(-4, 4, length = 100)
# Kurva plot untuk normal standar
plot(x, dnorm(x), type = "l", lty = 1, lwd = 2,
xlab = "x", ylab = "Density", main = "Normal vs t Distributions")
# Tambahkan kurva untuk distribusi t
lines(x, dt(x, df = 1), lty = 3, lwd = 2, col = "darkorange")
lines(x, dt(x, df = 10), lty = 2, lwd = 2, col = "dodgerblue")
# Tambahkan legenda
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"))
## 8.4 Interval Kepercayaan untuk Slope dan Intercept
Ingatlah bahwa interval kepercayaan untuk rata-rata sering kali berbentuk:
\[\text {EST} \pm \text {CRIT} \cdot \text {SE}\]
atau
\[\text {EST} \pm \text {MARGIN}\]
di mana \(\text {EST}\) adalah perkiraan untuk parameter yang diminati (interest), \(\text {SE}\) adalah kesalahan standar dari perkiraan, 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}}}\]
di mana \(t_{\alpha / 2, n - 2}\) adalah nilai kritis sehingga \(P(t_ {n-2}> t_{\alpha / 2, n - 2}) = \alpha / 2.\)
## 8.5 Tes Hipotesis
Ronald Aylmer Fisher
Ingat bahwa statistik pengujian \(\text {(TS)}\) untuk pengujian rata-rata sering kali berbentuk:
\[\text {TS} = \frac {\text {EST} - \text {HYP}} {\text {SE}}\]
di mana \(\text {EST}\) adalah perkiraan untuk parameter interest, \(\text {HYP}\) adalah nilai hipotesis dari parameter, dan \(\text {SE}\) adalah kesalahan standar dari perkiraan.
Jadi, untuk mengujinya
\[H_0: \beta_0 = \beta_{00} \quad \text{vs} \quad H_1: \beta_0 \neq \beta_{00}\]
Kami menggunakan statistik uji
\[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}}}}\]
yang, di bawah hipotesis nol, mengikuti distribusi \(t\) dengan \(n - 2\) derajat kebebasan. Kami menggunakan \(\beta_{00}\) untuk menunjukkan nilai hipotesis dari \(\beta_0\).
Begitu pula untuk menguji
\[H_0: \beta_1 = \beta_ {10} \quad \text {vs} \quad H_1: \beta_1 \neq \beta_ {10}\]
kami menggunakan statistik uji
\[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, di bawah hipotesis nol, mengikuti distribusi \(t\) dengan \(n - 2\) derajat kebebasan. Kami sekarang menggunakan \(\beta_{10}\) untuk menunjukkan nilai hipotesis dari \(\beta_1\) .
## 8.6 Contoh cars
Sekarang kita kembali ke contoh cars dari bab sebelumnya untuk mengilustrasikan konsep ini. Pertama-tama kami menyesuaikan model menggunakan lm() kemudian menggunakan summary() untuk melihat hasil secara lebih detail.
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
### 8.6.1 Tes di 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 informasi apa yang tersedia, dan kemudian kita menggunakan $ operator dan coefficients untuk mengekstrak informasi yang kita minati. Dua nilai di sini harusnya sudah kita kenal.
\[\hat{\beta}_0 = -17.5790949\]
dan
\[\hat{\beta}_1 = 3.9324088\]
yang merupakan perkiraan kami untuk parameter model \(\beta_0\) dan \(\beta_1\).
Sekarang mari kita fokus pada 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 error \(\hat {\beta}_1\),
\[\text {SE} [\hat {\beta}_1] = \frac {s_e} {\sqrt {S_{xx}}} = 0,4155128.\]
Nilai ketiga, t value, adalah nilai statistik uji 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 tes itu.
\[\text {p-value} = 1,4898365 \times 10 ^ {- 12}\]
Perhatikan di sini, kami secara khusus menguji apakah \(\beta_1 = 0\) atau tidak.
Baris pertama mengeluarkan 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, kode berikut menyimpan informasi summary(stop_dist_model)$coefficients dalam variabel baru stop_dist_model_test_info, lalu mengekstrak setiap elemen ke dalam variabel baru yang menjelaskan informasi yang dikandungnya.
stop_dist_model_test_info = summary(stop_dist_model)$coefficients
beta_0_hat = stop_dist_model_test_info[1, 1] # Estimasi
beta_0_hat_se = stop_dist_model_test_info[1, 2] # Std. Error
beta_0_hat_t = stop_dist_model_test_info[1, 3] # nilai t atau 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] # Estimasi
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 atau nilai t
beta_1_hat_pval = stop_dist_model_test_info[2, 4] # Pr(>|t|)
Kemudian kita dapat memverifikasi beberapa ekspresi penyetaraan: uji statistik \(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
### 8.6.2 Signifikansi Regresi, Uji-t
Kami 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 pikirkan tentang pengujian ini dalam kaitannya dengan model regresi linier sederhana,
\[Y_i = \beta_0 + \beta_1 x_i + \epsilon_i.\]
Jika kami mengasumsikan hipotesis nol benar, maka \(\beta_1 = 0\) dan kami memiliki model,
\[Y_i = \beta_0 + \epsilon_i.\]
Dalam model ini, respon tidak bergantung pada prediktor. Jadi kita bisa memikirkan tes ini dengan cara berikut,
Di bawah \(H_0\) tidak ada hubungan linier yang signifikan antara \(x\) dan \(y\).
Di bawah \(H_1\) ada hubungan linier yang signifikan antara \(x\) dan \(y\).
Untuk contoh cars
*Dibawah \(H_0\) tidak ada hubungan linier yang signifikan antara kecepatan dan jarak berhenti.
*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, kami akan menolak hipotesis nol pada tingkat \(\alpha\) yang masuk akal, misalnya \(\alpha = 0,01\). Jadi kami katakan ada hubungan linier yang signifikan antara kecepatan dan jarak berhenti. Perhatikan bahwa kami menekankan linier.
knitr::include_graphics("C:/Users/user/Downloads/plot.png")
Dalam plot data simulasi ini, kita melihat hubungan yang jelas antara \(x\) dan \(y\), namun ini bukan hubungan linier. Jika kita memasukkan garis ke data ini, itu sangatlah datar. Tes yang dihasilkan untuk \(H_0: \beta_1 = 0 \ vs \ H_1: \beta_1 \neq 0\) memberikan nilai-p yang besar, dalam hal ini \(0.7564548\), jadi kami 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”, tetapi untuk saat ini, sadari bahwa pengujian \(H_0: \beta_1 = 0 \ vs \ H_1: \beta_1 \neq 0\) ) hanya dapat mendeteksi hubungan garis lurus.
### 8.6.3 Interval Kepercayaan di R
Menggunakan R kita bisa mendapatkan interval kepercayaan untuk \(\beta_0\) and \(\beta_1\) dengan mudah
confint(stop_dist_model, level = 0.99)
## 0.5 % 99.5 %
## (Intercept) -35.706610 0.5484205
## speed 2.817919 5.0468988
Ini secara otomatis menghitung 99% interval keyakinan untuk \(\beta_0\) dan \(\beta_1\), baris pertama untuk \(\beta_0\), baris kedua untuk \(\beta_1\).
Untuk contoh cars saat menginterpretasikan interval ini, kami katakan, kami 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 keyakinan 99% ini tidak berisi 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.
Anda harus mencurigai interval keyakinan untuk \(\beta_0\), karena mencakup nilai negatif, yang berkoresponden dengan jarak henti yang negatif. Secara teknis, interpretasinya adalah bahwa kami yakin 99% bahwa jarak berhenti rata-rata sebuah mobil yang menempuh jarak 0 mil per jam adalah antara -35.7066103 dan 0.5484205 kaki, tetapi kami tidak begitu percaya, karena kami benar-benar yakin bahwa itu akan menjadi non-negatif.
Catatan, kita dapat mengekstrak nilai tertentu dari output ini dengan beberapa cara. Kode ini tidak berjalan, dan sebagai gantinya, anda harus memeriksa hubungannya dengan keluaran kode di atas.
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
Kami juga dapat memverifikasi bahwa kalkulasi yang dilakukan R untuk interval \(\beta_1\).
# perkiraan store
beta_1_hat = coef(stop_dist_model)[2]
# standard error store
beta_1_hat_se = summary(stop_dist_model)$coefficients[2, 2]
# Menghitung nilai kritis untuk 99% CI dua sisi
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
## 8.7 Interval Kepercayaan untuk Rata-Rata Respon
Selain interval keyakinan untuk \(\beta_0\) dan \(\beta_1\), ada dua perkiraan interval umum lainnya yang digunakan dengan regresi. Yang pertama disebut interval kepercayaan untuk respons rata-rata. Seringkali, kita menginginkan estimasi interval untuk rata-rata, \(E [Y \mid X = x]\) untuk nilai tertentu dari \(x\).
Dalam situasi ini kami menggunakan \(\hat {y} (x)\) sebagai perkiraan kami untuk \(E [Y \mid X = x]\). Kami mengubah sedikit notasi kami 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 baik karena tidak bias:
\[\text {E} [\hat {y} (x)] = \beta_0 + \beta_1 x.\]
Kami 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. Karena \(\hat {\beta}_0\) dan \(\hat {\beta}_1\) adalah kombinasi linier dari variabel acak normal, \(\hat {y} (x)\) juga.
\[\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, karena kita perlu memperkirakan varians ini, kita sampai pada kesalahan standar (standard error)dari perkiraan,
\[\text {SE} [\hat {y} (x)] = s_e \sqrt {\frac {1} {n} + \frac {(x- \bar {x})^2} {S_{xx }}}.\]
Kami 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 menggunakan R, kami menggunakan fungsi predict(). Kami memberikan fungsi model yang pas serta data baru, disimpan sebagai bingkai data. (Ini penting, agar R mengetahui nama variabel prediktor.) Di sini, kita menemukan interval kepercayaan untuk jarak henti rata-rata saat mobil berjalan 5 mil per jam dan saat mobil berjalan 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
## 8.8 Interval Prediksi untuk Pengamatan Baru
Terkadang kami menginginkan estimasi interval untuk observasi baru, \(Y\), untuk nilai tertentu dari \(x\). 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 kami 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 mengenai garis regresi yang sebenarnya menurut distribusi \(N (0, \sigma^2)\). Karena itu, kami menambahkan faktor tambahan \(\sigma^2\) ke variabilitas perkiraan kami untuk menjelaskan variabilitas pengamatan tentang garis regresi.
\[\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}}}\]
Kami 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 pemberitahuan R 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 respons rata-rata.
## 8.9 Band Keyakinan dan Prediksi
Seringkali kami ingin memplot kedua interval kepercayaan untuk respons rata-rata dan interval prediksi untuk semua kemungkinan nilai \(x\). Kami menyebutnya band kepercayaan dan 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:
Kami menggunakan argumen ylim untuk meregangkan sumbu \(y\) dari plot, karena pita atau band diperpanjang lebih jauh daripada titik.
Kami menambahkan titik di titik \((\bar {x}, \bar {y})\).
Ini adalah titik yang akan selalu dilewati garis regresi. (Pikirkan alasannya.)
Ini adalah titik di mana band atau pita keyakinan dan prediksi menjadi yang tersempit. Lihatlah kesalahan standar keduanya untuk memahami mengapa.
Pita prediksi (titik-titik biru) kurang melengkung dibandingkan dengan pita kepercayaan (biru putus-putus). Ini adalah hasil dari faktor tambahan \(\sigma^2\) yang ditambahkan ke varian pada nilai \(x\).
## 8.10 Signifikansi Regresi, Uji-F
Dalam kasus regresi linier sederhana, uji \(t\) untuk signifikansi regresi setara dengan uji lain, uji \(F\) untuk signifikansi regresi. Kesetaraan ini hanya akan benar untuk regresi linier sederhana, dan di bagian selanjutnya kita hanya akan menggunakan uji \(F\) untuk signifikansi regresi.
Ingat dari bagian terakhir dekomposisi varians yang kita lihat sebelum menghitung (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 tes \(F\), kami akan menyusun informasi ini dalam tabel ANOVA,
knitr::include_graphics("C:/Users/user/Downloads/tabel.png")
ANOVA, atau Analisis Varians akan menjadi konsep yang sering kita bahas dalam kursus ini. Untuk saat ini, kita akan fokus pada hasil tabel, yaitu statistik \(F\),
\[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\) di bawah hipotesis nol. Distribusi \(F\) adalah distribusi kontinu yang hanya mengambil nilai positif dan memiliki dua parameter, yaitu dua derajat kebebasan.
Ingat, dalam signifikansi uji regresi, \(Y\) tidak bergantung pada \(x\) dalam hipotesis nol.
\[H_0: \beta_1 = 0 \quad \quad Y_i = \beta_0 + \epsilon_i\]
Sedangkan ada kemungkinan pada hipotesis alternatif \(Y\) bergantung pada \(x\).
\[H_1: \beta_1 \neq 0 \quad \quad Y_i = \beta_0 + \beta_1 x_i + \epsilon_i\]
Kita dapat menggunakan statistik \(F\) 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, kami akan menolak hipotesis nol ketika statistik \(F\) adalah besar, yaitu, ketika ada kemungkinan rendah bahwa pengamatan bisa datang dari model hipotesis nol secara kebetulan. Kami akan membiarkan R menghitung nilai-p.
Untuk melakukan tes \(F\) di R anda dapat melihat baris terakhir dari output dari summary() yang disebut F-statistic yang memberikan nilai uji statistik, derajat kebebasan yang relevan, serta uji nilai p.
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, anda dapat menggunakan fungsi anova() untuk menampilkan informasi dalam 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 untuk pengujian. Anda harus memperhatikan bahwa nilai-p dari uji \(t\) adalah sama. Anda mungkin juga memperhatikan bahwa pengujian nilai statistik untuk uji \(t\),\(9.46399\), dapat dikuadratkan untuk mendapatkan nilai \(F\) statistik, \(89.5671065.\)
Perhatikan bahwa ada cara lain yang setara untuk melakukan ini di R, yang akan sering kita lihat 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 cars. Perhatikan bahwa \(\hat {y} = \bar {y}\) saat \(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.\)
Kami kemudian dapat menganggap penggunaan anova() ini sebagai perbandingan langsung kedua model. (Perhatikan bahwa kita mendapatkan nilai p yang sama lagi.)