Week 5 Analisis Regresi untuk data lingkungan
Aulia Putri Andini
2025-09-29
Pendahuluan
Hubungan antara variabel lingkungan sering kali bersifat kompleks dan tidak selalu mengikuti pola linear. Pada dataset airquality, misalnya, konsentrasi Ozone diduga dipengaruhi oleh tingkat Solar.R. Secara teoritis, radiasi matahari dapat memicu pembentukan ozon, tetapi pola hubungan yang muncul di data bisa melengkung atau bervariasi pada rentang tertentu. Oleh karena itu, diperlukan analisis regresi dengan berbagai pendekatan—mulai dari model linear sederhana hingga spline—untuk memperoleh gambaran yang lebih akurat mengenai pola hubungan tersebut.
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.1 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(splines2)
df_airquality <- datasets::airquality
head(df_airquality)## Ozone Solar.R Wind Temp Month Day
## 1 41 190 7.4 67 5 1
## 2 36 118 8.0 72 5 2
## 3 12 149 12.6 74 5 3
## 4 18 313 11.5 62 5 4
## 5 NA NA 14.3 56 5 5
## 6 28 NA 14.9 66 5 6Cek apakah ada missing value
colSums(is.na(airquality))## Ozone Solar.R Wind Temp Month Day
## 37 7 0 0 0 0Dari hasil pemeriksaan, dataset airquality memiliki missing value pada variabel Ozone dan Solar.R. Informasi ini penting sebelum melakukan analisis, karena kita perlu menentukan strategi untuk menangani nilai yang hilang, misalnya dengan menghapus baris, mengisi dengan rata-rata/median, atau menggunakan metode imputasi lainnya.
# Ganti NA pada kolom Ozone dengan median
airquality$Ozone[is.na(airquality$Ozone)] <- median(airquality$Ozone, na.rm = TRUE)
# Ganti NA pada kolom Solar.R dengan median
airquality$Solar.R[is.na(airquality$Solar.R)] <- median(airquality$Solar.R, na.rm = TRUE)
# Cek lagi apakah masih ada missing value
colSums(is.na(airquality))## Ozone Solar.R Wind Temp Month Day
## 0 0 0 0 0 0median(…, na.rm = TRUE) menghitung median dengan mengabaikan NA. Baris yang NA kemudian diisi dengan nilai median dari kolom tersebut. Setelah diganti, hasil colSums(is.na()) akan menunjukkan 0 → artinya tidak ada lagi missing value.
Data Air Quality
ggplot(df_airquality,aes(x=Solar.R, y=Ozone)) +
geom_point(alpha=0.55, color="black") +
theme_bw() ## Warning: Removed 42 rows containing missing values or values outside the scale range
## (`geom_point()`).
- Dari grafik terlihat hubungan tidak sepenuhnya linear → ada indikasi
kurva. - Hal ini menjadi dasar kita perlu coba model non-linear (tangga,
spline).
Regresi Linear
mod_linear = lm(Ozone~Solar.R,data=df_airquality)
summary(mod_linear)##
## Call:
## lm(formula = Ozone ~ Solar.R, data = df_airquality)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.292 -21.361 -8.864 16.373 119.136
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.59873 6.74790 2.756 0.006856 **
## Solar.R 0.12717 0.03278 3.880 0.000179 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.33 on 109 degrees of freedom
## (42 observations deleted due to missingness)
## Multiple R-squared: 0.1213, Adjusted R-squared: 0.1133
## F-statistic: 15.05 on 1 and 109 DF, p-value: 0.0001793- Model regresi linear sederhana: Ozone = β0 + β1*Solar.R.
- summary() menampilkan koefisien, nilai p, R², dll.
- Interpretasi: jika β1 signifikan, maka ada hubungan linear antara Solar.R dan Ozone.
- Kelemahan: hubungan bisa saja non-linear → model ini bisa terlalu kaku.
ggplot(df_airquality,aes(x=Solar.R, y=Ozone)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~x,lty = 1, col = "red",se = F)+
theme_bw()## Warning: Removed 42 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 42 rows containing missing values or values outside the scale range
## (`geom_point()`).
- Dari grafik terlihat bahwa linear fit cukup oke di bagian tengah, tapi
di ekor (nilai kecil & besar) tidak mengikuti pola dengan baik →
indikasi perlu model lebih fleksibel.
Regresi Fungsi Tangga
mod_tangga = lm(Ozone ~ cut(Solar.R,5),data=df_airquality)
summary(mod_tangga)##
## Call:
## lm(formula = Ozone ~ cut(Solar.R, 5), data = df_airquality)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.743 -20.647 -6.437 14.853 112.257
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.167 7.071 2.145 0.034254 *
## cut(Solar.R, 5)(72.4,138] 10.271 10.308 0.996 0.321338
## cut(Solar.R, 5)(138,203] 37.214 9.637 3.862 0.000194 ***
## cut(Solar.R, 5)(203,269] 40.576 8.702 4.663 9.13e-06 ***
## cut(Solar.R, 5)(269,334] 29.690 9.637 3.081 0.002629 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30 on 106 degrees of freedom
## (42 observations deleted due to missingness)
## Multiple R-squared: 0.2167, Adjusted R-squared: 0.1871
## F-statistic: 7.331 on 4 and 106 DF, p-value: 2.984e-05- Cut(Solar.R, 3) membagi Solar.R ke dalam 3 interval sama lebar.
- Artinya model mengestimasi rata-rata Ozone berbeda-beda untuk tiap interval Solar.R.
- Cocok kalau kita hanya ingin melihat perbedaan rata-rata antar kelompok.
- Kelemahan: tidak halus (stepwise), prediksi bisa “loncat-loncat”.
ggplot(df_airquality,aes(x=Solar.R, y=Ozone)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,5),
lty = 1, col = "red",se = F)+
theme_bw()## Warning: Removed 42 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 42 rows containing missing values or values outside the scale range
## (`geom_point()`).
- Membagi Solar.R jadi 5 kategori (lebih detail). - Garis biru
menunjukkan rata-rata Ozone per interval. - Hasilnya “berundak” (fungsi
tangga). - Mudah dipahami, tapi kehilangan informasi variasi halus.
Regresi Spline
mod_spline3 <- lm(Ozone ~ bSpline(Solar.R, knots = c(10, 50, 100), degree = 3),
data = df_airquality)# bikin basis spline dulu
X <- bSpline(df_airquality$Solar.R, knots = c(10, 50, 100), degree = 3)
# masukin ke model
mod_spline3 <- lm(df_airquality$Ozone ~ X)
summary(mod_spline3)##
## Call:
## lm(formula = df_airquality$Ozone ~ X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.372 -18.820 -4.303 16.001 109.628
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.1447 28.9398 0.420 0.6756
## X1 2.8255 43.7722 0.065 0.9487
## X2 -3.1503 35.0715 -0.090 0.9286
## X3 7.3278 33.9333 0.216 0.8295
## X4 25.3733 37.4921 0.677 0.5001
## X5 81.0412 36.4927 2.221 0.0285 *
## X6 -0.3257 32.3370 -0.010 0.9920
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 29.7 on 104 degrees of freedom
## (42 observations deleted due to missingness)
## Multiple R-squared: 0.2468, Adjusted R-squared: 0.2033
## F-statistic: 5.678 on 6 and 104 DF, p-value: 3.895e-05Spline memecah rentang Solar.R menjadi beberapa segmen dengan knot (titik belok), lalu pasang polinomial di tiap segmen. bs() → B-spline: fleksibel, bisa melengkung mengikuti data. ns() → Natural spline: mirip, tapi di luar knot paling luar dipaksa linear (lebih stabil untuk ekstrapolasi). quantile() dipakai supaya knot sesuai distribusi data (tidak sembarang angka). Hasil summary() menunjukkan koefisien basis spline (tidak mudah diinterpretasi langsung, tapi penting untuk prediksi).
- Garis biru = cubic spline, merah = natural spline.
- Keduanya mengikuti pola data lebih baik daripada regresi linear.
- Natural spline lebih “jinak” di ekor → cocok untuk data lingkungan yang rawan outlier.
- Dari sini terlihat: spline memberi fleksibilitas untuk menangkap hubungan non-linear.