Pengenalan

Regresi linier sederhana digunakan untuk memprediksi hasil kuantitatif y berdasarkan satu variabel prediktor tunggal x. Tujuannya adalah untuk membangun model matematika (atau rumus) yang mendefinisikan y sebagai fungsi dari variabel x. Setelah kita membangun model yang signifikan secara statistik, model ini dapat digunakan untuk memprediksi hasil di masa mendatang berdasarkan nilai x baru.

Untuk memeriksa asumsi regresi, kami akan memeriksa distribusi residual. Regresi linier membuat beberapa asumsi tentang data, seperti:

• Linearitas data. Hubungan antara prediktor \((x)\) dan hasil \((y)\) diasumsikan linier.
• Normalitas residu. Residual error diasumsikan berdistribusi normal.
• Homogenitas varians residual. Residu diasumsikan memiliki varians yang konstan (homoskedastisitas)
• Independensi kesalahan residual.

Anda harus memeriksa apakah asumsi ini benar atau tidak. Potensi masalah meliputi:

• Non-linearitas hasil - hubungan prediktor
• Heteroskedastisitas: Varians istilah kesalahan yang tidak konstan.
• Adanya nilai-nilai berpengaruh dalam data yang dapat berupa: ** Pencilan: nilai ekstrim dalam variabel hasil \((y)\) ** Poin leverage tinggi: nilai ekstrem dalam variabel prediktor \((x)\)

Catatan: Semua asumsi dan masalah potensial ini dapat diperiksa dengan membuat beberapa plot diagnostik yang memvisualisasikan residual error.

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

Prediksi vs Meramal

Dengan asumsi bahwa nilai ekspektasi dari ketentuan error dalam Persamaan 1 adalah nol, Persamaan 2 memberikan \(\hat{y}_i\), nilai prediksi dari ekspektasi \(y_i\) yang diberikan \(x_i\), di mana \(b1\) dan \(b2\) adalah estimasi (kuadrat terkecil) dari parameter regresi \(\beta_1\) dan \(\beta_2\).

\[\begin{align*} \hat{y}_i = b_1 + b_2xi \end{align*}\] \((2)\)

Nilai prediksi \(\hat{y}_i\) adalah variabel acak karena bergantung pada sampel; oleh karena itu, kita dapat menghitung interval kepercayaan dan menguji hipotesisnya, kemudian kita dapat menentukan distribusi dan variansnya. Prediksi yang memiliki distribusi normal, menjadi kombinasi linier dari dua variabel acak terdistribusi normal \(b1\) dan \(b2\), dan variansnya diberikan oleh 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-s\) dalam kumpulan data.

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

yang dapat direduksi menjadi

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

Catatan: Varians pada Persamaan 3 tidak sama dengan varians pada Persamaan 4; yang pertama adalah varians dari perkiraan ekspektasi \(y\), sedangkan yang terakhir adalah varians dari kejadian tertentu \(y\). Mari kita sebut yang terakhir varians dari kesalahan perkiraan. Tidak mengherankan, varians kesalahan ramalan lebih besar dari varians dari prediksi \(E(y|x)\).

Mari kita tentukan kesalahan standar untuk persamaan makanan untuk pendapatan rumah tangga (income) $2000 seminggu, yaitu, pada \(x=x0=20\), menggunakan Persamaan 4 ; untuk melakukannya, kita perlu mengambil \(var(b2)\) dan \(\sigma^2\), kesalahan standar regresi dari hasil regresi.

library('PoEdata')                      # load library PoEdata
data("food")                            # memanggil data set dari PoEdata
alpha<-0.05                             # masukan nilai signifikan
x<-20                                   # masukan nilai baru untuk income
xbar<-mean(food$income)                 # rata-rata dari income
linmod<-lm(food_exp~income, data=food)  # linear model untuk income
b1<-coef(linmod)[[1]]                   # untuk mengestimasi \beta 1
b2<-coef(linmod)[[2]]                   # untuk mengestimasi \beta 2
yhatx<-b1+b2*x                          # untuk mengestimasi regresi linier sederhana
slinmod<-summary(linmod)                # summary dari linier sederhana
df<-df.residual(linmod)                 # residual
tcr<-qt(1-alpha/2, df)                  # nilai kritis t
N<-nobs(linmod)                         # angka dari observasi, N
N<-NROW(food)                           # cara lain untuk menemukan N
varb2<-vcov(linmod)[2,2]                # menghitung kofarians
sighat2<-slinmod$sigma^2                # estimasi variansi
varf<-sighat2+sighat2/N+(x-xbar)^2*varb2 # memperkirakan variasnsi
sef<-sqrt(varf)                         # standar error dari perkiraan
lb<-yhatx-tcr*sef                       # batas bawah
ub<-yhatx+tcr*sef                       # batas atas
print(c(lb,ub))                         # print hasil

## [1] 104.1323 471.0854

Hasilnya adalah interval kepercayaan untuk forecast [104.1322769, 471.0854459], 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)

##        fit      lwr      upr
## 1 287.6089 258.9069 316.3108

predict(linmod, newdata=incomex, interval="prediction",level=0.95)

##        fit      lwr      upr
## 1 287.6089 104.1323 471.0854

Sekarang mari kita gunakan fungsi predict() untuk mereplikasi pita estimasi interval, titik-titik dalam kumpulan data. (kita akan menciptakan nilai baru untuk pendapatan hanya untuk tujuan merencanakan.)

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("black", "red", "red", "blue", "blue"),
        ylab="food expenditure", xlab="income")
points(food$income, food$food_exp)
legend("topleft", 
        legend=c("E[y|x]", "lwr_pred", "upr_pred","lwr_forcst","upr_forcst"),
        lty=c(1, 2, 2, 3, 3), 
        col=c("black", "red", "red", "blue", "blue"))

Gambar di atas menyajikan pita yang diprediksi dan yang diperkirakan pada grafik yang sama, untuk menunjukkan bahwa mereka memiliki perkiraan titik yang sama (garis hitam, solid) dan bahwa pita yang diperkirakan jauh lebih besar daripada yang diprediksi. Dengan kata lain, Anda mungkin berpikir tentang perbedaan antara dua jenis interval yang kami sebut prediksi dan ramalan sebagai berikut:

• Interval prediksi tidak seharusnya mencakup, katakanlah, 95 persen 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.

Koefisien Determinasi

Tujuan penting dari analisis data adalah untuk menjelaskan variabilitas dalam variabel dependen, \(y\). Untuk membandingkan model, kita mempertimbangkan situasi dasar di mana kami 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 rata-rata 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 rata-rata dari \(y\), Sum of Squares due to Segression (SSR), seperti yang ditunjukkan oleh Persamaan 5.

\[\begin{align*} SST=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[2,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('Tabel 1: Output yang dihasilkan oleh fungsi anova.')),
        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          # print summary dari model regresi linmod

## 
## Call:
## lm(formula = food_exp ~ income, data = food)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -223.025  -50.816   -6.324   67.879  212.044 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   83.416     43.410   1.922   0.0622 .  
## income        10.210      2.093   4.877 1.95e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 89.52 on 38 degrees of freedom
## Multiple R-squared:  0.385,  Adjusted R-squared:  0.3688 
## F-statistic: 23.79 on 1 and 38 DF,  p-value: 1.946e-05

(rsq<-slinmod$adj.r.squared)      # atau print hasil dari adj.r.squared

## [1] 0.3688181

Catatan

• Dalam model regresi sederhana kita, jumlah kuadrat karena regresi hanya mencakup pendapatan variabel. Dalam contoh makanan kita, \(R^2=0,3688\). yang menunjukkan model kami tidak cocok dengan data yang diamati. Untuk memastikan apakah ada hubungan linier yang signifikan secara statistik antara variabel dependen dan independen, pada bagian selanjutnya kita akan mempertimbangkan uji-T: uji signifikansi kemiringan.

• Dalam model regresi berganda, yaitu model dengan lebih dari satu variabel bebas, jumlah kuadrat akibat regresi sama dengan jumlah kuadrat semua variabel bebas, kita akan menggunakan Overall F Test: evaluasi model global.

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 H_a : \beta_1 \neq 0 \end{align*}\] \((7)\)

Statistik uji: \(t=\frac{b_1}{se_b_1}\), di mana seb1 adalah kesalahan standar b1. Tes 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.

Untuk melakukannya, kita bisa mendapatkan nilai p yang sesuai dengan nilai t yang dihitung, kita bisa merujuk ke situs web bernama ‘nilai p dari kalkulator skor T atau kalkulator skor F. Dalam hal ini, saya memberikan ’skor T’ sebagai 4,877, dan ‘DF’ sebagai 38, yang merupakan df dari MSE, n-2, dan akhirnya kita mendapatkan nilai p sebagai 1,946e-05, yang merupakan nilai p yang sama dari F-test di atas.

Jika Anda ingin menemukan nilai-p di R, jalankan kode ini di konsol Anda:

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

## [1] 1.945862e-05

Dari hasil tersebut diketahui bahwa \(p\)-value = 1.945862e-05 lebih kecil dari taraf signifikansi \(\alpha=0.05\). Oleh karena itu, kami menolak hipotesis nol bahwa kemiringannya nol dan menyimpulkan bahwa ada hubungan linier yang signifikan antara \(y\) dan \(x\).

Nilai dan Residual yang Dipasang

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 Anda lebih banyak tentang data Anda. Di R, Anda dapat dengan mudah menambah data Anda 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 distriburion 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)

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)

Di bagian berikut, kami akan menjelaskan, secara rinci, cara menggunakan grafik dan metrik ini untuk memeriksa asumsi regresi dan untuk mendiagnosis masalah potensial dalam model.

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 = "grey")

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

jarque.bera.test(ehat)

## 
##  Jarque Bera Test
## 
## data:  ehat
## X-squared = 0.06334, df = 2, p-value = 0.9688

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

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)

Quadratic Models

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 \(income^2\) yang merupakan kuadrat dari pendapatan variabel.

\[\begin{align*} food_{exp} =b_1 + b_2 x_{i}^2+e_{i} \end{align*}\] \((8)\)

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

## 
## Call:
## lm(formula = food_exp ~ income2, data = food)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -219.142  -56.559   -6.527   72.755  198.143 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 166.45749   29.22934   5.695 1.49e-06 ***
## income2       0.27232    0.05907   4.610 4.46e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 91.41 on 38 degrees of freedom
## Multiple R-squared:  0.3587, Adjusted R-squared:  0.3418 
## F-statistic: 21.25 on 1 and 38 DF,  p-value: 4.455e-05

Catatan: sintaks yang terlibat dalam pemasangan model linier dengan dua atau lebih prediktor. Kami menyertakan setiap prediktor dan memberi tanda plus di antara mereka.

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. Misalkan kita tertarik untuk memahami hubungan antara jumlah jam kerja dan kebahagiaan yang dilaporkan. Kami memiliki data berikut tentang jumlah jam kerja per minggu dan tingkat kebahagiaan yang dilaporkan (pada skala 0-100) untuk 11 orang yang berbeda:

# buat data
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))

# buat variabel baru untuk hours2
data$hours2 <- data$hours^2

# sesuaikan kuadrat model regresi
quadmod2<-lm(happiness~hours+hours^2, data=data)
summary(quadmod2)

## 
## Call:
## lm(formula = happiness ~ hours + hours^2, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -39.34 -21.99  -2.03  23.50  35.11 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  48.4531    17.3288   2.796   0.0208 *
## hours         0.2981     0.4599   0.648   0.5331  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.72 on 9 degrees of freedom
## Multiple R-squared:  0.0446, Adjusted R-squared:  -0.06156 
## F-statistic: 0.4201 on 1 and 9 DF,  p-value: 0.5331

Polinomial Model

Mari kita coba menerapkan model polinomial ke kumpulan data food:

\[\begin{align*} food_{exp} =b_1 + b_2 x_{i}^3+e_{i} \end{align*}\] \((9)\)

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

## 
## Call:
## lm(formula = food_exp ~ income + I(income2), data = food)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -223.037  -49.992   -5.946   67.605  212.458 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 81.982928  73.089590   1.122    0.269
## income      10.405746   8.263941   1.259    0.216
## I(income2)  -0.005607   0.228375  -0.025    0.981
## 
## Residual standard error: 90.72 on 37 degrees of freedom
## Multiple R-squared:  0.385,  Adjusted R-squared:  0.3518 
## F-statistic: 11.58 on 2 and 37 DF,  p-value: 0.0001242

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 (makanan).

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. Model linier diberikan

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

## 
## Call:
## lm(formula = greenough ~ time + I(time^2), data = wa_wheat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47116 -0.15002  0.00728  0.16748  0.31939 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.719e-01  5.069e-02  17.200   <2e-16 ***
## time        1.572e-04  1.339e-04   1.173   0.2468    
## I(time^2)   1.402e-07  6.347e-08   2.209   0.0323 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1848 on 45 degrees of freedom
## Multiple R-squared:  0.7551, Adjusted R-squared:  0.7442 
## F-statistic: 69.39 on 2 and 45 DF,  p-value: 1.782e-14

Linear-Log Model

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. Dalam model pengeluaran makanan, misalnya, masuk akal untuk meyakini bahwa jumlah yang dibelanjakan untuk makanan meningkat lebih cepat pada pendapatan yang lebih rendah daripada pada pendapatan yang lebih tinggi. 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})+e_{i} \end{align*}\] \((10)\)

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*}\] \((11)\)

Terkait dengan model log linier, ukuran lain yang menarik dalam ekonomi adalah semi-elastisitas \(y\) terhadap \(x\), yang diberikan oleh Persamaan 12 . Semi-elastisitas menunjukkan bahwa perubahan \(x\) sebesar \(1%\) mengubah \(y\) sebesar \(\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*}\] \((12)\)

Besaran lain yang mungkin menarik adalah elastisitas \(y\) terhadap \(x\), yang diberikan oleh Persamaan 13 dan menunjukkan bahwa satu persen peningkatan \(x\) menghasilkan \((β_2/y)\) persen perubahan \(y\).

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

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))
datatable(tbl,
          caption = htmltools::tags$caption(
            style = 'caption-side: bottom; text-align: center;',
            htmltools::em('Table 3: Linear-Log Model Output for the food example.')),
          options=list(dom = 't'))

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 Expendicture")
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")

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

## [1] 13.216584  1.321658  0.638061

Hasil untuk pendapatan sebesar $1000 adalah sebagai berikut: \(dy/dx=13.217\), yang menunjukkan bahwa peningkatan pendapatan sebesar $100 (yaitu, satu unit x ) meningkatkan pengeluaran sebesar $13.217 ; untuk peningkatan pendapatan sebesar 1%, yaitu peningkatan sebesar $10, pengeluaran meningkat sebesar $1,322 ; dan, akhirnya, untuk peningkatan 1% dalam pengeluaran pendapatan meningkat sebesar 0,638%.

Log-Linear Model

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. Model log-linier umum diberikan dalam Persamaan 14.

\[\begin{align*} log(y_{i}) = \beta_{1}+ \beta_{2}x_{i}+e_{i} \end{align*}\] \((14)\)

Rumus berikut dengan mudah diturunkan dari persamaan log-linier 14. 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 100\(b_2\) persen dalam \(y\). Untuk perubahan kecil pada \(x\), jumlah 100\(b_2\) dalam model log-linier juga dapat diinterpretasikan sebagai laju pertumbuhan dalam \(y\) (sesuai dengan peningkatan unit dalam \(x\)). Misalnya, jika \(x\) adalah waktu, maka 100\(b_2\) adalah laju pertumbuhan dalam y per satuan waktu.

• Prediksi \(\hat{y}_{n}=exp(b_{1}+b_{2}x)\), atau \(\hat{y}_{c}=exp(b_{1}+b_{2}x+\frac{\hat{\sigma}^2}{2}\) dengan prediktor “natural” \(\hat{y}_n\) untuk digunakan dalam sampel kecil dan prediktor “dikoreksi”, \(\hat{y}_c\), dalam sampel besar

• Marginal effect (slope): \(\frac{dy}{dx}=b_{2}y\)

• Semi-elasticity: \(\%\Delta y = 100b_{2}\Delta x\)

Mari kita lakukan perhitungan ini terlebih dahulu untuk persamaan hasil menggunakan dataset wa_wheat.

library(xtable)
library(kableExtra)
loglinmod1 <- lm(log(greenough)~time, data=wa_wheat)
sloglinmod1 <- summary(loglinmod1)
tbl <- data.frame (xtable(sloglinmod1))
datatable(tbl, caption = htmltools::tags$caption(
  style = 'caption-side: botton; text-align: center;',
  htmltools::em("Table 4: Log-linear model for the yield equation.")),
  options = list(dom = 't'))

Tabel 4 memberikan \(b_2=0,017844\) , yang menunjukkan bahwa laju pertumbuhan produksi gandum meningkat rata-rata sekitar 1,78 persen per tahun.

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 educ=12 tahun.

data("cps4_small", package="PoEdata")
xeduc = 12
loglinmod2 <- lm(log(wage)~educ, data=cps4_small)
sloglinmod2 <- summary(loglinmod2)
tabl <- data.frame(xtable(sloglinmod2))
datatable(tabl,
          caption = htmltools::tags$caption( 
            style = 'caption-side: bottom, text-align: center;',
            htmltools::em('Table 5: Log-linear wage regression output.')), 
          options = list(dom = 't'))

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)) 

## [1]  9.040825 14.795801 16.996428  1.337662

Berikut adalah hasil dari perhitungan tersebut: prediksi “natural” \(\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 “natural” 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="grey") 
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"))

Catatan: \(R^2\) reguler tidak dapat digunakan untuk membandingkan dua model regresi yang memiliki variabel dependen yang berbeda seperti model log-linier dan log-linier; ketika perbandingan seperti itu diperlukan, seseorang dapat menggunakan \(R^2\) umum, yaitu \(R_{g}^2=[corr(y,\hat{y}]^2\). Mari kita hitung \(R^2\) umum untuk model upah kuadratik dan log-linier.

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) 
rgquadmod3 <- cor(cps4_small$wage, yhat3)^2
rgloglinmod3 <- cor(cps4_small$wage, yhat4)^2
print(c(rgquadmod3, rgloglinmod3))

## [1] 0.1881762 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 \(\hat{y}_n\), kemudian mengambil antilog dari batas interval. Kesalahan peramalan sama seperti sebelumnya, diberikan dalam Persamaan 4. Perhitungan berikut menggunakan tingkat pendidikan sama dengan 12 dan \(\alpha= 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))

## [1]  5.260398 41.615811

Hasilnya adalah selang kepercayaan (5.26,41.62) . Selanjutnya, mari kita pertimbangkan band kepercayaan 95% untuk model wage 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", ylab="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"))

The Log-Log Model

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\).

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)) 
datatable (tbl, 
           caption=htmltools::tags$caption(
           style='caption-side: bottom; text-align: center;',
           htmltools::em ("Table 6: The log-log poultry regression equation.")),
           options=list(dom='t'))

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

ngrid <-20 
xmin <- min(newbroiler$p)
xmax <-max(newbroiler$p)
step <- (xmax-xmin)/ngrid 
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")

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

Refrensi

https://rpubs.com/dsciencelabs/econometrics4