# MEMUAT DATA
df <- read.csv("asuransi.csv", sep = ";", dec = ".")
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("1. DESKRIPSI DATA\n")
## 1. DESKRIPSI DATA
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
# Lihat struktur data (nama kolom, tipe data)
cat("Struktur Data\n")
## Struktur Data
str(df)
## 'data.frame': 1338 obs. of 7 variables:
## $ age : int 19 18 28 33 32 31 46 37 37 60 ...
## $ sex : chr "female" "male" "male" "male" ...
## $ bmi : num 27.9 33.8 33 22.7 28.9 ...
## $ children: int 0 1 3 0 0 0 1 3 2 0 ...
## $ smoker : chr "yes" "no" "no" "no" ...
## $ region : chr "southwest" "southeast" "southeast" "northwest" ...
## $ charges : num 16885 1726 4449 21984 3867 ...
cat("\n")
# Lihat 6 baris pertama data
cat("Data Awal\n")
## Data Awal
head(df)
## age sex bmi children smoker region charges
## 1 19 female 27.900 0 yes southwest 16884.924
## 2 18 male 33.770 1 no southeast 1725.552
## 3 28 male 33.000 3 no southeast 4449.462
## 4 33 male 22.705 0 no northwest 21984.471
## 5 32 male 28.880 0 no northwest 3866.855
## 6 31 female 25.740 0 no southeast 3756.622
cat("\n")
cat("Dimensi Data\n")
## Dimensi Data
cat("Jumlah baris:", nrow(df), "\n")
## Jumlah baris: 1338
cat("Jumlah kolom:", ncol(df), "\n")
## Jumlah kolom: 7
cat("\nCek Missing Value\n")
##
## Cek Missing Value
print(colSums(is.na(df)))
## age sex bmi children smoker region charges
## 0 0 0 0 0 0 0
cat("Total missing value:", sum(is.na(df)), "\n")
## Total missing value: 0
cat("\nCek Data Duplikat\n")
##
## Cek Data Duplikat
cat("Jumlah baris duplikat:", sum(duplicated(df)), "\n")
## Jumlah baris duplikat: 1
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("2. STATISTIKA DESKRIPTIF\n")
## 2. STATISTIKA DESKRIPTIF
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
summary(df[, c("age", "bmi", "children", "smoker", "charges")])
## age bmi children smoker
## Min. :18.00 Min. :15.96 Min. :0.000 Length:1338
## 1st Qu.:27.00 1st Qu.:26.30 1st Qu.:0.000 Class :character
## Median :39.00 Median :30.40 Median :1.000 Mode :character
## Mean :39.21 Mean :30.66 Mean :1.095
## 3rd Qu.:51.00 3rd Qu.:34.69 3rd Qu.:2.000
## Max. :64.00 Max. :53.13 Max. :5.000
## charges
## Min. : 1122
## 1st Qu.: 4740
## Median : 9382
## Mean :13270
## 3rd Qu.:16640
## Max. :63770
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("3. PRAPROSES DATA (MENGUBAH DAN MEMILIH VARIABEL)\n")
## 3. PRAPROSES DATA (MENGUBAH DAN MEMILIH VARIABEL)
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
# Mengubah variabel kategorik jadi Numerik
# Ubah smoker: yes = 1, no = 0
df$smoker <- ifelse(df$smoker == "yes", 1, 0)
# Cek hasil
head(df)
## age sex bmi children smoker region charges
## 1 19 female 27.900 0 1 southwest 16884.924
## 2 18 male 33.770 1 0 southeast 1725.552
## 3 28 male 33.000 3 0 southeast 4449.462
## 4 33 male 22.705 0 0 northwest 21984.471
## 5 32 male 28.880 0 0 northwest 3866.855
## 6 31 female 25.740 0 0 southeast 3756.622
str(df[, c("age", "bmi", "children", "smoker", "charges")])
## 'data.frame': 1338 obs. of 5 variables:
## $ age : int 19 18 28 33 32 31 46 37 37 60 ...
## $ bmi : num 27.9 33.8 33 22.7 28.9 ...
## $ children: int 0 1 3 0 0 0 1 3 2 0 ...
## $ smoker : num 1 0 0 0 0 0 0 0 0 0 ...
## $ charges : num 16885 1726 4449 21984 3867 ...
# Memilih Variabel independen (X)
X <- df[, c("age", "bmi", "children", "smoker")]
# Memilih Variabel dependen (Y)
y <- df$charges
# Cek 6 data pertama variabel independen
head(X)
## age bmi children smoker
## 1 19 27.900 0 1
## 2 18 33.770 1 0
## 3 28 33.000 3 0
## 4 33 22.705 0 0
## 5 32 28.880 0 0
## 6 31 25.740 0 0
# Cek 6 data pertama variabel dependen
cat("\n 6 Data pertama variabel Dependen (Charges) ")
##
## 6 Data pertama variabel Dependen (Charges)
head(y)
## [1] 16884.924 1725.552 4449.462 21984.471 3866.855 3756.622
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("4. EKSPLORASI DATA DAN MATRIKS KORELASI\n")
## 4. EKSPLORASI DATA DAN MATRIKS KORELASI
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("\nMatriks Korelasi:\n")
##
## Matriks Korelasi:
cor_matrix <- cor(df[, c("age", "bmi", "children", "smoker", "charges")])
print(round(cor_matrix, 3))
## age bmi children smoker charges
## age 1.000 0.109 0.042 -0.025 0.299
## bmi 0.109 1.000 0.013 0.004 0.198
## children 0.042 0.013 1.000 0.008 0.068
## smoker -0.025 0.004 0.008 1.000 0.787
## charges 0.299 0.198 0.068 0.787 1.000
# Visualisasi matriks korelasi
corrplot(cor_matrix, method = "color", type = "upper",
addCoef.col = "black", tl.col = "black", tl.srt = 45,
title = "Matriks Korelasi Antar Variabel", mar = c(0,0,2,0))
# Scatter plot matrix
pairs(df[, c("age", "bmi", "children", "smoker", "charges")], main = "Scatter Plot Matrix", pch = 19, col = "steelblue")
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("5. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION\n")
## 5. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
# Fit model MLR
model <- lm(charges ~ age + bmi + children + smoker, data = df)
# Ringkasan model
cat("\nRingkasan Model:\n")
##
## Ringkasan Model:
print(summary(model))
##
## Call:
## lm(formula = charges ~ age + bmi + children + smoker, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11897.9 -2920.8 -986.6 1392.2 29509.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -12102.77 941.98 -12.848 < 2e-16 ***
## age 257.85 11.90 21.675 < 2e-16 ***
## bmi 321.85 27.38 11.756 < 2e-16 ***
## children 473.50 137.79 3.436 0.000608 ***
## smoker 23811.40 411.22 57.904 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6068 on 1333 degrees of freedom
## Multiple R-squared: 0.7497, Adjusted R-squared: 0.7489
## F-statistic: 998.1 on 4 and 1333 DF, p-value: < 2.2e-16
# Ekstrak koefisien
cat("\nKoefisien Regresi:\n")
##
## Koefisien Regresi:
coef_df <- data.frame(
Variabel = names(coef(model)),
Koefisien = round(coef(model), 4),
Std_Error = round(summary(model)$coefficients[, 2], 4),
t_value = round(summary(model)$coefficients[, 3], 4),
p_value = round(summary(model)$coefficients[, 4], 4)
)
print(coef_df)
## Variabel Koefisien Std_Error t_value p_value
## (Intercept) (Intercept) -12102.7694 941.9839 -12.8482 0e+00
## age age 257.8495 11.8964 21.6746 0e+00
## bmi bmi 321.8514 27.3776 11.7560 0e+00
## children children 473.5023 137.7917 3.4364 6e-04
## smoker smoker 23811.3998 411.2197 57.9043 0e+00
# Interval kepercayaan 95%
cat("\nInterval Kepercayaan 95%:\n")
##
## Interval Kepercayaan 95%:
ci <- confint(model, level = 0.95)
print(round(ci, 4))
## 2.5 % 97.5 %
## (Intercept) -13950.7019 -10254.8369
## age 234.5118 281.1872
## bmi 268.1435 375.5593
## children 203.1902 743.8145
## smoker 23004.6915 24618.1082
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("6. UJI ASUMSI KLASIK\n")
## 6. UJI ASUMSI KLASIK
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("\n--- 6.1. UJI NORMALITAS RESIDUAL ---")
##
## --- 6.1. UJI NORMALITAS RESIDUAL ---
# Ambil residual dari model (nama variabel dibedakan dari fungsi residuals())
residuals_model <- residuals(model)
# Kolmogorov-Smirnov Test
ks_test <- ks.test(residuals_model, "pnorm", mean(residuals_model), sd(residuals_model))
## Warning in ks.test.default(residuals_model, "pnorm", mean(residuals_model), :
## ties should not be present for the one-sample Kolmogorov-Smirnov test
cat("\nKolmogorov-Smirnov Test:\n")
##
## Kolmogorov-Smirnov Test:
cat("D =", round(ks_test$statistic, 4), "\n")
## D = 0.1627
cat("p-value =", round(ks_test$p.value, 4), "\n")
## p-value = 0
cat("Kesimpulan:", ifelse(ks_test$p.value > 0.05,
"Residual berdistribusi normal",
"Residual TIDAK berdistribusi normal"), "\n")
## Kesimpulan: Residual TIDAK berdistribusi normal
# Visualisasi: Histogram residual
hist(residuals_model, main = "Histogram Residual", col = "lightblue",
xlab = "Residual", breaks = 30)
# Visualisasi: Q-Q Plot
qqnorm(residuals_model, main = "Q-Q Plot Residual")
qqline(residuals_model, col = "red", lwd = 2)
cat("\n--- 6.2. UJI HETEROSKEDASTISITAS ---")
##
## --- 6.2. UJI HETEROSKEDASTISITAS ---
# Uji Breusch-Pagan
bp_test <- bptest(model)
cat("\nUji Breusch-Pagan:\n")
##
## Uji Breusch-Pagan:
print(bp_test)
##
## studentized Breusch-Pagan test
##
## data: model
## BP = 117.16, df = 4, p-value < 2.2e-16
cat("\nKesimpulan:", ifelse(bp_test$p.value > 0.05,
"Tidak terjadi heteroskedastisitas (homoskedastisitas terpenuhi)",
"Terjadi heteroskedastisitas (varians residual tidak konstan)"), "\n")
##
## Kesimpulan: Terjadi heteroskedastisitas (varians residual tidak konstan)
# Visualisasi: Residual vs Fitted Values
plot(fitted(model), residuals_model,
main = "Residual vs Fitted Values",
xlab = "Fitted Values (Predicted Charges)",
ylab = "Residual",
pch = 19, col = "steelblue")
abline(h = 0, col = "red", lwd = 2, lty = 2)
# Scale-Location Plot (tambahan untuk heteroskedastisitas)
standardized_resid <- rstandard(model)
plot(fitted(model), sqrt(abs(standardized_resid)),
xlab = "Fitted Values", ylab = "√|Standardized Residuals|",
main = "Scale-Location Plot",
pch = 19, col = "steelblue")
lines(lowess(fitted(model), sqrt(abs(standardized_resid))),
col = "green", lwd = 2)
cat("\n--- 6.3. UJI MULTIKOLINEARITAS ---")
##
## --- 6.3. UJI MULTIKOLINEARITAS ---
# Hitung VIF
vif_values <- vif(model)
cat("\nNilai VIF:\n")
##
## Nilai VIF:
print(vif_values)
## age bmi children smoker
## 1.014498 1.012194 1.001950 1.000745
# Buat tabel rapi
vif_df <- data.frame(
Variabel = names(vif_values),
VIF = round(vif_values, 4))
print(vif_df)
## Variabel VIF
## age age 1.0145
## bmi bmi 1.0122
## children children 1.0019
## smoker smoker 1.0007
cat("\nInterpretasi VIF:\n")
##
## Interpretasi VIF:
for (i in 1:length(vif_values)) {
vif_val <- vif_values[i]
var_name <- names(vif_values)[i]
status <- ifelse(vif_val < 5, "Tidak ada multikolinearitas",
ifelse(vif_val < 10, "Multikolinearitas moderat",
"Multikolinearitas serius"))
cat(sprintf("%-10s: VIF = %7.4f → %s\n", var_name, vif_val, status))}
## age : VIF = 1.0145 → Tidak ada multikolinearitas
## bmi : VIF = 1.0122 → Tidak ada multikolinearitas
## children : VIF = 1.0019 → Tidak ada multikolinearitas
## smoker : VIF = 1.0007 → Tidak ada multikolinearitas
cat("\n--- 6.4. UJI LINEARITAS ---")
##
## --- 6.4. UJI LINEARITAS ---
# Ramsey RESET Test
reset_test <- resettest(model, power = 2:3, type = "fitted")
cat("\nRamsey RESET Test:\n")
##
## Ramsey RESET Test:
cat("F =", round(reset_test$statistic, 4), "\n")
## F = 71.2883
cat("p-value =", round(reset_test$p.value, 4), "\n")
## p-value = 0
cat("Kesimpulan:", ifelse(reset_test$p.value > 0.05,
"Model linear (spesifikasi benar)",
"Model TIDAK linear (spesifikasi salah)"), "\n")
## Kesimpulan: Model TIDAK linear (spesifikasi salah)
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("7. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS\n")
## 7. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
n <- nrow(df)
threshold_leverage <- 2 * length(coef(model)) / n
# Cook's Distance
cooks_d <- cooks.distance(model)
cat("\nCook's Distance:\n")
##
## Cook's Distance:
cat("Nilai maksimum:", round(max(cooks_d), 4), "\n")
## Nilai maksimum: 0.0304
cat("Jumlah observasi dengan Cook's D > 1:", sum(cooks_d > 1), "\n")
## Jumlah observasi dengan Cook's D > 1: 0
# Leverage (Hat Values)
leverage <- hatvalues(model)
cat("\nLeverage (Hat Values):\n")
##
## Leverage (Hat Values):
cat("Nilai maksimum:", round(max(leverage), 4), "\n")
## Nilai maksimum: 0.0151
cat("Threshold 2(k+1)/n:", round(threshold_leverage, 4), "\n")
## Threshold 2(k+1)/n: 0.0075
cat("Jumlah observasi dengan leverage > threshold:",
sum(leverage > threshold_leverage), "\n")
## Jumlah observasi dengan leverage > threshold: 52
# Studentized Residuals
std_resid <- rstudent(model)
cat("\nStudentized Residuals:\n")
##
## Studentized Residuals:
cat("Nilai maksimum absolut:", round(max(abs(std_resid)), 4), "\n")
## Nilai maksimum absolut: 4.9164
cat("Jumlah observasi dengan |std_resid| > 3:",
sum(abs(std_resid) > 3), "\n")
## Jumlah observasi dengan |std_resid| > 3: 27
# DFBETAS
dfbetas_val <- dfbetas(model)
cat("\nDFBETAS:\n")
##
## DFBETAS:
cat("Jumlah observasi dengan |DFBETAS| > 1:",
sum(abs(dfbetas_val) > 1), "\n")
## Jumlah observasi dengan |DFBETAS| > 1: 0
# Residual vs Leverage Plot
standardized_resid <- rstandard(model)
plot(leverage, standardized_resid,
xlab = "Leverage", ylab = "Standardized Residuals",
main = "Residual vs Leverage",
pch = 19, col = "steelblue")
abline(h = 0, col = "red", lty = 2)
abline(v = threshold_leverage, col = "red", lty = 2)
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("8. KOEFISIEN DETERMINASI (R² dan Adjusted R²)\n")
## 8. KOEFISIEN DETERMINASI (R² dan Adjusted R²)
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
r_squared <- summary(model)$r.squared
adj_r_squared <- summary(model)$adj.r.squared
resid_se <- summary(model)$sigma
cat("\nR-squared:", round(r_squared, 4), "\n")
##
## R-squared: 0.7497
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.7489
cat("Residual Standard Error:", round(resid_se, 4), "\n")
## Residual Standard Error: 6067.787
cat("Interpretasi: Model mampu menjelaskan",
round(r_squared * 100, 2), "% variasi pada Y\n")
## Interpretasi: Model mampu menjelaskan 74.97 % variasi pada Y
cat("Rata-rata, prediksi model meleset sekitar", round(resid_se, 2),
"dari nilai charges yang sebenarnya\n")
## Rata-rata, prediksi model meleset sekitar 6067.79 dari nilai charges yang sebenarnya
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("9. UJI SIGNIFIKANSI SIMULTAN (UJI F)\n")
## 9. UJI SIGNIFIKANSI SIMULTAN (UJI F)
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
# Ekstrak informasi F-test
f_stat <- summary(model)$fstatistic
f_value <- f_stat[1]
df1 <- f_stat[2]
df2 <- f_stat[3]
f_pvalue <- pf(f_value, df1, df2, lower.tail = FALSE)
cat("\nHipotesis:\n")
##
## Hipotesis:
cat("H0: β1 = β2 = β3 = β4 = 0 (Model tidak signifikan)\n")
## H0: β1 = β2 = β3 = β4 = 0 (Model tidak signifikan)
cat("H1: Minimal ada satu βj ≠ 0 (Model signifikan)\n")
## H1: Minimal ada satu βj ≠ 0 (Model signifikan)
cat("\nHasil Uji F:\n")
##
## Hasil Uji F:
cat("F-statistic:", round(f_value, 4), "\n")
## F-statistic: 998.1232
cat("df1:", df1, "\n")
## df1: 4
cat("df2:", df2, "\n")
## df2: 1333
cat("p-value:", format(f_pvalue, scientific = TRUE), "\n")
## p-value: 0e+00
cat("\nKesimpulan:\n")
##
## Kesimpulan:
if (f_pvalue < 0.05) {
cat("Tolak H0 → Model signifikan pada α = 5%\n")
cat("Artinya: Variabel independen (age, bmi, children, smoker) secara SIMULTAN",
"berpengaruh signifikan terhadap charges\n")
} else {cat("Gagal Tolak H0 → Model tidak signifikan\n")
cat("Artinya: Variabel independen secara simultan TIDAK berpengaruh signifikan terhadap charges\n")}
## Tolak H0 → Model signifikan pada α = 5%
## Artinya: Variabel independen (age, bmi, children, smoker) secara SIMULTAN berpengaruh signifikan terhadap charges
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("10. UJI SIGNIFIKANSI PARSIAL (UJI T)\n")
## 10. UJI SIGNIFIKANSI PARSIAL (UJI T)
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("\nHipotesis untuk setiap βj:\n")
##
## Hipotesis untuk setiap βj:
cat("H0: βj = 0 (Variabel tidak signifikan)\n")
## H0: βj = 0 (Variabel tidak signifikan)
cat("H1: βj ≠ 0 (Variabel signifikan)\n")
## H1: βj ≠ 0 (Variabel signifikan)
cat("\nHasil Uji T:\n")
##
## Hasil Uji T:
for (i in 2:nrow(coef_df)) {
var_name <- coef_df$Variabel[i]
t_val <- coef_df$t_value[i]
p_val <- coef_df$p_value[i]
sig <- ifelse(p_val < 0.001, "***",
ifelse(p_val < 0.01, "**",
ifelse(p_val < 0.05, "*",
ifelse(p_val < 0.1, ".", "ns"))))
cat(sprintf("%-10s: t = %7.4f, p = %8.4f %s\n",
var_name, t_val, p_val, sig))
}
## age : t = 21.6746, p = 0.0000 ***
## bmi : t = 11.7560, p = 0.0000 ***
## children : t = 3.4364, p = 0.0006 ***
## smoker : t = 57.9043, p = 0.0000 ***
cat("\nKeterangan: *** p<0.001, ** p<0.01, * p<0.05, . p<0.1, ns tidak signifikan\n")
##
## Keterangan: *** p<0.001, ** p<0.01, * p<0.05, . p<0.1, ns tidak signifikan
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("11. RINGKASAN HASIL ANALISIS\n")
## 11. RINGKASAN HASIL ANALISIS
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("\n--- MODEL AKHIR ---\n")
##
## --- MODEL AKHIR ---
cat("Persamaan Regresi:\n")
## Persamaan Regresi:
cat(sprintf("Charges = %.4f + %.4f*Age + %.4f*BMI + %.4f*Children + %.4f*Smoker\n",
coef(model)[1], coef(model)[2], coef(model)[3],
coef(model)[4], coef(model)[5]))
## Charges = -12102.7694 + 257.8495*Age + 321.8514*BMI + 473.5023*Children + 23811.3998*Smoker
cat("\n--- UJI SIGNIFIKANSI ---\n")
##
## --- UJI SIGNIFIKANSI ---
cat("Uji F (Simultan): p-value =", format(f_pvalue, scientific = TRUE), "\n")
## Uji F (Simultan): p-value = 0e+00
cat("Uji t (Parsial):\n")
## Uji t (Parsial):
for (i in 2:nrow(coef_df)) {
cat(sprintf(" %s: p-value = %.4f\n",
coef_df$Variabel[i], coef_df$p_value[i]))
}
## age: p-value = 0.0000
## bmi: p-value = 0.0000
## children: p-value = 0.0006
## smoker: p-value = 0.0000
cat("\n--- UJI ASUMSI KLASIK ---\n")
##
## --- UJI ASUMSI KLASIK ---
cat("Normalitas (Kolmogorov-Smirnov): p-value =",
round(ks_test$p.value, 4), "\n")
## Normalitas (Kolmogorov-Smirnov): p-value = 0
cat("Heteroskedastisitas (Breusch-Pagan): p-value =",
round(bp_test$p.value, 4), "\n")
## Heteroskedastisitas (Breusch-Pagan): p-value = 0
cat("Multikolinearitas (VIF max):",
round(max(vif_values), 4), "\n")
## Multikolinearitas (VIF max): 1.0145
cat("\n--- KESIMPULAN ---\n")
##
## --- KESIMPULAN ---
cat("R-squared:", round(r_squared, 4), "\n")
## R-squared: 0.7497
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.7489
cat("===== ANALISIS SELESAI =====\n")
## ===== ANALISIS SELESAI =====