Analisis Multiple Linear Regression (MLR) Data Balita
Jamottan Fierdianto Sihite
2026-09-22
1 1. Persiapan
Analisis ini menggunakan data balita dengan variabel Umur (bulan), Jenis Kelamin, Tinggi Badan (cm), dan Status Gizi.
Untuk Multiple Linear Regression, variabel respons yang digunakan adalah Tinggi Badan (cm). Variabel prediktornya adalah Umur (bulan) dan Jenis Kelamin.
Catatan:
Status Gizitidak dimasukkan sebagai prediktor karena status gizi pada data ini berkaitan langsung dengan ukuran pertumbuhan seperti tinggi badan. Memasukkannya sebagai prediktor tinggi badan dapat menyebabkan hubungan yang tidak tepat/leakage. Status gizi tetap digunakan untuk eksplorasi dan visualisasi.
# Jika package belum tersedia, jalankan:
# install.packages(c("readxl", "dplyr", "ggplot2", "car", "lmtest",
# "nortest", "corrplot", "broom", "knitr"))
library(readxl)
library(dplyr)
library(ggplot2)
library(car)
library(lmtest)
library(nortest)
library(corrplot)
library(broom)
library(knitr)
set.seed(123)2 2. Import Data
Simpan file data_balita_rapi(1).xlsx pada folder yang
sama dengan file R Markdown.
data <- read_excel("C:/Users/user/Downloads/data_balita_rapi.xlsx")
# Melihat struktur data
str(data)## tibble [120,999 × 4] (S3: tbl_df/tbl/data.frame)
## $ Umur (bulan) : num [1:120999] 0 0 0 0 0 0 0 0 0 0 ...
## $ Jenis Kelamin : chr [1:120999] "laki-laki" "laki-laki" "laki-laki" "laki-laki" ...
## $ Tinggi Badan (cm): num [1:120999] 44.6 56.7 46.9 47.5 42.7 ...
## $ Status Gizi : chr [1:120999] "stunted" "tinggi" "normal" "normal" ...# Enam baris pertama
head(data)## # A tibble: 6 × 4
## `Umur (bulan)` `Jenis Kelamin` `Tinggi Badan (cm)` `Status Gizi`
## <dbl> <chr> <dbl> <chr>
## 1 0 laki-laki 44.6 stunted
## 2 0 laki-laki 56.7 tinggi
## 3 0 laki-laki 46.9 normal
## 4 0 laki-laki 47.5 normal
## 5 0 laki-laki 42.7 severely stunted
## 6 0 laki-laki 44.3 stunted# Ukuran data
dim(data)## [1] 120999 4# Nama variabel
names(data)## [1] "Umur (bulan)" "Jenis Kelamin" "Tinggi Badan (cm)"
## [4] "Status Gizi"3 3. Pembersihan dan Persiapan Data
data <- data %>%
rename(
umur_bulan = `Umur (bulan)`,
jenis_kelamin = `Jenis Kelamin`,
tinggi_badan = `Tinggi Badan (cm)`,
status_gizi = `Status Gizi`
) %>%
mutate(
jenis_kelamin = factor(jenis_kelamin),
status_gizi = factor(status_gizi)
)
# Cek missing value
colSums(is.na(data))## umur_bulan jenis_kelamin tinggi_badan status_gizi
## 0 0 0 0# Cek duplikasi
sum(duplicated(data))## [1] 82647# Ringkasan struktur
str(data)## tibble [120,999 × 4] (S3: tbl_df/tbl/data.frame)
## $ umur_bulan : num [1:120999] 0 0 0 0 0 0 0 0 0 0 ...
## $ jenis_kelamin: Factor w/ 2 levels "laki-laki","perempuan": 1 1 1 1 1 1 1 1 1 1 ...
## $ tinggi_badan : num [1:120999] 44.6 56.7 46.9 47.5 42.7 ...
## $ status_gizi : Factor w/ 4 levels "normal","severely stunted",..: 3 4 1 1 2 3 4 2 3 4 ...4 4. Statistik Deskriptif
summary(data)## umur_bulan jenis_kelamin tinggi_badan status_gizi
## Min. : 0.00 laki-laki:59997 Min. : 40.01 normal :67755
## 1st Qu.:15.00 perempuan:61002 1st Qu.: 77.00 severely stunted:19869
## Median :30.00 Median : 89.80 stunted :13815
## Mean :30.17 Mean : 88.66 tinggi :19560
## 3rd Qu.:45.00 3rd Qu.:101.20
## Max. :60.00 Max. :128.00# Statistik numerik
data %>%
summarise(
n = n(),
rata_rata_umur = mean(umur_bulan, na.rm = TRUE),
sd_umur = sd(umur_bulan, na.rm = TRUE),
min_umur = min(umur_bulan, na.rm = TRUE),
max_umur = max(umur_bulan, na.rm = TRUE),
rata_rata_tinggi = mean(tinggi_badan, na.rm = TRUE),
sd_tinggi = sd(tinggi_badan, na.rm = TRUE),
min_tinggi = min(tinggi_badan, na.rm = TRUE),
max_tinggi = max(tinggi_badan, na.rm = TRUE)
)## # A tibble: 1 × 9
## n rata_rata_umur sd_umur min_umur max_umur rata_rata_tinggi sd_tinggi
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 120999 30.2 17.6 0 60 88.7 17.3
## # ℹ 2 more variables: min_tinggi <dbl>, max_tinggi <dbl># Distribusi jenis kelamin
table(data$jenis_kelamin)##
## laki-laki perempuan
## 59997 61002# Distribusi status gizi
table(data$status_gizi)##
## normal severely stunted stunted tinggi
## 67755 19869 13815 19560# Persentase status gizi
prop.table(table(data$status_gizi)) * 100##
## normal severely stunted stunted tinggi
## 55.99633 16.42080 11.41745 16.165425 5. Eksplorasi Data
5.1 5.1 Hubungan Umur dan Tinggi Badan
ggplot(data, aes(x = umur_bulan, y = tinggi_badan)) +
geom_point(alpha = 0.15) +
geom_smooth(method = "lm", se = TRUE) +
labs(
title = "Hubungan Umur dengan Tinggi Badan Balita",
x = "Umur (bulan)",
y = "Tinggi Badan (cm)"
) +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
5.2 5.2 Tinggi Badan Berdasarkan Jenis Kelamin
ggplot(data, aes(x = jenis_kelamin, y = tinggi_badan)) +
geom_boxplot() +
labs(
title = "Distribusi Tinggi Badan Berdasarkan Jenis Kelamin",
x = "Jenis Kelamin",
y = "Tinggi Badan (cm)"
) +
theme_minimal()
5.3 5.3 Status Gizi
ggplot(data, aes(x = status_gizi)) +
geom_bar() +
labs(
title = "Distribusi Status Gizi Balita",
x = "Status Gizi",
y = "Jumlah Balita"
) +
theme_minimal()
6 6. Korelasi Variabel Numerik
Karena korelasi Pearson hanya digunakan untuk variabel numerik,
bagian ini menggunakan umur_bulan dan
tinggi_badan.
cor_data <- data %>%
select(umur_bulan, tinggi_badan)
cor_matrix <- cor(cor_data, use = "complete.obs", method = "pearson")
round(cor_matrix, 3)## umur_bulan tinggi_badan
## umur_bulan 1.000 0.843
## tinggi_badan 0.843 1.000corrplot(
cor_matrix,
method = "color",
type = "upper",
addCoef.col = "black",
tl.col = "black",
title = "Matriks Korelasi",
mar = c(0, 0, 2, 0)
)
7 7. Estimasi Multiple Linear Regression
Model yang digunakan:
\[ Y_i = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \varepsilon_i \]
dengan:
- \(Y\) = Tinggi Badan (cm)
- \(X_1\) = Umur (bulan)
- \(X_2\) = Jenis Kelamin
- \(\beta_0\) = intersep
- \(\beta_1,\beta_2\) = koefisien regresi
- \(\varepsilon\) = error
Jenis kelamin dimasukkan sebagai variabel faktor, sehingga R akan membentuk variabel indikator secara otomatis.
model <- lm(
tinggi_badan ~ umur_bulan + jenis_kelamin,
data = data
)
summary(model)##
## Call:
## lm(formula = tinggi_badan ~ umur_bulan + jenis_kelamin, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.8102 -7.4067 0.1104 7.7100 18.7888
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.000483 0.059756 1071.03 <2e-16 ***
## umur_bulan 0.829728 0.001521 545.62 <2e-16 ***
## jenis_kelaminperempuan -0.755908 0.053455 -14.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.296 on 120996 degrees of freedom
## Multiple R-squared: 0.7113, Adjusted R-squared: 0.7113
## F-statistic: 1.49e+05 on 2 and 120996 DF, p-value: < 2.2e-168 8. Koefisien Regresi
coef_table <- tidy(model, conf.int = TRUE, conf.level = 0.95)
kable(
coef_table,
digits = 4,
caption = "Koefisien Model Regresi Linear Berganda"
)| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 64.0005 | 0.0598 | 1071.0307 | 0 | 63.8834 | 64.1176 |
| umur_bulan | 0.8297 | 0.0015 | 545.6206 | 0 | 0.8267 | 0.8327 |
| jenis_kelaminperempuan | -0.7559 | 0.0535 | -14.1411 | 0 | -0.8607 | -0.6511 |
8.1 Persamaan Regresi
coef_model <- coef(model)
cat("Persamaan model:\n\n")## Persamaan model:cat(
"Tinggi_Badan = ",
round(coef_model[1], 4),
" + (",
round(coef_model[2], 4),
" × Umur_Bulan)"
)## Tinggi_Badan = 64.0005 + ( 0.8297 × Umur_Bulan)if ("jenis_kelaminperempuan" %in% names(coef_model)) {
cat(
" + (",
round(coef_model["jenis_kelaminperempuan"], 4),
" × Jenis_Kelamin_Perempuan)"
)
}## + ( -0.7559 × Jenis_Kelamin_Perempuan)cat("\n")Interpretasi koefisien dilakukan berdasarkan hasil
summary(model):
- Koefisien
umur_bulanmenunjukkan perubahan rata-rata tinggi badan untuk setiap tambahan 1 bulan umur, dengan jenis kelamin dianggap tetap. - Koefisien jenis kelamin menunjukkan perbedaan rata-rata tinggi badan antara kategori jenis kelamin yang dibandingkan dengan kategori referensi.
- Intersep adalah nilai prediksi ketika seluruh prediktor bernilai nol; interpretasi substantif intersep perlu hati-hati karena umur 0 bulan adalah batas khusus pada data.
9 9. Uji Signifikansi Simultan (Uji F)
Hipotesis:
- \(H_0\): seluruh koefisien slope sama dengan 0.
- \(H_1\): minimal terdapat satu koefisien slope yang tidak sama dengan 0.
anova_model <- anova(model)
anova_model## Analysis of Variance Table
##
## Response: tinggi_badan
## Df Sum Sq Mean Sq F value Pr(>F)
## umur_bulan 1 25743750 25743750 297886.52 < 2.2e-16 ***
## jenis_kelamin 1 17282 17282 199.97 < 2.2e-16 ***
## Residuals 120996 10456635 86
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1f_stat <- summary(model)$fstatistic
F_value <- unname(f_stat[1])
df1 <- unname(f_stat[2])
df2 <- unname(f_stat[3])
p_value_F <- pf(
F_value,
df1,
df2,
lower.tail = FALSE
)
cat("F-statistic =", round(F_value, 4), "\n")## F-statistic = 149043.2cat("df1 =", df1, "\n")## df1 = 2cat("df2 =", df2, "\n")## df2 = 120996cat("p-value =", format.pval(p_value_F, digits = 5), "\n")## p-value = < 2.22e-16if (p_value_F < 0.05) {
cat("Keputusan: Tolak H0 pada taraf signifikansi 5%.\n")
cat("Model secara simultan signifikan.\n")
} else {
cat("Keputusan: Gagal menolak H0 pada taraf signifikansi 5%.\n")
}## Keputusan: Tolak H0 pada taraf signifikansi 5%.
## Model secara simultan signifikan.10 10. Uji Signifikansi Parsial (Uji t)
Hipotesis untuk setiap koefisien:
- \(H_0: \beta_j = 0\)
- \(H_1: \beta_j \neq 0\)
uji_t <- tidy(model) %>%
mutate(
keputusan = ifelse(
p.value < 0.05,
"Signifikan",
"Tidak signifikan"
)
)
kable(
uji_t,
digits = 4,
caption = "Hasil Uji t Parsial"
)| term | estimate | std.error | statistic | p.value | keputusan |
|---|---|---|---|---|---|
| (Intercept) | 64.0005 | 0.0598 | 1071.0307 | 0 | Signifikan |
| umur_bulan | 0.8297 | 0.0015 | 545.6206 | 0 | Signifikan |
| jenis_kelaminperempuan | -0.7559 | 0.0535 | -14.1411 | 0 | Signifikan |
11 11. Koefisien Determinasi
model_summary <- summary(model)
R2 <- model_summary$r.squared
Adj_R2 <- model_summary$adj.r.squared
RSE <- model_summary$sigma
cat("R-squared =", round(R2, 4), "\n")## R-squared = 0.7113cat("Adjusted R-squared =", round(Adj_R2, 4), "\n")## Adjusted R-squared = 0.7113cat("Residual Standard Error =", round(RSE, 4), "\n")## Residual Standard Error = 9.2963cat(
"Persentase variasi tinggi badan yang dijelaskan model =",
round(R2 * 100, 2),
"%\n"
)## Persentase variasi tinggi badan yang dijelaskan model = 71.13 %12 12. Uji Asumsi Klasik
12.1 12.1 Normalitas Residual
Karena jumlah observasi sangat besar, uji Shapiro-Wilk tidak diterapkan langsung pada seluruh residual. Digunakan sampel residual maksimum 5000 observasi untuk uji tersebut.
res <- residuals(model)
set.seed(123)
if (length(res) > 5000) {
res_test <- sample(res, 5000)
} else {
res_test <- res
}
# Shapiro-Wilk
shapiro_result <- shapiro.test(res_test)
shapiro_result##
## Shapiro-Wilk normality test
##
## data: res_test
## W = 0.97759, p-value < 2.2e-16# Anderson-Darling
ad_result <- ad.test(res_test)
ad_result##
## Anderson-Darling normality test
##
## data: res_test
## A = 26.26, p-value < 2.2e-16# QQ Plot
qqnorm(res_test, main = "QQ Plot Residual")
qqline(res_test)
Kriteria:
- p-value > 0,05: tidak terdapat bukti yang cukup untuk menolak normalitas residual.
- p-value ≤ 0,05: terdapat bukti residual menyimpang dari normalitas.
Dengan jumlah data yang sangat besar, keputusan uji formal perlu dibaca bersama QQ-plot karena uji normalitas sangat sensitif terhadap penyimpangan kecil.
12.2 12.2 Uji Heteroskedastisitas
Digunakan Breusch-Pagan Test.
bp_test <- bptest(model)
bp_test##
## studentized Breusch-Pagan test
##
## data: model
## BP = 1205.3, df = 2, p-value < 2.2e-16cat("p-value =", format.pval(bp_test$p.value, digits = 5), "\n")## p-value = < 2.22e-16if (bp_test$p.value > 0.05) {
cat("Tidak terdapat bukti heteroskedastisitas pada taraf 5%.\n")
} else {
cat("Terdapat bukti heteroskedastisitas pada taraf 5%.\n")
}## Terdapat bukti heteroskedastisitas pada taraf 5%.12.3 12.3 Uji Multikolinearitas
vif_result <- vif(model)
vif_result## umur_bulan jenis_kelamin
## 1.000101 1.000101Karena model hanya memiliki dua prediktor, VIF digunakan untuk melihat apakah prediktor saling berkorelasi secara kuat.
Kriteria umum:
- VIF < 5: tidak ada indikasi multikolinearitas yang kuat.
- VIF ≥ 5: perlu pemeriksaan lebih lanjut.
12.4 12.4 Uji Autokorelasi
Karena data bukan merupakan data runtun waktu, autokorelasi bukan asumsi utama seperti pada regresi time series. Namun, jika ingin mengikuti format Hands On, Durbin-Watson dapat ditampilkan sebagai pemeriksaan tambahan.
dw_result <- dwtest(model)
dw_result##
## Durbin-Watson test
##
## data: model
## DW = 1.7742, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0cat("Durbin-Watson =", round(as.numeric(dw_result$statistic), 4), "\n")## Durbin-Watson = 1.7742cat("p-value =", format.pval(dw_result$p.value, digits = 5), "\n")## p-value = < 2.22e-1613 13. Diagnostik Residual
13.1 13.1 Residual vs Fitted
plot(
fitted(model),
residuals(model),
pch = 19,
cex = 0.4,
xlab = "Fitted Values",
ylab = "Residuals",
main = "Residual vs Fitted"
)
abline(h = 0, lty = 2)
Plot ini digunakan untuk melihat apakah residual menyebar secara acak di sekitar nol.
13.2 13.2 QQ Plot
qqnorm(
residuals(model),
pch = 19,
cex = 0.4,
main = "Normal Q-Q Plot Residual"
)
qqline(residuals(model))
13.3 13.3 Plot Diagnostik Bawaan
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))14 14. Diagnostik Outlier dan Observasi Berpengaruh
14.1 14.1 Cook’s Distance
cook <- cooks.distance(model)
cat("Nilai maksimum Cook's Distance =",
round(max(cook), 6), "\n")## Nilai maksimum Cook's Distance = 8.7e-05cat("Jumlah observasi dengan Cook's D > 1 =",
sum(cook > 1), "\n")## Jumlah observasi dengan Cook's D > 1 = 0plot(
cook,
type = "h",
main = "Cook's Distance",
xlab = "Observasi",
ylab = "Cook's Distance"
)
abline(h = 1, lty = 2)
14.2 14.2 Leverage
lev <- hatvalues(model)
threshold <- 2 * length(coef(model)) / nrow(data)
cat("Nilai leverage maksimum =",
round(max(lev), 6), "\n")## Nilai leverage maksimum = 4.1e-05cat("Threshold =",
round(threshold, 6), "\n")## Threshold = 5e-05cat(
"Jumlah observasi dengan leverage > threshold =",
sum(lev > threshold),
"\n"
)## Jumlah observasi dengan leverage > threshold = 014.3 14.3 Studentized Residual
student_resid <- rstudent(model)
cat(
"Jumlah observasi dengan |studentized residual| > 3 =",
sum(abs(student_resid) > 3),
"\n"
)## Jumlah observasi dengan |studentized residual| > 3 = 015 15. Prediksi
Sebagai contoh, prediksi tinggi badan berdasarkan umur dan jenis kelamin.
data_baru <- data.frame(
umur_bulan = c(12, 24, 36),
jenis_kelamin = factor(
c("laki-laki", "perempuan", "laki-laki"),
levels = levels(data$jenis_kelamin)
)
)
prediksi <- predict(
model,
newdata = data_baru,
interval = "confidence",
level = 0.95
)
hasil_prediksi <- cbind(
data_baru,
as.data.frame(prediksi)
)
kable(
hasil_prediksi,
digits = 3,
caption = "Contoh Prediksi Tinggi Badan"
)| umur_bulan | jenis_kelamin | fit | lwr | upr |
|---|---|---|---|---|
| 12 | laki-laki | 73.957 | 73.865 | 74.050 |
| 24 | perempuan | 83.158 | 83.082 | 83.234 |
| 36 | laki-laki | 93.871 | 93.794 | 93.947 |
16 16. Ringkasan Hasil Model
hasil_model <- data.frame(
R_squared = summary(model)$r.squared,
Adjusted_R_squared = summary(model)$adj.r.squared,
F_statistic = unname(summary(model)$fstatistic[1]),
p_value_F = p_value_F,
Residual_Standard_Error = summary(model)$sigma
)
kable(
hasil_model,
digits = 4,
caption = "Ringkasan Model MLR"
)| R_squared | Adjusted_R_squared | F_statistic | p_value_F | Residual_Standard_Error |
|---|---|---|---|---|
| 0.7113 | 0.7113 | 149043.2 | 0 | 9.2963 |
17 17. Kesimpulan
Berdasarkan analisis Multiple Linear Regression, model digunakan untuk menganalisis hubungan Umur (bulan) dan Jenis Kelamin terhadap Tinggi Badan (cm) pada data balita.
Kesimpulan akhir mengenai signifikansi model, signifikansi masing-masing prediktor, besarnya \(R^2\), serta terpenuhi atau tidaknya asumsi regresi harus dituliskan berdasarkan output yang dihasilkan setelah kode dijalankan di R/RStudio.
Secara umum, interpretasi dilakukan dengan melihat:
- Uji F untuk mengetahui apakah prediktor secara bersama-sama berhubungan dengan tinggi badan.
- Uji t untuk mengetahui pengaruh masing-masing prediktor.
- R-squared dan Adjusted R-squared untuk melihat kemampuan model menjelaskan variasi tinggi badan.
- Uji normalitas residual untuk memeriksa distribusi residual.
- Breusch-Pagan Test untuk memeriksa heteroskedastisitas.
- VIF untuk memeriksa multikolinearitas.
- Plot diagnostik dan Cook’s Distance untuk memeriksa observasi yang berpotensi berpengaruh.