Pengenalan

Regresi linier sederhana digunakan untuk memprediksi hasil kuantitatif y berdasarkan satu variabel prediktor tunggal x. Tujuannya adalah untuk membangun model matematika yang mendefinisikan y sebagai fungsi dari variabel x.Model ini dapat digunakan untuk memprediksi hasil di masa mendatang berdasarkan nilai x baru.

Untuk memeriksa asumsi regresi, kita perlu memeriksa distribusi residual. Reregsi linier memuat beberapa asumsi tentang data sebagai berikut:
- Liniearitas Data : hubungan antara prediktor(x) dan y yang diasumsikan linier.
- Normalitas Linier : kesalahan residual diasumsikan terdistribusi normal.
- Homogenitas Varians Residual : residu diasumsikan memiliki varians yang konstan (homoskedastisitas)
- Independensi yang memiliki istilah kesalahan residual.

Kita perlu memeriksa apakah asumsi di atas benar atau tidak. Dan potensi masalah meliputi:
- Non-linearitas : hasil hubungan prediktor
- Heteroskedastisitas : varians istilah kesalahan yang tidak konstab
- Adanya nilai-nilai berpengaruh dalam data yang dapat berupa:
Pencilan :nilai ekstrim dalam variabel hasil (y)
Poin Leverage Tinggi : nilai ekstrim dalam variabel prediktor (x)

Catatan: Kita bisa memeriksa semua asumsi dan masalah potensial ini dengan membuat beberapa plot diagnostik yang memvisualisasikan kesalahan residual.

\[\begin{align*} y_i= \beta_1 + \beta_2xi + \epsilon_i \end{align*}\]

Prediksi dan Meramalkan

Dengan asumsi bahwa nilai ekspetasi dari error term dalam persamaan 1 adalah nol, persamaan 2 diberikan \(\hat{y}i\), nilai prediksi dari ekspetasi \(y_i\) yang diberikan \(x_i\), dimana \(b_1\) dan \(b_2\) adalah estimasi (kuadrat terkecil) dari parameter regresi \(\beta_1\) dan \(\beta_2\).

\[\begin{align*} \hat{y}i= \beta_1 + \beta_2xi \end{align*}\]

Nilai prediksi \(\hat{y}i\) adalah variabel acak tergantung pada sampel, oleh karena itu, kita bisa menghitung interval kepercayaan dan menguji hipotesis dan dapat menentukan distribusi dan variansnya. Prediksi memiliki distribusi normal, kombinasi linier dari dua variabel acak berdistribusi normal \(b_1\) dan \(b_2\), dan variansnya diberikan pada persamaan 3.
seperti sebelumnya, karena kita perlu menggunakan estimasi varians, kita menggunakan distribusi \(t-\) dan bukan distribusi normal. Persamaan 3 berlaku untuk setiap \(x\) yang diberikan, katakanlah \(x_0\), tidak hanya untuk x dalam kumpulan data.

\[\begin{align*} \hat{var(f_i)}= \hat{\sigma}^2[1+\frac{1}{N}+\frac{(x_i-\overline{x})^2}{\sum N_j=1 (x_j-\overline{x})^2}] \end{align*}\]

Yang akan direduksi menjadi,

\[\begin{align*} \hat{var(f_i)}= \hat{\sigma}^2 +\frac{\hat{\sigma}^2}{N}+(x_i-\overline{x})^2 * \hat{var(b_2)} \end{align*}\]

Catatan: Varians pada persamaan 3 tidak sama dengan varians pada persamaan 4. Pada yang pertama, varians dari perkiraan ekspetasi \(y\), sedangkan yang terakhir adalah varians varians dari kejadian tertentu \(y\). Bisa kita sebut varians terakhir sebagai varians dari kesalahan perkiraan. varians kesalahan memiliki ramalah yang lebih besar dari varians prediksi \(\ E(y \mid x)\)
Mari kita tentukan kesalahan standar untuk persamaan makanan untuk pendapata rumah tangga, (pendapatan) 2000 dollar seminggu yaitu pada \(x=x_0=20\), dengan persamaan 4. Untuk melakukannya, kita perlu mengambil \(var(b_2)\) dam \(\hat{\sigma}\), kesalahan standar regresi dari keluaran regresi.

library(PoEdata)
data("food")
alpha <- 0.05
x <- 20
xbar <- mean(food$income)
linmod <- lm(food_exp~income, data=food)
b1 <- coef(linmod)[[1]]
b2 <- coef(linmod)[[2]]
yhatx <- b1+b2*x
slinmod <- summary(linmod)
df <- df.residual(linmod)
tcr <- qt(1-alpha/2, df)
N <- nobs(linmod)
N <- NROW(food)
varb2 <- vcov(linmod)[2,2]
sighat2 <- slinmod$sigma^2
varf <- sighat2+sighat2/N+(x-xbar)^2*varb2
sef <- sqrt(varf)
lb <- yhatx-tcr*sef
ub <- yhatx+tcr*sef
print(c(lb,ub))

Hasilnya adalah interval kepercayaan untuk ramalan [104.1323, 471.0854], yang seperti yang diharapkan, lebih besar dari interval kepercayaan dari perkiraan nilai y yang diharapkan berdasarkan Persamaan 4.
Mari kita hitung interval kepercayaan ramalan untuk semua pengamatan dalam sampel dan menggambar batas atas dan bawah bersama-sama dengan garis regresi. Gambar berikut menunjukkan pita interval kepercayaan tentang garis regresi.

sef <- sqrt(sighat2+sighat2/N+(food$income-xbar)^2*varb2)
yhatv <- fitted.values(linmod)
lbv <- yhatv-tcr*sef
ubv <- yhatv+tcr*sef
xincome <- food$income
dplot <- data.frame (xincome,yhatv, lbv, ubv)
dplotord <- dplot[order(xincome), ]
xmax <- max(dplotord$xincome)
xmin <- min(dplotord$xincome)
ymax <- max(dplotord$ubv)
ymin <- min(dplotord$lbv)
plot(dplotord$xincome, dplotord$yhatv,
     xlim=c(xmin, xmax),
     ylim=c(ymin, ymax),
     xlab="Income", ylab="Food Expenditure",
     type="l")
lines(dplotord$ubv~dplotord$xincome, lty=2)
lines(dplotord$lbv~dplotord$xincome, lty=2)

Cara lain untuk menemukan taksiran titik dan interval untuk prediksi \(E(y|x)\) dan prakiraan \(y\) (silakan lihat perbedaan yang saya sebutkan di atas) adalah dengan menggunakan fungsi **predict()** di R. Fungsi ini mensyaratkan bahwa nilai-nilai independen variabel di mana prediksi (atau ramalan) dimaksudkan memiliki struktur kerangka data. Contoh berikutnya menunjukkan pada titik paralel dan perkiraan interval pengeluaran makanan yang diprediksi dan diperkirakan untuk pendapatan adalah $2000. Seperti yang telah saya tunjukkan sebelumnya, perkiraan titik sama untuk prediksi dan ramalan, tetapi perkiraan interval sangat berbeda.

incomex=data.frame(income=20)
predict(linmod, newdata=incomex, interval="confidence", level=0.95)
predict(linmod, newdata=incomex, interval="prediction", level=0.95)

hasil

Sekarang mari kita gunakan fungsi **predict()** untuk mereplikasi garis estimasi interval, titik-titik dalam kumpulan data.

xmin <- min(food$income)
xmax <- max(food$income)
income <- seq(from=xmin, to=xmax)
ypredict <- predict(linmod, newdata=data.frame(income), interval="confidence")
yforecast <- predict(linmod, newdata=data.frame(income), interval="predict")
matplot(income, cbind(ypredict[,1], ypredict[,2], ypredict[,3],     yforecast[,2], yforecast[,3]),
        type="l", lty=c(1, 2, 2, 3, 3),
        col=c("red", "light blue", "light blue", "orange", "orange"),
        ylab="Food Expenditure", xlab="Income")
points(food$income, food$food_exp)
legend("topleft",
       legend=c("E[y|x]", "lwr_pred", "upr_pred", "lwr_fercst", "upr_forcst"),
       lty=c(1, 2, 2, 3, 3),
       col=c("red", "light blue", "light blue", "orange", "orange"))

Gambar di atas menyajikan garis yang diprediksi dan yang diperkirakan pada grafik yang sama, untuk menunjukkan bahwa mereka memiliki perkiraan titik yang sama (garis merah, garis solid) dan bahwa garis yang diperkirakan jauh lebih besar daripada yang diprediksi. Dengan kata lain, kita mungkin berpikir tentang perbedaan antara dua jenis interval yang kami sebut prediksi dan ramalan sebagai berikut:
- Interval prediksi tidak seharusnya mencakup,bisa kita katakan 95% poin tetapi untuk memasukkan garis regresi, \(E(y|x)\), dengan probabilitas 95 persen.
- Interval yang diperkirakan, di sisi lain, harus mencakup setiap titik sebenarnya dengan probabilitas 95 persen.

Koefesien Determinasi

Tujuan penting dari analisis data adalah untuk menjelaskan variabilitas dalam variabel dependen \(y\). Untuk membandingkan model, kita perlu mempertimbangkan situasi dasar di mana kita memiliki informasi tentang variabel dependen saja. Dalam situasi ini, kita perlu memodelkan mean dari \(y\), \(\hat{y}\), untuk menjelaskan variabilitas, yang dianggap sebagai model dasar atau baseline. Namun, ketika kita memiliki informasi lain yang terkait dengan variabel dependen, seperti variabel independen kontinu \(x\), kita dapat mencoba mengurangi variabilitas variabel dependen (tidak dapat dijelaskan) dengan menambahkan \(x\) ke model dan kita memodelkan garis regresi linier yang dibentuk oleh nilai prediksi oleh \(x\), \(\hat{y}\).
Diharapkan, kita dapat mengharapkan pengurangan variabilitas y yang efektif sebagai hasil dari penambahan \(x\) ke dalam model. Oleh karena itu, perlu adanya evaluasi terhadap garis regresi, seperti “seberapa jauh lebih baik daripada model dasar yang hanya menggunakan rata-rata variabel terikat”.
Variasi total \(y\) terhadap mean sampelnya, Sum of Squares of Total (SST), dapat diuraikan menjadi variasi terhadap garis regresi, Sum of Squares of Errors (SSE), dan variasi garis regresi terhadap mean dari \(y\), Sum of Squares due to Segression (SSR), seperti yang ditunjukkan oleh Persamaan 5.

\[\begin{align*} STT=SSR+SSE \end{align*}\]

Jika Anda membutuhkan jumlah kesalahan kuadrat(SSE) atau jumlah kuadrat karena regresi(SSR) gunakan fungsi anova, yang strukturnya ditunjukkan pada Tabel 1.

library(stats)
library(kableExtra)
library(DT)
anov <- anova(linmod)
SSE <- anov[1,2]
SSR <- anov[1,2]
SST <- anov[1,2]+anov[2,2]
dfr <- data.frame(anov)
datatable(dfr,
      caption = htmltools::tags$caption(
            style = 'caption-side: bottom; text-align: center;',
            htmltools::em('Table 1:Output generated by the anova function.')),
          options = list(dom = 't'))

Hasil ANOVA pada Tabel 1 mencakup informasi berguna lainnya: jumlah derajat kebebasan, \(anov[2,1]\) dan varians yang diperkirakan \(\sigma^2=anov[2,3]\). Hasil ini juga menginformasikan, koefisien determinasi \(R^2\), didefinisikan sebagai proporsi varians dalam \(y\) yang dijelaskan oleh regresi, (SSR), dalam variasi total \(y\), (SST). Membagi kedua ruas Persamaan 5 dengan (SST) dan menyusun ulang suku-sukunya memberikan rumus untuk menghitung \(R^2\), seperti yang ditunjukkan pada Persamaan 6.

\[\begin{align*} R^2=\frac{SSR}{SST}=1-\frac{SSE}{SST} \end{align*}\]

\(R^2\) mengambil nilai antara 0 dan 1, dengan nilai yang lebih tinggi menunjukkan kecocokan yang lebih dekat dari garis regresi ke data (nilai R-kuadrat yang besar menunjukkan kecocokan yang baik). Dalam R, nilai \(R^2\) dapat diambil dari ringkasan model regresi dengan nama adj.r.squared; misalnya, dalam contoh makanan kita, \(R^2=0,3688\). Itu juga dicetak sebagai bagian dari ringkasan model regresi, seperti yang ditunjukkan oleh urutan kode berikut.

slinmod
(rsq <-slinmod$adj.r.squared)

hsl

Evaluasi Kemiringan Individu

Pada regresi linier sederhana dengan hanya satu \(x\), hasil uji-F global dan Uji-T (signifikansi kemiringan) memiliki kesimpulan yang sama. Uji formal signifikansi kemiringan adalah sebagai berikut:

\[\begin{align*} H_0:\beta_1=0 \end{align*}\] dan \[\begin{align*} H_a:\beta_1\neq 0 \end{align*}\]

Uji statistik: \(t=\frac{b_1}{s\epsilon_b1}\), dimana \(se_b1\) adalah standar error \(b_1\). Test ini menanyakan apakah ada hubungan linier yang signifikan secara statistik antara variabel dependen dan independen. Jika jawabannya tidak, maka model regresi tidak baik karena asumsi dasarnya adalah hubungan linier antara dua variabel. Jika tidak ada hubungan linier, maka kemiringan mendekati nol dan scatterplot y dan x diharapkan menunjukkan sebaran acak, yang berarti tidak ada korelasi.
Cara menemukan P-value di R, kita bisa menjalankan konsol sebagai berikut:

anova(linmod)$'Pr(>F)'[1]

Dari hasil di atas, kita dapat mengetahui bahwa P-value = 1.945862e-05, lebih kecil dari taraf signifikansi nya yaitu 0.05. Oleh karena itu, Hipotesis 0 ditolak atau \(H_0\) ditolak bahwa kemiringan nol dan ada hubugan linier yang signifikan antara y dan x.

Nilai dan Residual yang Terpasang

Model regresi linier membuat asumsi bahwa hubungan antara prediktor (\(x\)) dan variabel hasil adalah linier. Ini mungkin tidak benar. Hubungannya bisa polinomial atau logaritmik. Selain itu, data mungkin berisi beberapa pengamatan yang berpengaruh, seperti outlier (atau nilai ekstrim), yang dapat mempengaruhi hasil regresi. Oleh karena itu, Anda harus mendiagnosa model regresi yang Anda bangun untuk mendeteksi potensi masalah dan untuk memeriksa apakah asumsi yang dibuat oleh model regresi linier terpenuhi atau tidak.
Untuk melakukannya, kami biasanya memeriksa distribusi kesalahan residual, yang dapat memberi tahu kita lebih banyak tentang data kita. Di R, kita dapat dengan mudah menambah data kita untuk menambahkan nilai dan residu yang sesuai dengan menggunakan fungsi augment() [broom package]. Sebut saja model keluaran.diag.metrics karena berisi beberapa metrik yang berguna untuk diagnostik regresi. Kami akan menjelaskan tema nanti.

library(broom)
library(DT)
model.diag.metrics <- augment(linmod)
datatable(model.diag.metrics,
          caption = htmltools::tags$caption(
            style='caption-side: bottom; text-align: center;',
            htmltools::em('Table 2: The distribution of residuals errors.')),
          extensions = 'FixedColumns',
          options = list(scrollX = TRUE, fixedColumns = TRUE)
          )

Selanjutnya, mari kita visualisasikan kesalahan residual (dalam warna merah) antara nilai yang diamati dan garis regresi yang dipasang. Setiap segmen merah vertikal mewakili kesalahan residual antara nilai penjualan yang diamati dan nilai prediksi (yaitu pas) yang sesuai.

library(ggplot2)
ggplot(model.diag.metrics, aes(income, food_exp)) +
  geom_point() +
  stat_smooth(method = lm, se= FALSE) +
  geom_segment(aes(xend = income, yend = .fitted), color = "red", size = 0.3)

grafikk

Diagnostik Regresi

Diagnostik Regresi

Dalam statistik, diagnostik regresi adalah salah satu dari serangkaian prosedur yang tersedia untuk analisis regresi yang berusaha menilai validitas model dalam beberapa cara yang berbeda. Plot diagnostik menunjukkan residu dalam empat cara berbeda:
- Residuals vs Fitted: Digunakan untuk memeriksa asumsi hubungan linier. Garis horizontal, tanpa pola yang berbeda merupakan indikasi untuk hubungan linier.
- Normal Q-Q: Digunakan untuk memeriksa apakah residual berdistribusi normal. Ada baiknya jika titik residu mengikuti garis putus-putus lurus. Scale-Location (atau Spread-Location): Digunakan untuk memeriksa homogenitas varians dari residual (homoscedasticity). Garis horizontal dengan titik penyebaran yang sama merupakan indikasi yang baik dari homoskedastisitas.
- Residuals vs Leverage: Digunakan untuk mengidentifikasi kasus-kasus yang berpengaruh, yaitu nilai-nilai ekstrim yang mungkin mempengaruhi hasil regresi ketika dimasukkan atau dikeluarkan dari analisis. Plot ini akan dijelaskan lebih lanjut di bagian selanjutnya.
Plot diagnostik regresi dapat dibuat menggunakan fungsi dasar R plot() atau fungsi autoplot() [paket ggfortify], yang membuat grafik berbasis ggplot2:

library(dplyr)
library(ggfortify)
par(mfrow = c(2, 2))
autoplot(linmod)

grafik

Disini kita akan belajar cara menggunakan grafik dan metriks ini untuk memeriksa asumsi regresi dan untuk mendiagnosis masalah potensial dalam model.

Residual vs Fitted

Idealnya, plot residual tidak akan menunjukkan pola yang pas. Artinya, garis merah harus kira-kira horizontal di nol. Kehadiran pola dapat menunjukkan masalah dengan beberapa aspek model linier.
Perhatikan bahwa, jika plot residual menunjukkan hubungan non-linier dalam data, maka pendekatan sederhana adalah dengan menggunakan transformasi non-linier dari prediktor, seperti \(log(x)\), \((√x)\) dan \(x^2\), dalam model regresi .
Contoh berikut mengilustrasikan bagaimana residual terlihat ketika bentuk fungsional linier digunakan ketika hubungan sebenarnya adalah kuadrat. Gambar berikut menunjukkan residual dari estimasi model ekonometrik linier yang tidak ditentukan dengan benar ketika spesifikasi yang benar harus kuadrat.

set.seed(12345)
x <- runif(1000, -2.5, 2.5)
e <- rnorm(1000, 0, 4)
y <- 15-4*x^2+e
quadmod <- lm(y~x)
ehat <- resid(quadmod)
ymi <- min(ehat)
yma <- max(ehat)
plot(x, ehat, ylim=c(ymi, yma), xlab="x", ylab="residuals", col="yellow")

grafik

Normal QQ

Plot QQ residual dapat digunakan untuk memeriksa asumsi normalitas secara visual. Plot probabilitas normal dari residual kira-kira harus mengikuti garis lurus. Dalam contoh kita, semua titik jatuh kira-kira di sepanjang garis referensi ini, sehingga kita dapat mengasumsikan normalitas.
Asumsi lain yang ingin kami uji adalah normalitas residual, yang menjamin pengujian hipotesis yang andal dan interval kepercayaan bahkan dalam sampel kecil. Asumsi ini dapat dinilai dengan memeriksa histogram residu di bawah ini:

library(quantmod)
library(tseries)
ehat <- resid(linmod)
ebar <- mean(ehat)
sde <- sd(ehat)
hist(ehat, col="yellow", freq=FALSE, main="", ylab="density", xlab="ehat")
curve(dnorm(x, ebar, sde), col=2, add=TRUE, ylab="density", xlab="ehat")

grafik

Sementara histogram di masa depan di atas mungkin tidak mendukung satu kesimpulan atau lainnya tentang normalitas ehat, kita dapat melakukan uji Jarque-Bera, yang hipotesis nolnya adalah “Deret terdistribusi normal”.
Jika p-value \(<\alpha\) , maka hipotesis nol ditolak, yang berarti deret tersebut gagal dalam uji normalitas. Tes Jarque-Bera memerlukan penginstalan dan pemuatan paket tseries di R.

jarque.bera.test(ehat)

Scale Location

Plot ini menunjukkan jika residual tersebar merata di sepanjang rentang prediktor. Ada baiknya jika kita melihat garis horizontal dengan titik penyebaran yang sama. Dalam contoh kita, ini tidak terjadi. Dapat dilihat bahwa variabilitas (varians) dari titik residual meningkat dengan nilai variabel hasil yang dipasang, menunjukkan varians yang tidak konstan dalam kesalahan residual (atau heteroskedastisitas).

linlogmod <- lm(log(food_exp)~income, data=food)
plot(linlogmod, 3)

grafik

Residu vs Leverage

• Pencilan : Pencilan adalah titik yang memiliki nilai variabel hasil ekstrem. Kehadiran outlier dapat mempengaruhi interpretasi model, karena meningkatkan RSE. Pencilan dapat diidentifikasi dengan memeriksa sisa standar (atau sisa siswa), yang merupakan sisa dibagi dengan kesalahan standar yang diperkirakan. Residu terstandarisasi dapat diartikan sebagai jumlah kesalahan standar yang menjauhi garis regresi. Pengamatan yang residu standarnya lebih besar dari 3 dalam nilai absolut adalah kemungkinan outlier (James et al. 2014).
  
• Poin leverage tinggi: Titik data memiliki leverage tinggi, jika memiliki nilai \(x\) prediktor ekstrim. Ini dapat dideteksi dengan memeriksa statistik leverage atau nilai topi. Nilai statistik ini di atas \(\frac{2(p+1)}{n}\) menunjukkan observasi dengan leverage yang tinggi (P. Bruce dan Bruce 2017); dimana, \(p\) adalah jumlah prediktor dan n adalah jumlah observasi.
  Dalam contoh kami, plot di atas menyoroti 3 titik paling ekstrem teratas (#31, #38 dan #39). Namun, tidak ada outlier yang melebihi 3 standar deviasi, yang menurut (James et al. 2014).
  Sekarang mari kita tunjukkan contoh lain, di mana data berisi dua nilai ekstrem dengan pengaruh potensial pada hasil regresi:

high_lev <- data.frame(
          x = c(food$income, 500, 700),
          y = c(food$food_exp, 8000, 10000)
)
high_lev_mod <- lm(y~x, high_lev)
par(mfrow = c(1, 2))
plot(high_lev_mod, 4)
plot(high_lev_mod, 5)

grafik

Pada contoh di atas, dua titik data berada jauh di luar garis jarak Cook. Residu lainnya tampak berkerumun di sebelah kiri. Plot mengidentifikasi pengamatan yang berpengaruh sebagai #41 dan #42.

Model Kuadrat

Model pada bagian sebelumnya terlihat bagus, tetapi kita dapat melihat bahwa plot memiliki kelengkungan yang tidak dijelaskan dengan baik oleh model linier. Sekarang kita cocok dengan model pendapatan kuadrat. Kami membuat variabel yang disebut pendapatan2 yang merupakan kuadrat dari pendapatan variabel.

\[\begin{align*} food_exp = b_1 + b_2x^2_i+e_i \end{align*}\]

library(PoEdata)
data("food")
food$income2 <- food$income^2
quadmod1 <- lm(food_exp~income+income2, data=food)
summary(quadmod1)

result

Hasil dari model kuadratik kami di atas ternyata tidak sesuai dengan yang kami harapkan. Sekarang, mari kita lihat ketika kita memiliki contoh yang baik untuk model kuadrat.

data <- data.frame(hours=c(6, 9, 12, 14, 30, 35, 40, 47, 51, 55, 60),
                   happiness=c(14, 28, 50, 70, 89, 94, 90, 75, 59, 44, 27))
data$hours2 <- data$hours^2
quadmod2 <- lm(happiness ~ hours +hours2, data=data)
summary(quadmod2)

result

Model Polinomial

Mari kita coba menerapkan model polinomial ke kumpulan data makanan:

\[\begin{align*} food_exp=b_1+b_2x^3_i+e_i \end{align*}\]

food$income2 <- food$income^2
polymod1 <- lm(food_exp~income+I(income2), data=food)
summary(polymod1)

result

Dari ringkasan polinomial di atas, kita tahu bahwa model tidak sesuai dengan yang kita harapkan. Oleh karena itu, kami tidak tersedia untuk menerapkan model ini ke kumpulan data kami (food).
Selanjutnya, mari kita lihat model regresi lain yang mungkin menyertakan istilah kuadrat atau kubik untuk lebih menggambarkan sifat data. Berikut ini adalah contoh model kuadrat dan kubik untuk dataset wa_wheat, yang memberikan hasil gandum tahunan dalam ton per hektar di Greenough Shire di Australia Barat selama periode 48 tahun.

library(PoEdata)
data("wa_wheat")
wa_wheat$time <- wa_wheat$time^2
polymod2 <- lm(greenough~time+I(time^2), data=wa_wheat)
summary(polymod2)

result

Model Linear-log

Bentuk fungsional non-linier dari model regresi berguna ketika hubungan antara dua variabel tampaknya lebih kompleks daripada yang linier. Solusi yang mungkin untuk memecahkan masalah ini adalah dengan menggunakan bentuk fungsional non-linier berdasarkan model matematika, penalaran, atau hanya memeriksa plot pencar dari data. Fungsi logaritmik cocok dengan profil ini dan, ternyata, relatif mudah untuk ditafsirkan, yang membuatnya sangat populer dalam model ekonometrik.
Bentuk umum dari model ekonometrika log linier diberikan pada Persamaan 8.

\[\begin{align*} y_i=\beta_1+\beta_2*log(x_i)+\epsilon_i \end{align*}\]

Efek marjinal dari perubahan \(x\) pada \(y\) adalah kemiringan kurva regresi dan diberikan oleh Persamaan 11; tidak seperti dalam bentuk linier, itu tergantung pada \(x\) dan oleh karena itu, hanya berlaku untuk perubahan kecil di \(x\).

\[\begin{align*} \frac{dy}{dx}=\frac{\beta_2}{x} \end{align*}\]

Terkait dengan model log linier, ukuran lain yang menarik dalam ekonomi adalah semi-elastisitas \(y\) terhadap \(x\). Semi-elastisitas menunjukkan bahwa perubahan \(x\) sebesar 1% mengubah \(y\) sebesar \(\frac{\beta_2}{100%}\) unit $y%. Karena semi-elastisitas juga berubah ketika \(x\) berubah, itu seharusnya hanya ditentukan untuk perubahan kecil pada \(x\).

\[\begin{align*} dy=\frac{\beta_2}{100}(%\Delta x) \end{align*}\]

Besaran lain yang mungkin menarik adalah elastisitas \(y\) terhadap \(x\) dan menunjukkan bahwa satu persen peningkatan \(x\) menghasilkan \((\frac{\beta_2}{y})\) persen perubahan \(y\).

\[\begin{align*} dy=\frac{\beta_2}{100}(%\Delta x) \end{align*}\]

Mari kita memperkirakan model log linier untuk dataset makanan, menggambar kurva regresi, dan menghitung efek marjinal untuk beberapa nilai tertentu dari variabel dependen.

library(xtable)
linlogmod <- lm(food_exp~log(income), data=food)
tbl <- data.frame(xtable(linlogmod))
tbl

result

b1 <- coef(linlogmod)[[1]]
b2 <- coef(linlogmod)[[2]]
pred_linlogmod <- predict(linlogmod, newdata=data.frame(income), interval="confidence")
plot(food$income, food$food_exp, xlab="Income", ylab="Food Expenditure")
lines(pred_linlogmod[,1]~income, lty=1, col="black")
lines(pred_linlogmod[,2]~income, lty=2, col="red")
lines(pred_linlogmod[,3]~income, lty=2, col="red")

result

x <- 10
y <- b1+b2*log(x)
DyDx <- b2/x
DyPDx <- b2/100
PDyPDx <- b2/y
print(c(DyDx, DyPDx, PDyPDx))

Catatan: Hasil dari DyDx adalah 13.216584, hasil dari DyPDx adalah 1.321658, hasil dari PDyPDx adalah 0.638061.
Hasil untuk pendapatan sebesar 1000 dollar adalah sebagai berikut: \(dy/dx=13.217\), yang menunjukkan bahwa peningkatan pendapatan sebesar 100 dollar (yaitu, satu unit x ) meningkatkan pengeluaran sebesar 13.217 dollar ; untuk peningkatan pendapatan sebesar 1%, yaitu peningkatan sebesar 10 dollar, pengeluaran meningkat sebesar 1,322 dollar ; dan, akhirnya, untuk peningkatan 1% dalam pengeluaran pendapatan meningkat sebesar 0,638%.

Model Log-Linear

Mengubah variabel dependen dengan fungsi \(log()\) berguna ketika variabel memiliki distribusi yang miring, yang umumnya terjadi dengan jumlah yang tidak boleh negatif. Transformasi \(log()\) seringkali membuat distribusi mendekati normal.

\[\begin{align*} log(y_i)=\beta_1+\beta_2xi+\epsilon_i \end{align*}\]

Rumus berikut dengan mudah diturunkan dari persamaan log-linier di atas. Semi-elastisitas di sini memiliki interpretasi yang berbeda dari yang ada dalam model log-linier: di sini, peningkatan x sebesar satu unit (dari x) menghasilkan perubahan \(100b_2\) persen dalam y. Untuk perubahan kecil pada x, jumlah 100b2 dalam model log-linier juga dapat diinterpretasikan sebagai laju pertumbuhan dalam y (sesuai dengan peningkatan unit dalam x). Misalnya, jika x adalah waktu, maka \(100b_2\) adalah laju pertumbuhan dalam y per satuan waktu.
- Prediksi : \(\hat{y}n=exp(b1+b2x)\) atau \(\hat{y}_c=exp(b1+b2x+\frac{\sigma}{2})\) dengan priduktor alami \(\hat{y}n\) untuk digunakan dalam sampel kecil dan prediktor “dikoreksi”, \(\hat{y}_c\) dalam sampel besar
- Efek marginal(kemiringan) : \(\frac{dy}{dx}=b_2y\)
- Semi-elastisitas : \(%\Delta y=100b_2\Delta x\)
Menggunakan dataset wa_wheat untuk menghitung persamaan hasil.

library(xtable)
library(kableExtra)
loglinmod1 <- lm(log(greenough)~time, data=wa_wheat)
sloglinmod1 <- summary(loglinmod1)
tbl <- data.frame(xtable(sloglinmod1))
tbl

Persamaan log-linear upah memberikan contoh lain untuk menghitung tingkat pertumbuhan, tetapi kali ini variabel bebasnya bukan waktu , tetapi pendidikan . Prediksi dan kemiringan dihitung untuk pendidikan=12 tahun.

data("cps4_small", package="PoEdata")
xeduc <- 12
loglingmod2 <- lm(log(wage)~educ, data=cps4_small)
sloglinmod2 <- summary(loglingmod2)
tabl <- data.frame(xtable(sloglinmod2))
tabl

b1 <- coef(sloglinmod2)[[1]]
b2 <- coef(sloglinmod2)[[2]]
sighat2 <- sloglinmod2$sigma^2
g <- 100*b2
yhatn <- exp(b1+b2*xeduc)
yhatc <- exp(b1+b2*xeduc+sighat2/2)
DyDx <- b2*yhatn
print(c(g, yhatn, yhatc, DyDx))

Berikut adalah hasil dari perhitungan tersebut: prediksi “alami” \(\hat{y}_n=14.796\); prediksi terkoreksi, \(\hat{y}_c=16.996\); tingkat pertumbuhan \(g=9.041\); dan efek marginal \(\frac{dy}{dx}=1.34\). Laju pertumbuhan tersebut menunjukkan bahwa peningkatan pendidikan sebesar satu satuan (lihat uraian data menggunakan ?cps4_small) meningkatkan upah per jam sebesar 9,041 persen.
Gambar di bawah menyajikan garis regresi “alami” dan “dikoreksi” untuk persamaan upah, bersama dengan titik data aktual.

education=seq(0,22,2)
yn <- exp(b1+b2*education)
yc <- exp(b1+b2*education+sighat2/2)
plot(cps4_small$educ, cps4_small$wage, xlab="education", ylab="wage", col="yellow")
lines(yn~education, lty=2, col="black")
lines(yc~education, lty=1, col="blue")
legend("topleft", legend=c("yc","yn"), lty=c(1,2), col=c("blue","black"))

result

quadmod3 <- lm(wage~I(educ^2), data=cps4_small)
yhat3 <- predict(quadmod3)
loglinmod3 <- lm(log(wage)~educ, data=cps4_small)
sloglinmod3 <- summary(loglinmod3)
b1 <- coef(sloglinmod3)[[1]]
b2 <- coef(sloglinmod3)[[2]]
sighat2 <- sloglinmod3$sigma^2
yhat4 <- exp(b1+b2*cps4_small$educ+sighat2/2)
rquadmod3 <- cor(cps4_small$wage, yhat3)^2
rgloglinmod3 <- cor(cps4_small$wage, yhat4)^2
print(c(rquadmod3, rgloglinmod3))

Catatan: Hasil rquadmod3 adalah 0.1881762 dan hasil rloglinmod3 adalah 0.1859307.
Model kuadrat menghasilkan \(R^2=0,188\), dan model log-linier menghasilkan \(R^2=0,186\); karena yang pertama lebih tinggi, kami menyimpulkan bahwa model kuadrat lebih cocok dengan data daripada yang log-linear. (Namun, pengujian lain tentang bagaimana kedua model memenuhi asumsi regresi linier dapat mencapai kesimpulan yang berbeda; \(R^2\) hanyalah salah satu kriteria pemilihan model.)
Untuk menentukan perkiraan interval perkiraan dalam model log-linier, pertama-tama kita membangun interval dalam log menggunakan prediktor alami y^n, kemudian mengambil antilog dari batas interval. Perhitungan berikut menggunakan tingkat pendidikan sama dengan 12 dan = 0,05.

alpha <- 0.05
xeduc <- 12
xedbar <- mean(cps4_small$educ)
loglinmod4 <- lm(log(wage)~educ, data=cps4_small)
b1 <- coef(loglinmod4)[[1]]
b2 <- coef(loglinmod4)[[2]]
df5 <- loglinmod4$df.residual
N <- nobs(loglinmod4)
tcr <- qt(1-alpha/2, df=df5)
sloglinmod4 <- summary(loglinmod4)
varb2 <- vcov(loglinmod4)[2,2]
sighat2 <- sloglinmod4$sigma^2
varf <- sighat2+sighat2/N+(xeduc-xedbar)^2*varb2
sef <- sqrt(varf)
lnyhat <- b1+b2*xeduc
lowb <- exp(lnyhat-tcr*sef)
upb <- exp(lnyhat+tcr*sef)
print(c(lowb, upb))

Hasilnya adalah selang kepercayaan (5.26,41.62) . Selanjutnya, mari kita pertimbangkan garis kepercayaan 95% untuk model upah log-linear di bawah ini:

xmin <- min(cps4_small$educ)
xmax <- max(cps4_small$educ)+2
education <- seq(xmin, xmax, 2)
lnyhat <- b1+b2*education
yhat <- exp(lnyhat)
varf <- sighat2+sighat2/N+(education-xedbar)^2*varb2
sef <- sqrt(varf)
lowb <- exp(lnyhat+tcr*sef)
upb <- exp(lnyhat+tcr*sef)
plot(cps4_small$educ, cps4_small$wage, col="grey",
     xlab="education", lylab="wage", ylim=c(0,100))
lines(yhat~education, lty=1, col="black")
lines(lowb~education, lty=2, col="blue")
lines(upb~education, lty=2, col="blue")
legend("topleft", legend=c("yhat", "lowb", "upb"),
       lty=c(1, 2, 2), col=c("black", "blue", "blue"))

result

Model Log-log

Model log-log memiliki sifat yang diinginkan yaitu koefisien variabel bebas sama dengan elastisitas (konstan) y terhadap x. Oleh karena itu, model ini sering digunakan untuk mengestimasi persamaan penawaran dan permintaan. Bentuk standarnya diberikan dalam persamaan 15, di mana y,x, dan e adalah vektor \(N×1\).

\[\begin{align*} log(y)=\beta_1+\beta_2log(x)+\epsilon_i \end{align*}\]

data("newbroiler", package="PoEdata")
loglogmod <- lm(log(q)~log(p), data=newbroiler)
b1 <- coef(loglogmod)[[1]]
b2 <- coef(loglogmod)[[2]]
sloglogmod <- summary(loglogmod)
tbl <- data.frame(xtable(sloglogmod))
tbl

Tabel di atas memberikan output regresi log-log. Koefisien pada p menunjukkan bahwa kenaikan harga sebesar 1% mengubah jumlah yang diminta sebesar 1.121% .

ngird <- 20
xmin <- min(newbroiler$p)
xmax <- max(newbroiler$p)
step <- (xmax-xmin)/ngird
xp <- seq(xmin, xmax, step)
sighat2 <- sloglogmod$sigma^2
yhatc <- exp(b1+b2*log(newbroiler$p)+sighat2/2)
yc <- exp(b1+b2*log(xp)+sighat2/2)
plot(newbroiler$p, newbroiler$q, ylim=c(10,60), xlab="price", ylab="quantity")
lines(yc~xp, lty=1, col="black")

result

Untuk R-squared:

rgsq <- cor(newbroiler$q, yhatc)^2
rgsq

\(R^2\) umum, yang menggunakan nilai pas yang dikoreksi, sama dengan 0,8818.