Lampiran Kode Projek I BIOSTAT Kel. 2
Penulis.
1. Persiapan dan Import Data
Langkah pertama adalah memanggil package yang dibutuhkan untuk manipulasi data, visualisasi, uji regresi, dan pengolahan data spasial. Jika Anda belum menginstalnya, jalankan perintah install.packages() terlebih dahulu.
# Memuat library yang dibutuhkan
library(readxl)## Warning: package 'readxl' was built under R version 4.4.3library(dplyr)## Warning: package 'dplyr' was built under R version 4.4.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)## Warning: package 'ggplot2' was built under R version 4.4.3library(tidyr)## Warning: package 'tidyr' was built under R version 4.4.3library(car)## Warning: package 'car' was built under R version 4.4.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.4.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(lmtest)## Warning: package 'lmtest' was built under R version 4.4.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.4.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(sf)## Warning: package 'sf' was built under R version 4.4.3## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUESelanjutnya, kita impor dataset kesehatan tahun 2025 dan melakukan pemeriksaan awal untuk melihat struktur dan dimensi data.
# Mengimpor data dengan konfigurasi direktori komputer peneliti
datasehat <- read_excel("C:\\Users\\aethe\\Downloads\\biosatat\\Data_Kesehatan_2025.xlsx")
# Memeriksa struktur awal data
str(datasehat)## tibble [38 × 5] (S3: tbl_df/tbl/data.frame)
## $ Provinsi : chr [1:38] "ACEH" "SUMATERA UTARA" "SUMATERA BARAT" "RIAU" ...
## $ Y (Keluhan Kesehatan): num [1:38] 24.6 24.3 28.5 21.7 25.9 ...
## $ X1 (Merokok) : num [1:38] 28.2 26.4 30.5 28.9 30.3 ...
## $ X2 (Sanitasi Layak) : chr [1:38] "82.21" "87.47" "74.59" "91.21" ...
## $ X3 (Air Minum Layak) : num [1:38] 92.1 93.2 87.3 92.2 82.5 ...dim(datasehat)## [1] 38 5summary(datasehat)## Provinsi Y (Keluhan Kesehatan) X1 (Merokok) X2 (Sanitasi Layak)
## Length:38 Min. : 8.57 Min. :19.02 Length:38
## Class :character 1st Qu.:20.23 1st Qu.:25.23 Class :character
## Mode :character Median :24.46 Median :27.73 Mode :character
## Mean :24.23 Mean :27.44
## 3rd Qu.:29.32 3rd Qu.:29.99
## Max. :42.00 Max. :34.25
## X3 (Air Minum Layak)
## Min. :32.89
## 1st Qu.:83.29
## Median :89.94
## Mean :87.61
## 3rd Qu.:94.90
## Max. :99.982. Data Cleaning
Sebelum dianalisis, data perlu dibersihkan dengan menghapus baris kosong, mengubah tipe data ke numerik, dan menyederhanakan nama kolom agar lebih mudah dipanggil di fungsi-fungsi selanjutnya.
# Membersihkan data
datasehat <- datasehat %>%
# Menghapus baris yang tidak memiliki nama provinsi
filter(!is.na(Provinsi)) %>%
# Mengubah variabel menjadi numerik
mutate(
`Y (Keluhan Kesehatan)` = as.numeric(`Y (Keluhan Kesehatan)`),
`X1 (Merokok)` = as.numeric(`X1 (Merokok)`),
`X2 (Sanitasi Layak)` = as.numeric(`X2 (Sanitasi Layak)`),
`X3 (Air Minum Layak)` = as.numeric(`X3 (Air Minum Layak)`)
) %>%
# Mengganti nama variabel
rename(
Keluhan = `Y (Keluhan Kesehatan)`,
Merokok = `X1 (Merokok)`,
Sanitasi = `X2 (Sanitasi Layak)`,
Air_Minum = `X3 (Air Minum Layak)`
)
# Mengecek kembali ketersediaan data dan missing values
colSums(is.na(datasehat))## Provinsi Keluhan Merokok Sanitasi Air_Minum
## 0 0 0 0 0sum(duplicated(datasehat$Provinsi))## [1] 03. Statistika Deskriptif
Bagian ini merangkum ukuran pemusatan dan penyebaran data untuk melihat gambaran umum variabel kesehatan.
# Ringkasan data untuk kolom analisis
summary(datasehat[, c("Keluhan", "Merokok", "Sanitasi", "Air_Minum")])## Keluhan Merokok Sanitasi Air_Minum
## Min. : 8.57 Min. :19.02 Min. :16.34 Min. :32.89
## 1st Qu.:20.23 1st Qu.:25.23 1st Qu.:81.15 1st Qu.:83.29
## Median :24.46 Median :27.73 Median :85.99 Median :89.94
## Mean :24.23 Mean :27.44 Mean :82.96 Mean :87.61
## 3rd Qu.:29.32 3rd Qu.:29.99 3rd Qu.:90.12 3rd Qu.:94.90
## Max. :42.00 Max. :34.25 Max. :98.20 Max. :99.98# Standar deviasi dari masing-masing variabel
sapply(datasehat[, c("Keluhan","Merokok","Sanitasi","Air_Minum")], sd)## Keluhan Merokok Sanitasi Air_Minum
## 6.647625 3.330466 14.960211 11.6148754. Visualisasi Data
4.1 Peta Spasial Keluhan Kesehatan
Memetakan sebaran keluhan kesehatan di seluruh Indonesia menggunakan data geojson.
# Mengimpor data peta spasial
file_peta <- "C:\\Users\\aethe\\Downloads\\biosatat\\indonesia_38_provinces.geojson"
peta_prov <- st_read(file_peta, quiet = TRUE)
# Menyiapkan data agar nama provinsi cocok (huruf kapital)
data_peta <- datasehat %>%
mutate(PROVINSI = toupper(trimws(Provinsi))) %>%
select(PROVINSI, Keluhan)
peta_prov <- peta_prov %>%
mutate(PROVINSI = toupper(trimws(PROVINSI)))
# Menggabungkan data kesehatan dengan data peta spasial
peta_keluhan <- peta_prov %>% left_join(data_peta, by = "PROVINSI")
# Membuat kategorisasi rentang keluhan
peta_keluhan <- peta_keluhan %>%
mutate(
Kategori = case_when(
Keluhan < 18 ~ "Rendah",
Keluhan < 24 ~ "Menengah",
Keluhan < 30 ~ "Tinggi",
Keluhan >= 30 ~ "Sangat tinggi"
)
)
peta_keluhan$Kategori <- factor(
peta_keluhan$Kategori,
levels = c("Rendah", "Menengah", "Tinggi", "Sangat tinggi")
)
# Plot peta
ggplot(peta_keluhan) +
geom_sf(aes(fill = Kategori), color = "white", linewidth = 0.3) +
scale_fill_manual(
values = c("Rendah" = "#D7191C", "Menengah" = "#F58200",
"Tinggi" = "#FFD54F", "Sangat tinggi" = "#003B5C"),
na.value = "grey80"
) +
theme_void() +
theme(
legend.position = "none",
plot.background = element_rect(fill = "transparent", color = NA),
panel.background = element_rect(fill = "transparent", color = NA),
plot.margin = margin(0, 0, 0, 0)
)
4.2 Hubungan Air Minum dan Keluhan Kesehatan
Scatterplot berikut menunjukkan arah hubungan antara ketersediaan air minum layak dan keluhan kesehatan masyarakat.
ggplot(datasehat, aes(x = Air_Minum, y = Keluhan)) +
geom_point(size = 3, color = "#003049") +
geom_smooth(method = "lm", se = TRUE, color = "#D62828", fill = "#FCBF49") +
labs(x = "Air minum layak (%)", y = "Keluhan kesehatan (%)") +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.line = element_line(color = "#003049", linewidth = 0.5),
axis.ticks = element_line(color = "#003049"),
axis.text = element_text(color = "#003049"),
axis.title = element_text(color = "#003049", face = "bold")
)## `geom_smooth()` using formula = 'y ~ x'
4.3 Komparasi Ekstrem Sanitasi Layak
Grafik ini membandingkan satu provinsi dengan persentase sanitasi layak terendah dan tertinggi.
# Mengambil nilai ekstrim
sanitasi_terendah <- datasehat %>%
filter(!is.na(Sanitasi)) %>%
slice_min(Sanitasi, n = 1, with_ties = FALSE) %>%
mutate(Kategori = "Terendah")
sanitasi_tertinggi <- datasehat %>%
filter(!is.na(Sanitasi)) %>%
slice_max(Sanitasi, n = 1, with_ties = FALSE) %>%
mutate(Kategori = "Tertinggi")
data_dua_batang <- bind_rows(sanitasi_terendah, sanitasi_tertinggi) %>%
mutate(Kategori = factor(Kategori, levels = c("Terendah", "Tertinggi")))
# Membuat bar chart
ggplot(data_dua_batang, aes(x = Kategori, y = Sanitasi, fill = Kategori)) +
geom_col(width = 0.48, color = NA) +
geom_text(aes(label = paste0(round(Sanitasi, 1), "%")),
vjust = -0.6, size = 9, fontface = "bold", color = "#003049") +
geom_text(aes(y = 0, label = Provinsi),
vjust = 1.8, size = 3.5, fontface = "bold", color = "#003049") +
scale_fill_manual(values = c("Terendah" = "#F77F00", "Tertinggi" = "#003049")) +
scale_y_continuous(limits = c(0, 110), breaks = seq(0, 100, 20), expand = expansion(mult = c(0, 0.02))) +
labs(
title = "SANITASI LAYAK ANTARPROVINSI",
subtitle = "Perbandingan nilai terendah dan tertinggi | Indonesia 2025",
x = NULL, y = "Sanitasi layak (%)"
) +
theme_minimal() +
theme(
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.major.y = element_line(color = "#EAE2B7", linewidth = 0.4),
plot.title = element_text(color = "#003049", face = "bold", size = 15, hjust = 0.5),
plot.subtitle = element_text(color = "#003049", size = 9, hjust = 0.5),
axis.title.y = element_text(color = "#003049", face = "bold"),
axis.text.x = element_blank(),
axis.text.y = element_text(color = "#003049"),
legend.position = "none"
)
4.4 6 Provinsi dengan Keluhan Tertinggi
# Menyiapkan top 6 data
top6_keluhan <- datasehat %>%
arrange(desc(Keluhan)) %>%
slice_head(n = 6) %>%
mutate(Provinsi = factor(Provinsi, levels = rev(Provinsi)))
ggplot(top6_keluhan, aes(x = Keluhan, y = Provinsi)) +
geom_col(fill = "#D62828", width = 0.80) +
geom_text(aes(label = paste0(round(Keluhan, 1), "%")),
hjust = -0.15, size = 4, fontface = "bold", color = "#003B5C") +
scale_x_continuous(expand = expansion(mult = c(0, 0.15))) +
labs(
title = "6 Provinsi dengan Keluhan Kesehatan Tertinggi",
subtitle = "Persentase keluhan kesehatan tahun 2025",
x = "Keluhan kesehatan (%)", y = NULL
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 16, color = "#003B5C", hjust = 0),
plot.subtitle = element_text(color = "#666666", hjust = 0),
axis.title.x = element_text(face = "bold", color = "#003B5C"),
axis.text.y = element_text(face = "bold", color = "#003B5C"),
panel.grid.minor = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_line(color = "#EAE2B7", linewidth = 0.4)
)
4.5 Distribusi Kebiasaan Merokok
# Membuat ringkasan kategori
data_donut <- datasehat %>%
mutate(
Kategori_Merokok = case_when(
Merokok < 20 ~ "Rendah",
Merokok < 30 ~ "Sedang",
Merokok >= 30 ~ "Tinggi",
TRUE ~ NA_character_
)
) %>%
count(Kategori_Merokok) %>%
mutate(Persentase = n / sum(n) * 100)
# Membuat donut chart
ggplot(data_donut, aes(x = 2, y = n, fill = Kategori_Merokok)) +
geom_col(width = 1, color = "white") +
coord_polar(theta = "y") +
xlim(0.5, 2.5) +
geom_text(
aes(label = paste0(n, " provinsi\n", round(Persentase, 1), "%")),
position = position_stack(vjust = 0.5), size = 4
) +
scale_fill_manual(
values = c("Rendah" = "#003049", "Sedang" = "#F77F00", "Tinggi" = "#D62828")
) +
labs(
title = "Distribusi Provinsi Berdasarkan Tingkat Merokok",
subtitle = "Indonesia Tahun 2025",
fill = "Kategori Merokok"
) +
theme_void() +
theme(
plot.title = element_text(face = "bold", size = 15, hjust = 0.5),
plot.subtitle = element_text(size = 11, hjust = 0.5),
legend.position = "right"
)
5. Analisis Korelasi Pearson
Melihat matriks korelasi antarvariabel serta uji signifikansinya secara spesifik.
# Matriks Korelasi Keseluruhan
cor_matrix <- cor(datasehat[, c("Keluhan", "Merokok", "Sanitasi", "Air_Minum")], method = "pearson")
round(cor_matrix, 3)## Keluhan Merokok Sanitasi Air_Minum
## Keluhan 1.000 0.408 0.505 0.337
## Merokok 0.408 1.000 -0.061 -0.090
## Sanitasi 0.505 -0.061 1.000 0.758
## Air_Minum 0.337 -0.090 0.758 1.000# Uji Signifikansi satu per satu
cor.test(datasehat$Keluhan, datasehat$Merokok)##
## Pearson's product-moment correlation
##
## data: datasehat$Keluhan and datasehat$Merokok
## t = 2.6806, df = 36, p-value = 0.01102
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1014561 0.6436585
## sample estimates:
## cor
## 0.4079094cor.test(datasehat$Keluhan, datasehat$Sanitasi)##
## Pearson's product-moment correlation
##
## data: datasehat$Keluhan and datasehat$Sanitasi
## t = 3.514, df = 36, p-value = 0.001211
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2214681 0.7102996
## sample estimates:
## cor
## 0.5053714cor.test(datasehat$Keluhan, datasehat$Air_Minum)##
## Pearson's product-moment correlation
##
## data: datasehat$Keluhan and datasehat$Air_Minum
## t = 2.1442, df = 36, p-value = 0.03884
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.01886808 0.59246713
## sample estimates:
## cor
## 0.33652166. Pemodelan Regresi Linear Berganda
Kita membuat model regresi untuk memprediksi keluhan kesehatan berdasarkan persentase merokok, sanitasi layak, dan air minum layak.
# Membuat model
model <- lm(Keluhan ~ Merokok + Sanitasi + Air_Minum, data = datasehat)
summary(model)##
## Call:
## lm(formula = Keluhan ~ Merokok + Sanitasi + Air_Minum, data = datasehat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.4283 -2.1856 -0.4558 3.2951 12.4563
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.93388 9.93535 -1.805 0.07993 .
## Merokok 0.87368 0.25480 3.429 0.00160 **
## Sanitasi 0.25826 0.08663 2.981 0.00528 **
## Air_Minum -0.03698 0.11182 -0.331 0.74288
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.14 on 34 degrees of freedom
## Multiple R-squared: 0.4505, Adjusted R-squared: 0.402
## F-statistic: 9.292 on 3 and 34 DF, p-value: 0.0001251# Mengeluarkan interval kepercayaan
confint(model)## 2.5 % 97.5 %
## (Intercept) -38.12494874 2.2571882
## Merokok 0.35586775 1.3914914
## Sanitasi 0.08220624 0.4343054
## Air_Minum -0.26423824 0.1902716Persamaan Regresi Berdasarkan koefisien yang didapatkan, kita dapat menyusun persamaan matematisnya:
koef <- coef(model)
cat("Persamaan regresi:\nYhat = ",
round(koef[1], 3), " + ",
round(koef[2], 3), "X1(Merokok) + ",
round(koef[3], 3), "X2(Sanitasi) + ",
round(koef[4], 3), "X3(Air_Minum)\n", sep = "")## Persamaan regresi:
## Yhat = -17.934 + 0.874X1(Merokok) + 0.258X2(Sanitasi) + -0.037X3(Air_Minum)7. Pengujian Hipotesis Regresi
Uji Simultan (Uji F) Digunakan untuk melihat apakah minimal satu variabel independen berpengaruh secara signifikan.
anova(model)## Analysis of Variance Table
##
## Response: Keluhan
## Df Sum Sq Mean Sq F value Pr(>F)
## Merokok 1 272.06 272.06 10.2957 0.0029060 **
## Sanitasi 1 461.68 461.68 17.4718 0.0001928 ***
## Air_Minum 1 2.89 2.89 0.1094 0.7428808
## Residuals 34 898.43 26.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1F_stat <- summary(model)$fstatistic[1]
df1 <- summary(model)$fstatistic[2]
df2 <- summary(model)$fstatistic[3]
p_F <- pf(F_stat, df1, df2, lower.tail = FALSE)
cat("Nilai F =", round(F_stat, 3), "\n")## Nilai F = 9.292cat("P-value Uji F =", p_F, "\n")## P-value Uji F = 0.0001250685Koefisien Determinasi (\(R^2\)) Melihat seberapa besar varians dari keluhan kesehatan mampu dijelaskan oleh model ini.
cat("R-squared =", summary(model)$r.squared, "\n")## R-squared = 0.4505213cat("Adjusted R-squared =", summary(model)$adj.r.squared, "\n")## Adjusted R-squared = 0.40203788. Uji Asumsi Klasik Regresi
Agar hasil analisis regresi OLS valid, kita perlu menguji beberapa asumsi dasar.
# Uji normalitas residual (H0: Residual berdistribusi normal)
shapiro.test(residuals(model))##
## Shapiro-Wilk normality test
##
## data: residuals(model)
## W = 0.97521, p-value = 0.5498# Uji heteroskedastisitas (H0: Homoskedastisitas / varians konstan)
bptest(model)##
## studentized Breusch-Pagan test
##
## data: model
## BP = 1.2541, df = 3, p-value = 0.7401# Uji multikolinearitas (Batas toleransi VIF < 10)
vif(model)## Merokok Sanitasi Air_Minum
## 1.008317 2.351735 2.3621029. Menghitung Prediksi dan Residual
Langkah terakhir, kita menyematkan nilai fitted (prediksi model) dan sisaan (residual) ke dalam data utama untuk evaluasi lanjutan jika diperlukan.
# Menambahkan kolom Prediksi dan Residual
datasehat <- datasehat %>%
mutate(
Prediksi = predict(model),
Residual = residuals(model)
)
# Melihat head dari data baru
head(datasehat %>% select(Provinsi, Keluhan, Prediksi, Residual))## # A tibble: 6 × 4
## Provinsi Keluhan Prediksi Residual
## <chr> <dbl> <dbl> <dbl>
## 1 ACEH 24.6 24.5 0.0783
## 2 SUMATERA UTARA 24.3 24.3 0.0387
## 3 SUMATERA BARAT 28.5 24.7 3.74
## 4 RIAU 21.7 27.4 -5.70
## 5 JAMBI 26.0 27.7 -1.71
## 6 SUMATERA SELATAN 26.2 28.0 -1.87