Di sesi ini kita akan mempelajari bagaimana kita bisa menggunakan R untuk melakukan analisis regresi linear. Di akhir sesi, kita akan belajar juga beberapa metode untuk melaporkan hasil regresi dengan membawa hasil di dalam R ke dokumen eksternal.
Catatan: Sumber utama tutorial adalah buku Danielle Navarro: Learning Statistics with R (https://learningstatisticswithr.com/) dan buku Field, et. al: Discovering Statistics using R (https://www.discoveringstatistics.com/books/discovering-statistics-using-r/).
Untuk contoh, kita akan menggunakan dataset buatan mengenai usia, tinggi badan, berat badan, asupan protein per hari, dan catatan apakah tiap individu aktif berolahraga atau tidak.
head(df)Regresi adalah salah satu fungsi yang ada di dalam base R. Fungsi ini, lm(), pada dasarnya dapat digunakan paling tidak dengan menulis dua input, yaitu formula dan data.
Di dalam formula, kita masukkan rumus model regresi yang ingin kita uji dengan rumus kriterion ~ prediktor, di mana dalam contoh di bawah ini kriterionnnya adalah berat badan, weight, dan prediktornya adalah tinggi badan, height. Lalu kita tinggal menuliskan dataset yang digunakan.
model_wh <- lm(formula = weight ~ height,
data = df)Perhatikan bahwa di dalam perintah di atas, model ini kita simpan sebagai sebuah objek, model_wh. Kita bisa memanggil model_wh untuk mendapatkan intercept dan koefisien regresi dari variabel yang kita gunakan di dalam model.
model_whNamun, mengetahui hanya koefisien dan intercept dari model tidak cukup. Seringkali kita ingin mendapatkan informasi mengenai signifikansi dari masing-masing model, degrees of freedom, serta R-squared dan signifikansi dari model secara umum. Untuk itu, kita bisa menggunakan summary() dan memanggil model regresi tersebut.
summary(model_wh)Dari summary() kita bisa melihat bahwa height secara signifikan memprediksi weight, dengan P-value yang sangat rendah. Kita juga bisa melihat bahwa F-statistic dari model signifikan, indikasi bahwa model ini benar-benar dapat menjelaskan data, tidak sekadar kebetulan.
Model di atas cukup sederhana, satu kriterion dan satu prediktor. Ketika kita ingin melihat efek dari lebih dari dua prediktor, kita cukup tambahkan prediktor tersebut di dalam model. Contohnya di sini, kita membuat model untuk melihat tinggi badan dan protein harian dalam gram sebagai dua prediktor dari berat badan.
multimodel_w <- lm(weight ~ height + daily_protein + age,
df)
multimodel_wsummary(multimodel_w)Selain menambah efek langsung dari satu variabel, kita juga bisa menambah interaksi antara dua variabel. Interaksi dapat ditambah dengan notasi *. Menambahkan interaksi dengan cara ini akan memberikan kita baik efek interaksin dari variabel yang dispesifikasikan, dalam kasus ini daily_protein:age, berikut juga efek langsung dari masing-masing variabel.
lm(weight ~ height + daily_protein * age,
df)Jika kita tidak ingin mendapatkan efek langsung dari masing-masing variabel, dan hanya ingin melihat interaksinya saja, kita harus menggunakan notasi :.
lm(weight ~ height + daily_protein:age,
df)Untuk melihat apakah model regresi yang baru lebih baik dibanding model regresi yang sebelumnya, kita bisa menggunakan fungsi anova() dari base R.
Fungsi anova() secara otomatis melakukan komparasi antara dua atau lebih model yang berbeda. Komparasi menggunakan fungsi ini bisa dilakukan jika model-model yang dibandingkan memiliki komponen yang sama. Artinya misalnya, model kedua harus memiliki semua variabel prediktor dan kriterion dari model pertama, dan model ketiga harus memiliki semua variabel dari model kedua, dan seterusnya.
anova(model_wh, multimodel_w)Dari komparasi ini kita bisa melihat bahwa model multiple regression memiliki fit yang lebih baik.
Kita bisa melihat distribusi residual dalam data dengan menggunakan fungsi hist dan menspesifikasikan data yang digambar sebagai residual dari model regresi kita.
hist(x = residuals(model_wh),
xlab = "Value of residual",
main = "",
breaks = 20)Selain itu, kita juga dapat menggunakan fungsi plot() dengan menggunakan perintah which = 2 untuk menggambar theoretical quantile dengan Q-Q plot dari residual.
plot(model_wh, which = 2)Sejauh ini, baik dari histogram maupun QQ plot terlihat bahwa data residual di model model_wh cukup normal. Histogram dari residual berbentuk distribusi unimodal. Q-Q plot menunjukkan residual terletak umumnya dalam garis theoretical quantile.
Kita dapat menggunakan plot() lagi untuk mendapat bantuan visual dalam menentukan apakah hubungan antar variabel di dalam model kita linear atau tidak. Dengan memasukkan perintah which = 1, kita dapat melihat scatter plot yang menggambarkan hubungan antara fitted values dan average residual. Di dalam model yang linear, hubungan antara keduanya akan konsisten, yang artinya garis di dalam scatterplot ini akan kurang lebih lurus.
plot(model_wh, which = 1)Dari plot residual vs. fitted di atas, garis di dalam model ini terlihat tidak terlalu menyimpang dari garis tengah yang uniform. Konsistensi ini berarti kurang lebih kita bisa menganggap hubungan antar variabel bersifat linear.
Jika kita memiliki model dengan lebih dari satu prediktor seperti halnya dalam model multiple regression, kita bisa menggunakan Tukey test untuk masing masing prediktor dan menggambar plot residual vs. fitted untuk setiap prediktor. Kita dapat menggunakan residualPlots() untuk melakukan keduanya. Fungsi residualPlots() ada di dalam package car.
install.packages('car')
library(car)
residualPlots(multimodel_w)Untuk melihat homogeneity dari varians di dalam model, kita dapat menggunakan plot() lagi. Kali ini, plot() digunakan untuk membuat scatterplot yang membandingkan standardized residuals dengan fitted values, yang dipanggil dengan menggunakan perintah which = 3. Garis horizontal yang konsisten mengindikasikan tidak adanya permasalahan homogeneity.
plot(model_wh, which = 3)Kita juga bisa melakukan non-constant variance test, dengan menggunakan ncvTest() dari package car. Jika tidak signifikan, maka kemungkinan asumsi homogeneity tidak bermasalah.
ncvTest(model_wh)Collinearity pada dasarnya adalah asumsi bahwa prediktor di dalam model tidak terlalu berkorelasi antara satu sama lain. Cara sederhana untuk menguji hal ini adalah dengan membuat tes korelasi.
cor(df[c("age", "height", "daily_protein")])Dari matriks korelasi ini, sepertinya tidak ada yang perlu dikhawatirkan. Kita bisa menambah analisis collinearity dengan tes variance inflation factors (VIF). Tes ini dapat dilakukan dengan fungsi vif(), yang merupakan bagian dari car.
vif(multimodel_w)Nilai VIF yang lebih tinggi umumnya dianggap lebih memiliki resiko multicollinearity. Nilai tersebut dapat diiterpretasikan sebagai seberapa besar korelasi antara variabel tersebut dengan variabel lain di dalam model meningkatkan varians di dalam nilai. Semakin tinggi, maka hubungan antara korelasi dan varians dapat dianggap semakin besar. Maka dari itu, nilai VIF yang relatif rendah dianggap lebih ideal.
Regresi umumnya dilaporkan menggunakan tabel berisi koefisien-koefisien regresi, yang isinya adalah koefisien untuk masing-masing prediktor beserta standard errror dan p-valuenya. Nilai-nilai ini sudah kita lihat sebelumnya, menggunakan fungsi summary(). Namun, jika kita melaporkan regresi dengan format di luar R, seperti di dalam dokumen yang dibuat di software word processor, kita tidak bisa menggunakan tabel dari summary() secara langsung. Butuh beberapa langkah agar kita bisa memindahkan tabel koefisien regresi di dalam R ke luar dengan rapih.
Pertama, kita bisa mengambil tabel tersebut menggunakan as.data.frame. Perhatikan bahwa di dalam baris ini kita mengambil koefisien dari model multiple regression, dengan cara menggunakan coef() dengan objek summary(multimodel_w) di dalamnya.
tab_coeff_multimodel_w <- as.data.frame(coef(summary(multimodel_w)))
tab_coeff_multimodel_wTabel ini sudha berisi hampir semua yang umumnya dibutuhkan dalam melaporkan regresi, seprti estimate dan p-value dari setiap koefisien. Sekarang, karena tabel tersebut sudah disimpan menjadi objek, kita bisa menyimpan tabel ini dengan write.table untuk mengkopi ke clipboard untuk di-paste ke spreadsheet yang sudah ada, atau write.csv jika kita ingin membuat spreadsheet baru.
# Copy-paste
write.table(tab_coeff_multimodel_w, "clipboard", sep="\t")
# Membuat spreadsheet baru
write.csv(tab_coeff_multimodel_w, "./Session 4/Tabel Multiple Regression.csv")Namun, tentunya tabel tersebut belum benar-benar sesuai dengan format tabel yang biasanya diperlukan di dalam pelaporan regresi. Kita masih perlu mengedit dan memasukkan nilai-nilai seperti DF di dalam tabel, dan format tabel juga harus kita buat sendiri.
Sebagai alternatif yang dari metode di atas, ada package yang secara spesifik dibuat untuk menulis sebuah tabel regresi yang sudah publication ready, yaitu stargazer.
install.packages("stargazer")
library(stargazer)Package stargazer() dapat digunakan untuk membuat tabel yang kurang lebih sudah cukup lengkap. Di contoh di bawah ini adalah skrip untuk membuat sebuah tabel dari model multiple regression yang kita gunakan. Perintah type digunakan untuk memberikan spesifikasi format file yang dibuat, dengan pilihan ASCII, HTML, atau LATEX. out adalah judul dari file yang akan dibuat. star.cutoffs merupakan perintah untuk menentukan pada p-value mana saja kita akan memberikan tanda bintang.
stargazer(multimodel_w,
type = "html",
out = "Multiple Regression Model.html",
star.cutoffs = c(0.05, 0.01, 0.001))Regresi linear dapat direpresentasikan dengan plot yang berisi scatterplot untuk tiap data point dan sebuah garis yang merepresentasikan koefisien regresi. ggplot, package visualisasi data di dalam tidyverse, dapat digunakan untuk membuat visualisasi ini.
ggplot sendiri merupakan package yang memiliki banyak sekali kegunaan, dan membutuhkan penjelasan yang jauh lebih mendalam. Namun pada dasarnya, script di bawah memiliki logika sebagai berikut. pertama, dibuat ‘kanvas’ dengan ggplot() di mana data dan variabel kita panggil. Lalu kanvas tersebut kita tambah dengan menggambar beberapa objek, seperti poin data dalam scatterplot, dan garis regresi. Setiap kali kita menambahkan objek di dalam kanvas, kita harus menggunakan tanda +. Pertama, kita tambahkan scatterplot untuk data, dengan menggunakan perintah geom_point(). Lalu, garis regresi digambar menggunakan geom_smooth(method='lm').
plot_model_wh <-ggplot(df,aes(weight, height)) +
geom_point() +
geom_smooth(method='lm')
plot_model_whBagaimana dengan multiple regression? Scatterplot dan garis umumnya terlalu sederhana, dan butuh diubah menjadi plot tiga dimensi untuk bisa merepresentasikan lebih dari satu koefisien regresi. Namun jika variabel yang digunakan banyak, metode tersebut tidak terlalu intuitif untuk dibaca dan cukup rumit untuk dibuat.
Alternatif yang bisa digunakan adalah forest plot. Forest plot adalah visualisasi data yang umumnya digunakan untuk meta analisis, ketika analis butuh melihat efek dari beberapa studi sekaligus. Namun, forest plot bisa juga digunakan untuk membaca model regresi, dengan menjabarkan beberapa koefisien, alih-alih studi yang berbeda-beda.
Package yang digunakan untuk membuat forest plot di R salah satunya adalah jtools. Untuk contoh, mari kita coba buat forest plot untuk koefisien-koefisien dari multimodel_w.
install.packages("jtools")
install.packages("ggstance") #Package ini adalah salah satu dependency dari jtools untuk membuat Forest Plot
library(jtools)Dari jtools, kita gunakan fungsi plot_summs untuk membuat forest plot dengan semua koefisien (minus intercept dari model). Perintah scale = TRUE digunakan untuk menstandardisasi skala antar variabel. Di dalam forest plot yang terlihat di bawah terdapat titik dan garis, di mana titik merepresentasikan koefisien, dan garis di titik tersebut merepresentasikan confidence interval 95%.
plot_summs(multimodel_w, scale = FALSE)