Pada tugas mata kuliah Pengantar Data Sains ini, diberikan perintah untuk membuat laporan analisis menggunakan R Markdown sebagai latihan menggabungkan penjelasan, kode R, dan output analisis dalam satu dokumen yang kemudian dirender ke format HTML dan dipublikasikan melalui RPubs.
Dalam tugas ini, penulis menganalisis dataset bawaan R yaitu dataset airquality yang berisi data kualitas udara harian di New York pada tahun 1973. Dataset ini memuat variabel Ozone, Solar.R, Wind, dan Temp yang berkaitan dengan proses pembentukan ozon troposfer. Sehingga analisis kualitas udara penting dilakukan karena ozon permukaan dapat berdampak buruk pada kesehatan dan lingkungan.
Metode yang dipilih dan digunakan adalah regresi linear berganda untuk melihat apakah terdapat pengaruh beberapa faktor meteorologi secara simultan terhadap konsentrasi Ozone. Proses analisis juga mencakup EDA dan uji asumsi klasik guna memastikan model yang digunakan valid.
Tugas ini bertujuan untuk memahami langkah-langkah analisis data menggunakan R sekaligus melatih pembuatan laporan yang terstruktur dan otomatis melalui R Markdown.
library(datasets)
library(dplyr)
library(ggplot2)
library(car)
library(lmtest)
data("airquality")
df <- airquality
summary(df)
## Ozone Solar.R Wind Temp
## Min. : 1.00 Min. : 7.0 Min. : 1.700 Min. :56.00
## 1st Qu.: 18.00 1st Qu.:115.8 1st Qu.: 7.400 1st Qu.:72.00
## Median : 31.50 Median :205.0 Median : 9.700 Median :79.00
## Mean : 42.13 Mean :185.9 Mean : 9.958 Mean :77.88
## 3rd Qu.: 63.25 3rd Qu.:258.8 3rd Qu.:11.500 3rd Qu.:85.00
## Max. :168.00 Max. :334.0 Max. :20.700 Max. :97.00
## NA's :37 NA's :7
## Month Day
## Min. :5.000 Min. : 1.0
## 1st Qu.:6.000 1st Qu.: 8.0
## Median :7.000 Median :16.0
## Mean :6.993 Mean :15.8
## 3rd Qu.:8.000 3rd Qu.:23.0
## Max. :9.000 Max. :31.0
##
# Cek NA
t(colSums(is.na(df)))
## Ozone Solar.R Wind Temp Month Day
## [1,] 37 7 0 0 0 0
# Hapus NA
air <- na.omit(df)
t(colSums(is.na(air)))
## Ozone Solar.R Wind Temp Month Day
## [1,] 0 0 0 0 0 0
Pada dataset asli Airquality memiliki NA pada Ozone dan Solar.R, sehingga perlu dibersihkan/dibuang agar model dapat dijalankan dan analisi dapat dilakukan.
summary(air)
## Ozone Solar.R Wind Temp
## Min. : 1.0 Min. : 7.0 Min. : 2.30 Min. :57.00
## 1st Qu.: 18.0 1st Qu.:113.5 1st Qu.: 7.40 1st Qu.:71.00
## Median : 31.0 Median :207.0 Median : 9.70 Median :79.00
## Mean : 42.1 Mean :184.8 Mean : 9.94 Mean :77.79
## 3rd Qu.: 62.0 3rd Qu.:255.5 3rd Qu.:11.50 3rd Qu.:84.50
## Max. :168.0 Max. :334.0 Max. :20.70 Max. :97.00
## Month Day
## Min. :5.000 Min. : 1.00
## 1st Qu.:6.000 1st Qu.: 9.00
## Median :7.000 Median :16.00
## Mean :7.216 Mean :15.95
## 3rd Qu.:9.000 3rd Qu.:22.50
## Max. :9.000 Max. :31.00
str(air)
## 'data.frame': 111 obs. of 6 variables:
## $ Ozone : int 41 36 12 18 23 19 8 16 11 14 ...
## $ Solar.R: int 190 118 149 313 299 99 19 256 290 274 ...
## $ Wind : num 7.4 8 12.6 11.5 8.6 13.8 20.1 9.7 9.2 10.9 ...
## $ Temp : int 67 72 74 62 65 59 61 69 66 68 ...
## $ Month : int 5 5 5 5 5 5 5 5 5 5 ...
## $ Day : int 1 2 3 4 7 8 9 12 13 14 ...
## - attr(*, "na.action")= 'omit' Named int [1:42] 5 6 10 11 25 26 27 32 33 34 ...
## ..- attr(*, "names")= chr [1:42] "5" "6" "10" "11" ...
Statistik deskriptif digunakan untuk memberikan gambaran awal mengenai karakteristik data pada setiap variabel. Output fungsi summary() menampilkan nilai minimum, maximum, median, mean, serta kuartil, sehingga kita dapat melihat persebaran data dan potensi adanya pencilan (outlier). Sementara itu, fungsi str() menunjukkan struktur data, tipe variabel, dan jumlah observasi, sehingga membantu memastikan bahwa setiap variabel memiliki tipe data yang sesuai untuk analisis. Melalui statistik deskriptif ini, kita dapat memahami kondisi umum dataset sebelum melakukan analisis lebih lanjut.
par(mfrow=c(2,2))
hist(air$Ozone, main="Histogram Ozone", xlab="Ozone")
hist(air$Temp, main="Histogram Temp", xlab="Temp")
hist(air$Wind, main="Histogram Wind", xlab="Wind")
hist(air$Solar.R, main="Histogram Solar.R", xlab="Solar.R")
par(mfrow=c(1,1))
Histogram digunakan untuk melihat bentuk distribusi dari setiap variabel dalam dataset. Melalui grafik ini, kita dapat mengidentifikasi apakah data cenderung simetris, miring (skewed) ke kiri atau kanan, serta apakah terdapat outlier atau nilai ekstrem. Analisis distribusi ini penting dilakukan sebelum membangun model regresi, karena bentuk distribusi dapat memengaruhi interpretasi maupun pemenuhan asumsi model.
plot(air$Temp, air$Ozone, main="Ozone vs Temp", xlab="Temp", ylab="Ozone")
plot(air$Wind, air$Ozone, main="Ozone vs Wind", xlab="Wind", ylab="Ozone")
plot(air$Solar.R, air$Ozone, main="Ozone vs Solar.R", xlab="Solar.R", ylab="Ozone")
cor(air[,c("Ozone","Temp","Wind","Solar.R")])
## Ozone Temp Wind Solar.R
## Ozone 1.0000000 0.6985414 -0.6124966 0.3483417
## Temp 0.6985414 1.0000000 -0.4971897 0.2940876
## Wind -0.6124966 -0.4971897 1.0000000 -0.1271835
## Solar.R 0.3483417 0.2940876 -0.1271835 1.0000000
Scatterplot digunakan untuk melihat pola hubungan antara variabel prediktor dengan variabel respon (Ozone). Dari grafik:
Temp vs Ozone: terlihat kecenderungan hubungan positif, di mana peningkatan suhu berkaitan dengan meningkatnya konsentrasi Ozone.
Wind vs Ozone: menunjukkan pola negatif, yaitu semakin tinggi kecepatan angin, konsentrasi Ozone cenderung menurun.
Solar.R vs Ozone: tampak adanya hubungan positif meskipun dengan persebaran data yang lebih bervariasi.
Selain itu, perhitungan korelasi memperkuat pengamatan visual pada scatterplot. Korelasi yang bernilai positif menunjukkan hubungan searah, sementara korelasi negatif menunjukkan hubungan yang berlawanan arah. Informasi ini membantu menentukan prediktor mana yang berpotensi signifikan dalam model regresi berganda.
Kita membangun model untuk melihat pengaruh Temp, Wind, dan Solar.R terhadap Ozone.
model <- lm(Ozone ~ Temp + Wind + Solar.R, data = air)
summary(model)
##
## Call:
## lm(formula = Ozone ~ Temp + Wind + Solar.R, data = air)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.485 -14.219 -3.551 10.097 95.619
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -64.34208 23.05472 -2.791 0.00623 **
## Temp 1.65209 0.25353 6.516 2.42e-09 ***
## Wind -3.33359 0.65441 -5.094 1.52e-06 ***
## Solar.R 0.05982 0.02319 2.580 0.01124 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.18 on 107 degrees of freedom
## Multiple R-squared: 0.6059, Adjusted R-squared: 0.5948
## F-statistic: 54.83 on 3 and 107 DF, p-value: < 2.2e-16
Dari hasil summary(model):
Temp berpengaruh positif → semakin tinggi suhu, Ozone meningkat.
Wind berpengaruh negatif → angin membantu mengurangi konsentrasi Ozone.
Solar.R dilihat dari p-value → signifikan jika p-value < 0.05, tidak berpengaruh jika p-value > 0.05.
Model ini menunjukkan variabel mana yang paling berkontribusi terhadap perubahan Ozone.
Asumsi linearitas menyatakan bahwa hubungan antara variabel independen (Temp, Wind, Solar.R) dengan variabel dependen (Ozone) harus bersifat linear. Artinya, perubahan pada variabel prediktor diharapkan berdampak secara proporsional terhadap perubahan pada Ozone. Jika hubungan ini tidak linear, maka model regresi linear tidak dapat menangkap pola data dengan baik dan menghasilkan estimasi yang bias.
Untuk memeriksa linearitas, digunakan beberapa pendekatan:
Plot Residual vs Fitted Plot ini digunakan untuk melihat apakah residual menyebar secara acak di sekitar garis horizontal nol. Jika pola tertentu seperti lengkungan (curve) muncul, maka hubungan mungkin tidak linear.
Component + Residual Plot (Partial Residual Plot) Plot ini menampilkan hubungan antara setiap variabel prediktor dengan residual yang sudah disesuaikan. Jika titik-titik membentuk pola melengkung atau tidak konsisten, maka hubungan prediktor terhadap Ozone tidak linear.
Scatterplot Sederhana (EDA) Melihat pola awal hubungan antara Ozone dan tiap variabel meteorologi untuk mendeteksi indikasi awal ketidaklinearan.
Berikut kode untuk mengecek linearitas:
par(mfrow = c(2,2))
plot(model)
Hasil dari plot ini kemudian dianalisis untuk memastikan apakah terdapat
pola melengkung pada residual, titik-titik yang terstruktur, atau
indikasi ketidaklinearan lainnya.
resid <- model$residuals
shapiro.test(resid)
##
## Shapiro-Wilk normality test
##
## data: resid
## W = 0.91709, p-value = 3.618e-06
par(mfrow=c(1,2))
hist(resid, main="Histogram Residual")
qqnorm(resid); qqline(resid)
par(mfrow=c(1,1))
Uji Shapiro–Wilk digunakan untuk mengecek apakah residual berdistribusi normal. Jika p-value > 0.05, maka residual dapat dianggap normal. Grafik histogram dan QQ-plot juga membantu melihat apakah pola residual mengikuti garis normal.
bptest(model)
##
## studentized Breusch-Pagan test
##
## data: model
## BP = 5.0554, df = 3, p-value = 0.1678
Uji Breusch–Pagan digunakan untuk mengecek apakah varians residual konstan. Jika p-value < 0.05, berarti terdapat heteroskedastisitas (varians residual tidak homogen). Jika p-value ≥ 0.05, maka model dianggap tidak mengalami masalah heteroskedastisitas.
vif(model)
## Temp Wind Solar.R
## 1.431367 1.329070 1.095253
Uji VIF digunakan untuk melihat apakah antar variabel prediktor saling berkorelasi kuat. Jika VIF > 10, berarti terdapat multikolinearitas tinggi. Jika VIF masih di bawah 10, maka model dianggap aman dari masalah multikolinearitas.
dwtest(model)
##
## Durbin-Watson test
##
## data: model
## DW = 1.9355, p-value = 0.3347
## alternative hypothesis: true autocorrelation is greater than 0
Uji Durbin–Watson digunakan untuk mengecek apakah residual saling berkorelasi. Jika p-value < 0.05, berarti terdapat autokorelasi, sedangkan p-value ≥ 0.05 menunjukkan tidak ada masalah autokorelasi pada residual.
par(mfrow=c(2,2))
plot(model)
par(mfrow=c(1,1))
Plot diagnostik dari plot(model) menampilkan empat grafik utama yang membantu mengevaluasi asumsi regresi:
Residual vs Fitted → digunakan untuk memeriksa apakah hubungan antara variabel prediktor dan respons bersifat linear serta apakah varians residual konstan (homoskedastisitas).
Q-Q Plot → mengecek apakah residual mengikuti distribusi normal.
Scale-Location → menilai penyebaran residual agar tetap konsisten di seluruh rentang prediksi.
Residual vs Leverage → mengidentifikasi outlier atau titik yang memiliki pengaruh besar terhadap model.
Dengan keempat plot ini, kita dapat memastikan model regresi berganda yang dibangun valid dan memenuhi asumsi klasik.
Berdasarkan ringkasan model:
Temp signifikan → peningkatan suhu meningkatkan konsentrasi Ozone.
Wind signifikan → kecepatan angin cenderung menurunkan Ozone.
Solar.R dapat signifikan atau tidak tergantung pada p-value.
Adjusted R² menunjukkan seberapa besar variasi Ozone dijelaskan oleh semua variabel prediktor secara simultan.
Interpretasi tambahan:
Jika Adjusted R² rendah, masih banyak faktor lain di luar model yang mempengaruhi Ozone.
Tanda positif atau negatif pada koefisien sesuai dugaan fisik: suhu tinggi mempercepat reaksi pembentukan ozon, angin membantu menyebarkan polutan.
Berdasarkan analisis: