Email         : naufal3433@gmail.com
RPubs       : https://rpubs.com/muhammad_naufal/
Jurusan : Statistika
Address : Jalan Gunung Galunggung 5 Blok E9, No. 21
   .

1 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 kami 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. Kesalahan residual diasumsikan terdistribusi normal.
Homogenitas varians residual. Residu diasumsikan memiliki varians yang konstan (homoskedastisitas)
Independensi istilah kesalahan residual.

Selain itu, diperlukan pengecekan kembali terhadap asumsi apakah asumsi ini benar atau tidak melalui beberapa potensi masalah seperti:

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)) dan Poin leverage tinggi (nilai ekstrim dalam variabel prediktor (x)).

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

$Y_i= β_1+β_2x_i+e_i$

2 Prediksi VS Peramalan

Dengan asumsi bahwa nilai ekspektasi dari error term dalam Persamaan 1 adalah nol, Persamaan 2 memberikan y^i, nilai prediksi dari ekspektasi yi yang diberikan xi, di mana b1 dan b2 adalah estimasi (kuadrat terkecil) dari parameter regresi 1 dan 2.

$y^i=b_1+b_2x_i$

Nilai prediksi y^i adalah variabel acak karena bergantung pada sampel; oleh karena itu, kita dapat menghitung interval kepercayaan dan menguji hipotesis tentangnya, kemudian kita dapat menentukan distribusi dan variansnya. Prediksi memiliki distribusi normal, menjadi kombinasi linier dari dua variabel acak terdistribusi normal b_1 dan b_2, 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.

$var(fi)ˆ=σ^2[1+[1/N]+(x_i−x¯)^2/∑N_j=1(x_j−x¯)^2]$

yang mana dapat bereduksi menjadi

$var(fi)ˆ=σ^2+[σ^2/N]+(x_i−x¯)^2∗var(b2)ˆ$

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 (pendapatan) $2000 seminggu, yaitu, pada x=x0=20, menggunakan Persamaan 4 ; untuk melakukannya, kita perlu mengambil var(b2) dan σ, kesalahan standar regresi dari keluaran regresi..

library('remotes')
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)[1,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))

## [1] 104.2899 470.9279

Hasilnya adalah interval kepercayaan untuk ramalan [104.1322769, 471.0854459], yang, seperti yang diharapkan, lebih besar dari interval kepercayaan dari perkiraan nilai y yang diharapkan berdasarkan Persamaan 4 dan 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="1")
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)

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

xmax<-max(food$income)
xmin<-min(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="prediction")
matplot(income, cbind(ypredict[,1], ypredict[,2], ypredict[,3], yforecast[,2], yforecast[,3]),
        type="1",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.

3 Kofisien Determinasi

Tujuan penting dari analisis data adalah untuk menjelaskan variabilitas dalam variabel dependen, y. Untuk membandingkan model, kami mempertimbangkan situasi dasar di mana kami memiliki informasi tentang variabel dependen saja. Dalam situasi ini, kita perlu memodelkan mean dari y, 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, y. Diharapkan, kita dapat mengharapkan pengurangan variabilitas x 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.

$SST=SSR+SSE$

Jika membutuhkan jumlah kesalahan kuadrat, SSE, atau jumlah kuadrat karena regresi, SSR, gunakan fungsi anova yang strukturnya ditunjukkan pada Tabel dibawah ini.

library(stats)
#install.packages("hms")
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('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 σ^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.

$R^2=(SSR/SST)=1−[SSE/SST]$

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

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

## [1] 0.3688181

catatan

Dalam model regresi sederhana ini jumlah kuadrat karena regresi hanya mencakup pendapatan variabel. Dalam contoh makanan kita, R2=0,3688. yang menunjukkan model tersebut tidak cocok dengan data yang diamati. Untuk memastikan apakah ada hubungan linier yang signifikan secara statistik antara variabel dependen dan independen maka akan dilakukan dengan 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 uji F Keseluruhan: evaluasi model global.

4 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:

$H_0:β1=0$ $H_a:β1≠0.$

Statistik uji: b1/SE_(b1)\,di mana SE_(b1) 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.

Dengan begitu akan mendapatkan nilai p yang sesuai dengan nilai t yang dihitung dan merujuk ke situs web bernama ‘nilai p dari T score calculator atau F score calculator . Dalam hal ini,diberikan ’skor T’ sebagai 4,877, dan ‘DF’ sebagai 38, yang merupakan df dari MSE, n-2, dan akhirnya mendapatkan nilai p sebagai 1,946e-05, yang merupakan nilai p yang sama dari F-test di atas.

Jika ingin menemukan nilai-p di R, jalankan kode dibawah ini.

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

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

5 Fitted dan Nilai Residual

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, diperlukannya mendiagnosa model regresi yang akan bangun untuk mendeteksi potensi masalah dan untuk memeriksa apakah asumsi yang dibuat oleh model regresi linier terpenuhi atau tidak.

Dengan begitu dilakukan pemeriksaan distribusi kesalahan residual, yang dapat memberi tahu lebih banyak tentang data yang akan dianalisis. Di R 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.

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

Setelah meilhat data diatas maka 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)

## `geom_smooth()` using formula 'y ~ x'

6 Diagnosis Regresi

Diagnostik regresi adalah salah satu dari serangkaian prosedur yang tersedia untuk analisis regresi yang berusaha menilai validitas model dengan salah satu dari sejumlah 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)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:kableExtra':
## 
##     group_rows

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

#install.packages("ggfortify")
library(ggfortify)
par(mfrow=c(2,2))
autoplot(linmod)

Di bagian berikut akan menjelaskan secara rinci cara menggunakan grafik dan metrik ini untuk memeriksa asumsi regresi dan untuk mendiagnosis masalah potensial di masa depan.

6.1 Nilai 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. Dalam contoh kita, tidak ada pola di plot residual. Hal ini menunjukkan bahwa dapat mengasumsikan hubungan linier antara prediktor dan variabel hasil.

Perhatikan bahwa, jika plot residual menunjukkan hubungan non-linier dalam data, maka pendekatan sederhana adalah dengan menggunakan transformasi non-linier dari prediktor, seperti (x)\, √x\ dan ^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="grey")

6.2 QQ Normal

Plot QQ residual dapat digunakan untuk memeriksa asumsi normalitas secara visual. Plot probabilitas normal dari residual kira-kira harus mengikuti garis lurus. Pada contoh ini, semua titik jatuh kira-kira di sepanjang garis referensi ini, sehingga dapat mengasumsikan normalitas.

Asumsi lain yang ingin di 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)

## Loading required package: xts

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## Attaching package: 'xts'

## The following objects are masked from 'package:dplyr':
## 
##     first, last

## Loading required package: TTR

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

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

Sementara histogram di masa depan di atas mungkin tidak mendukung satu kesimpulan atau lainnya tentang normalitas ehat, sehingga dapat melakukan uji Jarque-Bera, yang hipotesis nolnya adalah “Deret terdistribusi normal”. Jika p-value < aplha, 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)

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

6.3 Skala - Laokasi

Plot ini menunjukkan jika residual tersebar merata di sepanjang rentang prediktor. Ada baiknya jika Anda melihat garis horizontal dengan titik penyebaran yang sama. Dalam contoh 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).

Note: Solusi yang mungkin untuk mengurangi masalah heteroskedastisitas adalah dengan menggunakan transformasi log atau akar kuadrat dari variabel hasil x.

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

6.4 Residual VS Laverage

Pencilan / Outliners: Pencilan adalah titik yang memiliki nilai variabel hasil ekstrem. Kehadiran outlier dapat mempengaruhi interpretasi model karena meningkatkan RSE. Outliner 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 \2(p+1)/n\ menunjukkan observasi dengan leverage yang tinggi (P. Bruce dan Bruce 2017); dimana, \ adalah jumlah prediktor dan \ adalah jumlah observasi.

Dalam contoh ini, 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).

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

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.

7 Model Kuadratik

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

$food_(exp)=b_1+b_2x^2_i+e_i$

#install.packages("remotes")
remotes::install_github("ccolonescu/PoEdata")

## Skipping install of 'PoEdata' from a github remote, the SHA1 (c41bee8f) has not changed since last install.
##   Use `force = TRUE` to force installation

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. Maka diatas disertakan setiap prediktor dan memberi tanda plus di antara mereka.

Hasil dari model kuadratik di atas ternyata tidak sesuai dengan yang diharapkan. Maka akan dilihat ketika memiliki contoh yang baik untuk model kuadrat. Misalkan tertarik untuk memahami hubungan antara jumlah jam kerja dan kebahagiaan yang dilaporkan. Maka disini menggunakan data berikut tentang jumlah jam kerja per minggu dan tingkat kebahagiaan yang dilaporkan (pada skala 0-100) untuk 11 orang yang berbeda.

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

8 Model Polimonial

$food_(exp)=b_1+b_2x^3_i+e_i$

Berikut adalah contoh penerapan model polinomial pada data food.

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

9 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. 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.

$y_i=β_1+β_2∗log(x_i)+e_i$

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.

$dy/dx=β_2/x$

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 sebesar β_2/100%\ unit y. Karena semi-elastisitas juga berubah ketika x berubah, itu seharusnya hanya ditentukan untuk perubahan kecil pada x.

$dy=β_2/100[%Δx]$/>

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

$dy=β_2/100[%Δx]$

Berikut adalah 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] 1.02096430 0.10209643 0.09548453

Hasil untuk pendapatan sebesar $1000\ adalah sebagai berikut: /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%\.

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

$log(y_i)=β_1+β_2x_i+e_i$

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 satu unit (dari x) menghasilkan perubahan \100b2\ persen dalam x. Untuk perubahan kecil pada x, jumlah \100b2\ dalam model log-linier juga dapat diinterpretasikan sebagai laju pertumbuhan dalam x (sesuai dengan peningkatan unit dalam x). Misalnya, jika x adalah waktu, maka \100b2\ adalah laju pertumbuhan dalam x per satuan waktu.

Prediksi : $y_n=exp(b_1+b_2x)$, or $y_c=exp(b_1+b_2x+σ^2/2)$, dengan prediktor “alami” _n\ untuk digunakan dalam sampel kecil dan prediktor “dikoreksi”, _c\, dalam sampel besar
Slope / efek marjinal :$dy/dx=b_2y$
Semi-elastisitas :$%Δy=100b_2 Δx$

Let us do these calculations first for the yield equation using the wa_wheat dataset.

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 telah meningkat dengan laju 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 “alami” ^n=14.796\; prediksi terkoreksi, _c=16.996\; tingkat pertumbuhan \; dan efek marginal /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="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^{2_g=[corr(y,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 dengan membangun interval dalam log menggunakan prediktor alami 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 α = 0,05.

# The wage Log-Linear model 
#Prediction interval for educ= 12

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, dengan pertimbangkan pita kepercayaan 95% untuk model upah log-linear di bawah ini:

#Drawing a confidence band for the Log-Linear
#wage equation

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

11 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 x,x, dan e adalah vektor N×1.

$log(y)=β_1+β_2log(x)+e_i$

# Calculating Log-Log demand for chicken
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% .

#Drawing the fitted values of the Log-Log equation

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

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

## [1] 0.8817758

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

12 Referensi

[1] Bruce, Peter, and Andrew Bruce. 2017. Practical Statistics for Data Scientists. O’Reilly Media.

[2] James, Gareth, Daniela Witten, Trevor Hastie, and Robert Tibshirani. 2014. An Introduction to Statistical Learning: With Applications in R. Springer Publishing Company, Incorporated.

EKONOMETRIK

Tugas 4

Muhammad Naufal Ardiansyah (20204920017)

February 22, 2022