# ============================================================
# LOAD LIBRARY
# ============================================================
library(car) # VIF dan diagnostik
library(lmtest) # Breusch-Pagan, Durbin-Watson, RESET
library(nortest) # Anderson-Darling
library(MASS) # Box-Cox
library(ggplot2) # Visualisasi
library(corrplot) # Matriks korelasi
# ============================================================
# 1. MEMBACA DATASET
# ============================================================
# Membaca dataset
data <- read.csv("C:/Users/Acer/Downloads/company_data.csv")
# Melihat struktur data
str(data)
## 'data.frame': 167 obs. of 11 variables:
## $ X : int 0 1 2 3 4 5 6 7 8 9 ...
## $ company : chr "Company_1" "Company_2" "Company_3" "Company_4" ...
## $ employee_turnover: num 90.2 16.6 27.3 119 10.3 14.5 18.1 4.8 4.3 39.2 ...
## $ revenue_growth : num 10 28 38.4 62.3 45.5 18.9 20.8 19.8 51.3 54.3 ...
## $ rd_investment : num 7.58 6.55 4.17 2.85 6.03 8.1 4.4 8.73 11 5.88 ...
## $ operational_cost : num 44.9 48.6 31.4 42.9 58.9 16 45.3 20.9 47.8 20.7 ...
## $ average_salary : int 1610 9930 12900 5900 19100 18700 6700 41400 43200 16000 ...
## $ market_volatility: num 9.44 4.49 16.1 22.4 1.44 20.9 7.77 1.16 0.873 13.8 ...
## $ average_tenure : num 56.2 76.3 76.5 60.1 76.8 75.8 73.3 82 80.5 69.1 ...
## $ growth_potential : num 5.82 1.65 2.89 6.16 2.13 2.37 1.69 1.93 1.44 1.92 ...
## $ net_profit : int 553 4090 4460 3530 12200 10300 3220 51900 46900 5840 ...
# Melihat nama kolom
names(data)
## [1] "X" "company" "employee_turnover"
## [4] "revenue_growth" "rd_investment" "operational_cost"
## [7] "average_salary" "market_volatility" "average_tenure"
## [10] "growth_potential" "net_profit"
# Melihat ukuran dataset
dim(data)
## [1] 167 11
# Melihat 6 baris pertama
head(data)
## X company employee_turnover revenue_growth rd_investment operational_cost
## 1 0 Company_1 90.2 10.0 7.58 44.9
## 2 1 Company_2 16.6 28.0 6.55 48.6
## 3 2 Company_3 27.3 38.4 4.17 31.4
## 4 3 Company_4 119.0 62.3 2.85 42.9
## 5 4 Company_5 10.3 45.5 6.03 58.9
## 6 5 Company_6 14.5 18.9 8.10 16.0
## average_salary market_volatility average_tenure growth_potential net_profit
## 1 1610 9.44 56.2 5.82 553
## 2 9930 4.49 76.3 1.65 4090
## 3 12900 16.10 76.5 2.89 4460
## 4 5900 22.40 60.1 6.16 3530
## 5 19100 1.44 76.8 2.13 12200
## 6 18700 20.90 75.8 2.37 10300
# ------------------------------------------------------------
# 2. STATISTIK DESKRIPTIF DAN EKSPLORASI DATA
# ------------------------------------------------------------
# Memilih variabel penelitian
data_model <- data[, c(
"net_profit",
"revenue_growth",
"rd_investment",
"operational_cost",
"average_salary",
"market_volatility"
)]
cat("\n", rep("=", 60), "\n", sep = "")
##
## ============================================================
cat("1. STATISTIK DESKRIPTIF\n")
## 1. STATISTIK DESKRIPTIF
cat(rep("=", 60), "\n", sep = "")
## ============================================================
# Statistik deskriptif
cat("\nSTATISTIK DESKRIPTIF:\n")
##
## STATISTIK DESKRIPTIF:
print(summary(data_model))
## net_profit revenue_growth rd_investment operational_cost
## Min. : 231 Min. : 0.109 Min. : 1.810 Min. : 0.0659
## 1st Qu.: 1330 1st Qu.: 23.800 1st Qu.: 4.920 1st Qu.: 30.2000
## Median : 4660 Median : 35.000 Median : 6.320 Median : 43.3000
## Mean : 12964 Mean : 41.109 Mean : 6.816 Mean : 46.8902
## 3rd Qu.: 14050 3rd Qu.: 51.350 3rd Qu.: 8.600 3rd Qu.: 58.7500
## Max. :105000 Max. :200.000 Max. :17.900 Max. :174.0000
## average_salary market_volatility
## Min. : 609 Min. : -4.210
## 1st Qu.: 3355 1st Qu.: 1.810
## Median : 9960 Median : 5.390
## Mean : 17145 Mean : 7.782
## 3rd Qu.: 22800 3rd Qu.: 10.750
## Max. :125000 Max. :104.000
# Matriks korelasi
cat("\nMATRIKS KORELASI:\n")
##
## MATRIKS KORELASI:
cor_matrix <- cor(data_model, use = "complete.obs")
print(round(cor_matrix, 3))
## net_profit revenue_growth rd_investment operational_cost
## net_profit 1.000 0.419 0.346 0.115
## revenue_growth 0.419 1.000 -0.114 0.737
## rd_investment 0.346 -0.114 1.000 0.096
## operational_cost 0.115 0.737 0.096 1.000
## average_salary 0.896 0.517 0.130 0.122
## market_volatility -0.222 -0.107 -0.255 -0.247
## average_salary market_volatility
## net_profit 0.896 -0.222
## revenue_growth 0.517 -0.107
## rd_investment 0.130 -0.255
## operational_cost 0.122 -0.247
## average_salary 1.000 -0.148
## market_volatility -0.148 1.000
# Visualisasi matriks korelasi
corrplot::corrplot(
cor_matrix,
method = "color",
type = "upper",
addCoef.col = "black",
tl.col = "black",
tl.srt = 45,
title = "Matriks Korelasi Antarvariabel",
mar = c(0, 0, 2, 0)
)
# Scatter plot matrix
pairs(
data_model,
main = "Scatter Plot Matrix",
pch = 19,
col = "steelblue"
)
# ============================================================
# 3. ESTIMASI MODEL REGRESI LINEAR BERGANDA
# ============================================================
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("2. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION\n")
## 2. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION
cat(strrep("=", 60), "\n")
## ============================================================
# Membentuk model regresi
model <- lm(
net_profit ~ revenue_growth +
rd_investment +
operational_cost +
average_salary +
market_volatility,
data = data_model
)
# Ringkasan model
cat("\nRINGKASAN MODEL:\n")
##
## RINGKASAN MODEL:
print(summary(model))
##
## Call:
## lm(formula = net_profit ~ revenue_growth + rd_investment + operational_cost +
## average_salary + market_volatility, data = data_model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23022 -4172 164 3188 33683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.017e+04 2.037e+03 -4.993 1.53e-06 ***
## revenue_growth 3.337e+01 4.349e+01 0.767 0.444
## rd_investment 1.579e+03 2.276e+02 6.940 9.10e-11 ***
## operational_cost -4.399e+01 4.250e+01 -1.035 0.302
## average_salary 7.981e-01 4.089e-02 19.517 < 2e-16 ***
## market_volatility -8.003e+01 5.579e+01 -1.434 0.153
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7012 on 161 degrees of freedom
## Multiple R-squared: 0.858, Adjusted R-squared: 0.8536
## F-statistic: 194.6 on 5 and 161 DF, p-value: < 2.2e-16
# Tabel koefisien
cat("\nKOEFISIEN REGRESI:\n")
##
## KOEFISIEN REGRESI:
hasil_koefisien <- summary(model)$coefficients
coef_df <- data.frame(
Variabel = rownames(hasil_koefisien),
Koefisien = hasil_koefisien[, 1],
Std_Error = hasil_koefisien[, 2],
t_value = hasil_koefisien[, 3],
p_value = hasil_koefisien[, 4],
row.names = NULL
)
# Membulatkan kolom numerik saja
coef_df[, -1] <- round(coef_df[, -1], 4)
# Menampilkan tabel
print(coef_df, row.names = FALSE)
## Variabel Koefisien Std_Error t_value p_value
## (Intercept) -10168.5967 2036.6481 -4.9928 0.0000
## revenue_growth 33.3675 43.4867 0.7673 0.4440
## rd_investment 1579.2823 227.5674 6.9398 0.0000
## operational_cost -43.9892 42.5001 -1.0350 0.3022
## average_salary 0.7981 0.0409 19.5174 0.0000
## market_volatility -80.0348 55.7939 -1.4345 0.1534
# ------------------------------------------------------------
# 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) -14190.5859 -6146.6075
## revenue_growth -52.5103 119.2454
## rd_investment 1129.8803 2028.6843
## operational_cost -127.9188 39.9404
## average_salary 0.7173 0.8788
## market_volatility -190.2170 30.1474
# ============================================================
# 4. UJI SIGNIFIKANSI SIMULTAN (UJI F)
# ============================================================
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("3. UJI SIGNIFIKANSI SIMULTAN (UJI F)\n")
## 3. UJI SIGNIFIKANSI SIMULTAN (UJI F)
cat(strrep("=", 60), "\n")
## ============================================================
cat("\nHIPOTESIS:\n")
##
## HIPOTESIS:
cat("H0: β1 = β2 = ... = β5 = 0\n")
## H0: β1 = β2 = ... = β5 = 0
cat("H1: Minimal terdapat satu βj ≠0\n")
## H1: Minimal terdapat satu βj ≠0
cat("Taraf signifikansi: α = 0.05\n")
## Taraf signifikansi: α = 0.05
# Mengambil informasi F-statistic
f_stat <- summary(model)$fstatistic
f_value <- unname(f_stat[1])
df1 <- unname(f_stat[2])
df2 <- unname(f_stat[3])
f_pvalue <- pf(f_value, df1, df2, lower.tail = FALSE)
cat("\nHASIL UJI F:\n")
##
## HASIL UJI F:
cat("F-statistic:", round(f_value, 4), "\n")
## F-statistic: 194.6336
cat("df1:", df1, "\n")
## df1: 5
cat("df2:", df2, "\n")
## df2: 161
cat("p-value:", format.pval(f_pvalue, digits = 5), "\n")
## p-value: < 2.22e-16
if (f_pvalue < 0.05) {
cat("Kesimpulan: Tolak H0.\n")
cat("Model signifikan secara simultan pada α = 5%.\n")
} else {
cat("Kesimpulan: Gagal menolak H0.\n")
cat("Model tidak signifikan secara simultan pada α = 5%.\n")
}
## Kesimpulan: Tolak H0.
## Model signifikan secara simultan pada α = 5%.
# ============================================================
# 5. UJI SIGNIFIKANSI PARSIAL (UJI T)
# ============================================================
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("4. UJI SIGNIFIKANSI PARSIAL (UJI T)\n")
## 4. UJI SIGNIFIKANSI PARSIAL (UJI T)
cat(strrep("=", 60), "\n")
## ============================================================
cat("\nHIPOTESIS UNTUK SETIAP KOEFISIEN:\n")
##
## HIPOTESIS UNTUK SETIAP KOEFISIEN:
cat("H0: βj = 0\n")
## H0: βj = 0
cat("H1: βj ≠0\n")
## H1: βj ≠0
cat("Taraf signifikansi: α = 0.05\n")
## Taraf signifikansi: α = 0.05
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]
keputusan <- ifelse(
p_val < 0.05,
"Tolak H0",
"Gagal menolak H0"
)
cat(sprintf("%-22s: t = %8.4f, p-value = %10.5f, %s\n",
var_name, t_val, p_val, keputusan))
}
## revenue_growth : t = 0.7673, p-value = 0.44400, Gagal menolak H0
## rd_investment : t = 6.9398, p-value = 0.00000, Tolak H0
## operational_cost : t = -1.0350, p-value = 0.30220, Gagal menolak H0
## average_salary : t = 19.5174, p-value = 0.00000, Tolak H0
## market_volatility : t = -1.4345, p-value = 0.15340, Gagal menolak H0
# ============================================================
# 6. KOEFISIEN DETERMINASI
# ============================================================
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("5. KOEFISIEN DETERMINASI\n")
## 5. KOEFISIEN DETERMINASI
cat(strrep("=", 60), "\n")
## ============================================================
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.858
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.8536
cat("Residual Standard Error:", round(resid_se, 4), "\n")
## Residual Standard Error: 7012.083
cat(
"Interpretasi awal: Model menjelaskan sekitar",
round(r_squared * 100, 2),
"% variasi pada net_profit.\n"
)
## Interpretasi awal: Model menjelaskan sekitar 85.8 % variasi pada net_profit.
# ============================================================
# 7. KRITERIA PEMILIHAN MODEL
# AIC, BIC, DAN AICC
# ============================================================
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("6. KRITERIA PEMILIHAN MODEL\n")
## 6. KRITERIA PEMILIHAN MODEL
cat(strrep("=", 60), "\n")
## ============================================================
# AIC
aic_value <- AIC(model)
# BIC
bic_value <- BIC(model)
# AICc
n_model <- nobs(model)
k <- length(coef(model))
aicc_value <- AIC(model) +
(2 * k * (k + 1)) / (n_model - k - 1)
cat("\nAIC :", round(aic_value, 4), "\n")
##
## AIC : 3439.515
cat("BIC :", round(bic_value, 4), "\n")
## BIC : 3461.341
cat("AICc:", round(aicc_value, 4), "\n")
## AICc: 3440.04
cat(
"\nCatatan: Nilai AIC, BIC, dan AICc digunakan untuk membandingkan model.\n",
"Nilai yang lebih kecil menunjukkan model yang lebih baik menurut kriteria tersebut,\n",
"dengan mempertimbangkan kecocokan model dan kompleksitasnya.\n"
)
##
## Catatan: Nilai AIC, BIC, dan AICc digunakan untuk membandingkan model.
## Nilai yang lebih kecil menunjukkan model yang lebih baik menurut kriteria tersebut,
## dengan mempertimbangkan kecocokan model dan kompleksitasnya.
# ============================================================
# 8. UJI ASUMSI RESIDUAL
# ============================================================
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("7. UJI ASUMSI RESIDUAL\n")
## 7. UJI ASUMSI RESIDUAL
cat(strrep("=", 60), "\n")
## ============================================================
# Mengambil komponen model
residual_model <- residuals(model)
fitted_values <- fitted(model)
standardized_resid <- rstandard(model)
# ------------------------------------------------------------
# 8.1. UJI NORMALITAS RESIDUAL
# ------------------------------------------------------------
cat("\n--- 7.1. UJI NORMALITAS RESIDUAL ---\n")
##
## --- 7.1. UJI NORMALITAS RESIDUAL ---
# Shapiro-Wilk
shapiro_test <- shapiro.test(residual_model)
cat("\nSHAPIRO-WILK TEST:\n")
##
## SHAPIRO-WILK TEST:
cat("W =", round(unname(shapiro_test$statistic), 4), "\n")
## W = 0.8982
cat("p-value =", format.pval(shapiro_test$p.value, digits = 5), "\n")
## p-value = 2.5201e-09
if (shapiro_test$p.value > 0.05) {
cat("Kesimpulan: Tidak terdapat bukti yang cukup untuk menyatakan residual tidak normal.\n")
} else {
cat("Kesimpulan: Terdapat bukti bahwa residual tidak berdistribusi normal.\n")
}
## Kesimpulan: Terdapat bukti bahwa residual tidak berdistribusi normal.
# Anderson-Darling
ad_test <- ad.test(residual_model)
cat("\nANDERSON-DARLING TEST:\n")
##
## ANDERSON-DARLING TEST:
cat("A =", round(unname(ad_test$statistic), 4), "\n")
## A = 2.9991
cat("p-value =", format.pval(ad_test$p.value, digits = 5), "\n")
## p-value = 1.4629e-07
# ------------------------------------------------------------
# 8.2. UJI HETEROSKEDASTISITAS
# ------------------------------------------------------------
cat("\n--- 7.2. UJI HETEROSKEDASTISITAS ---\n")
##
## --- 7.2. UJI HETEROSKEDASTISITAS ---
bp_test <- bptest(model)
cat("\nBREUSCH-PAGAN TEST:\n")
##
## BREUSCH-PAGAN TEST:
cat("BP =", round(unname(bp_test$statistic), 4), "\n")
## BP = 72.1064
cat("p-value =", format.pval(bp_test$p.value, digits = 5), "\n")
## p-value = 3.7327e-14
if (bp_test$p.value > 0.05) {
cat("Kesimpulan: Tidak terdapat bukti heteroskedastisitas.\n")
} else {
cat("Kesimpulan: Terdapat indikasi heteroskedastisitas.\n")
}
## Kesimpulan: Terdapat indikasi heteroskedastisitas.
# ------------------------------------------------------------
# 8.3. UJI AUTOKORELASI
# ------------------------------------------------------------
cat("\n--- 7.3. UJI AUTOKORELASI ---\n")
##
## --- 7.3. UJI AUTOKORELASI ---
dw_test <- dwtest(model)
cat("\nDURBIN-WATSON TEST:\n")
##
## DURBIN-WATSON TEST:
cat("DW =", round(unname(dw_test$statistic), 4), "\n")
## DW = 1.8828
cat("p-value =", format.pval(dw_test$p.value, digits = 5), "\n")
## p-value = 0.21739
if (dw_test$p.value > 0.05) {
cat("Kesimpulan: Tidak terdapat bukti autokorelasi.\n")
} else {
cat("Kesimpulan: Terdapat indikasi autokorelasi berdasarkan uji Durbin-Watson.\n")
}
## Kesimpulan: Tidak terdapat bukti autokorelasi.
# ------------------------------------------------------------
# 8.4. UJI MULTIKOLINEARITAS
# ------------------------------------------------------------
cat("\n--- 7.4. UJI MULTIKOLINEARITAS ---\n")
##
## --- 7.4. UJI MULTIKOLINEARITAS ---
vif_values <- vif(model)
cat("\nVIF:\n")
##
## VIF:
print(round(vif_values, 4))
## revenue_growth rd_investment operational_cost average_salary
## 4.7974 1.3192 3.5741 2.0979
## market_volatility
## 1.1743
cat("\nINTERPRETASI VIF:\n")
##
## INTERPRETASI VIF:
for (i in seq_along(vif_values)) {
vif_val <- vif_values[i]
var_name <- names(vif_values)[i]
status <- ifelse(
vif_val < 5,
"Tidak ada indikasi multikolinearitas serius",
ifelse(
vif_val < 10,
"Indikasi multikolinearitas moderat",
"Indikasi multikolinearitas tinggi"))
cat(sprintf("%-22s: VIF = %8.4f -> %s\n",
var_name, vif_val, status))
}
## revenue_growth : VIF = 4.7974 -> Tidak ada indikasi multikolinearitas serius
## rd_investment : VIF = 1.3192 -> Tidak ada indikasi multikolinearitas serius
## operational_cost : VIF = 3.5741 -> Tidak ada indikasi multikolinearitas serius
## average_salary : VIF = 2.0979 -> Tidak ada indikasi multikolinearitas serius
## market_volatility : VIF = 1.1743 -> Tidak ada indikasi multikolinearitas serius
# Tolerance
tolerance <- 1 / vif_values
cat("\nTOLERANCE:\n")
##
## TOLERANCE:
print(round(tolerance, 4))
## revenue_growth rd_investment operational_cost average_salary
## 0.2084 0.7581 0.2798 0.4767
## market_volatility
## 0.8515
# Condition number
condition_number <- kappa(model.matrix(model))
cat("\nCONDITION NUMBER:\n")
##
## CONDITION NUMBER:
cat(round(condition_number, 4), "\n")
## 107503.7
# ------------------------------------------------------------
# 8.5. UJI SPESIFIKASI MODEL
# ------------------------------------------------------------
cat("\n--- 7.5. UJI SPESIFIKASI MODEL ---\n")
##
## --- 7.5. UJI SPESIFIKASI MODEL ---
reset_test <- resettest(
model,
power = 2:3,
type = "fitted")
cat("\nRAMSEY RESET TEST:\n")
##
## RAMSEY RESET TEST:
cat("F =", round(unname(reset_test$statistic), 4), "\n")
## F = 35.1696
cat("p-value =", format.pval(reset_test$p.value, digits = 5), "\n")
## p-value = 2.2545e-13
if (reset_test$p.value > 0.05) {
cat("Kesimpulan: Tidak terdapat bukti yang cukup mengenai kesalahan spesifikasi model.\n")
} else {
cat("Kesimpulan: Terdapat indikasi kesalahan spesifikasi model.\n")
}
## Kesimpulan: Terdapat indikasi kesalahan spesifikasi model.
par(mfrow = c(1, 2))
hist(
residuals(model),
main = "Histogram Residual",
xlab = "Residual")
qqnorm(residuals(model))
qqline(residuals(model), col = "red")
par(mfrow = c(1, 1))
plot(
fitted(model),
residuals(model),
main = "Residual vs Fitted",
xlab = "Nilai Prediksi",
ylab = "Residual")
abline(h = 0, col = "red", lty = 2)
# ============================================================
# 9. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
# ============================================================
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("8. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS\n")
## 8. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
cat(strrep("=", 60), "\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: 1.186
cat(
"Jumlah observasi dengan Cook's D > 1:",
sum(cooks_d > 1),
"\n")
## Jumlah observasi dengan Cook's D > 1: 1
# Leverage
leverage <- hatvalues(model)
p_parameter <- length(coef(model))
leverage_threshold <- 2 * p_parameter / n_model
cat("\nLEVERAGE:\n")
##
## LEVERAGE:
cat("Nilai maksimum:", round(max(leverage), 4), "\n")
## Nilai maksimum: 0.5336
cat("Threshold:", round(leverage_threshold, 4), "\n")
## Threshold: 0.0719
cat("Jumlah observasi dengan leverage > threshold:", sum(leverage > leverage_threshold), "\n")
## Jumlah observasi dengan leverage > threshold: 13
# Studentized residual
studentized_resid <- rstudent(model)
cat("\nSTUDENTIZED RESIDUAL:\n")
##
## STUDENTIZED RESIDUAL:
cat("Nilai maksimum absolut:", round(max(abs(studentized_resid)), 4), "\n")
## Nilai maksimum absolut: 5.3424
cat("Jumlah observasi dengan |studentized residual| > 3:", sum(abs(studentized_resid) > 3), "\n")
## Jumlah observasi dengan |studentized residual| > 3: 5
# 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: 5
# ============================================================
# 10. VISUALISASI DIAGNOSTIK
# ============================================================
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("9. VISUALISASI DIAGNOSTIK\n")
## 9. VISUALISASI DIAGNOSTIK
cat(strrep("=", 60), "\n")
## ============================================================
par(mfrow = c(2, 2))
# Plot 1: Residual vs Fitted
plot(
fitted_values,
residual_model,
xlab = "Fitted Values",
ylab = "Residuals",
main = "Residual vs Fitted",
pch = 19, col = "steelblue")
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted_values, residual_model), col = "green", lwd = 2)
# Plot 2: Normal Q-Q
qqnorm(residual_model, pch = 19, col = "steelblue", main = "Normal Q-Q Plot")
qqline(residual_model, col = "red", lwd = 2)
# Plot 3: Scale-Location
plot(
fitted_values,
sqrt(abs(standardized_resid)),
xlab = "Fitted Values",
ylab = "√|Standardized Residuals|",
main = "Scale-Location",
pch = 19, col = "steelblue")
lines(
lowess(
fitted_values,
sqrt(abs(standardized_resid))),
col = "green", lwd = 2)
# Plot 4: Residual vs Leverage
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 = leverage_threshold, col = "red",
lty = 2)
par(mfrow = c(1, 1))
# ============================================================
# 11. UJI BOX-COX
# ============================================================
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("10. UJI BOX-COX\n")
## 10. UJI BOX-COX
cat(strrep("=", 60), "\n")
## ============================================================
# Box-Cox memerlukan respons bernilai positif
if (all(model.response(model.frame(model)) > 0)) {
bc <- boxcox(
model,
lambda = seq(-2, 2, by = 0.1)
)
lambda_opt <- bc$x[which.max(bc$y)]
cat("\nLambda optimal:", round(lambda_opt, 4), "\n")
cat("\nInterpretasi umum:\n")
cat("Lambda mendekati 1 : transformasi mungkin tidak diperlukan.\n")
cat("Lambda mendekati 0 : transformasi log dapat dipertimbangkan.\n")
cat("Lambda mendekati 0.5: transformasi akar kuadrat dapat dipertimbangkan.\n")
} else {
cat("\nBox-Cox tidak dijalankan karena terdapat nilai Y <= 0.\n")
}
##
## Lambda optimal: 0.3838
##
## Interpretasi umum:
## Lambda mendekati 1 : transformasi mungkin tidak diperlukan.
## Lambda mendekati 0 : transformasi log dapat dipertimbangkan.
## Lambda mendekati 0.5: transformasi akar kuadrat dapat dipertimbangkan.
# ============================================================
# 12. PEMILIHAN MODEL STEPWISE
# ============================================================
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("11. PEMILIHAN MODEL STEPWISE\n")
## 11. PEMILIHAN MODEL STEPWISE
cat(strrep("=", 60), "\n")
## ============================================================
# ------------------------------------------------------------
# 12.1. STEPWISE BOTH
# ------------------------------------------------------------
cat("\n--- STEPWISE BOTH ---\n")
##
## --- STEPWISE BOTH ---
step_model <- step(
model,
direction = "both",
trace = 1)
## Start: AIC=2963.59
## net_profit ~ revenue_growth + rd_investment + operational_cost +
## average_salary + market_volatility
##
## Df Sum of Sq RSS AIC
## - revenue_growth 1 2.8949e+07 7.9452e+09 2962.2
## - operational_cost 1 5.2675e+07 7.9689e+09 2962.7
## <none> 7.9163e+09 2963.6
## - market_volatility 1 1.0118e+08 8.0174e+09 2963.7
## - rd_investment 1 2.3681e+09 1.0284e+10 3005.3
## - average_salary 1 1.8730e+10 2.6646e+10 3164.3
##
## Step: AIC=2962.2
## net_profit ~ rd_investment + operational_cost + average_salary +
## market_volatility
##
## Df Sum of Sq RSS AIC
## - operational_cost 1 2.5307e+07 7.9705e+09 2960.7
## - market_volatility 1 8.6198e+07 8.0314e+09 2962.0
## <none> 7.9452e+09 2962.2
## + revenue_growth 1 2.8949e+07 7.9163e+09 2963.6
## - rd_investment 1 2.6256e+09 1.0571e+10 3007.9
## - average_salary 1 3.9934e+10 4.7879e+10 3260.2
##
## Step: AIC=2960.73
## net_profit ~ rd_investment + average_salary + market_volatility
##
## Df Sum of Sq RSS AIC
## - market_volatility 1 7.0227e+07 8.0407e+09 2960.2
## <none> 7.9705e+09 2960.7
## + operational_cost 1 2.5307e+07 7.9452e+09 2962.2
## + revenue_growth 1 1.5807e+06 7.9689e+09 2962.7
## - rd_investment 1 2.6138e+09 1.0584e+10 3006.1
## - average_salary 1 4.0061e+10 4.8032e+10 3258.7
##
## Step: AIC=2960.2
## net_profit ~ rd_investment + average_salary
##
## Df Sum of Sq RSS AIC
## <none> 8.0407e+09 2960.2
## + market_volatility 1 7.0227e+07 7.9705e+09 2960.7
## + operational_cost 1 9.3360e+06 8.0314e+09 2962.0
## + revenue_growth 1 2.2470e+05 8.0405e+09 2962.2
## - rd_investment 1 2.9983e+09 1.1039e+10 3011.1
## - average_salary 1 4.1051e+10 4.9091e+10 3260.3
cat("\nFORMULA MODEL STEPWISE BOTH:\n")
##
## FORMULA MODEL STEPWISE BOTH:
print(formula(step_model))
## net_profit ~ rd_investment + average_salary
cat("\nRINGKASAN MODEL STEPWISE BOTH:\n")
##
## RINGKASAN MODEL STEPWISE BOTH:
print(summary(step_model))
##
## Call:
## lm(formula = net_profit ~ rd_investment + average_salary, data = data_model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23663 -3953 421 3234 33531
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.178e+04 1.486e+03 -7.923 3.35e-13 ***
## rd_investment 1.560e+03 1.995e+02 7.820 6.08e-13 ***
## average_salary 8.227e-01 2.843e-02 28.936 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7002 on 164 degrees of freedom
## Multiple R-squared: 0.8558, Adjusted R-squared: 0.8541
## F-statistic: 486.7 on 2 and 164 DF, p-value: < 2.2e-16
# ------------------------------------------------------------
# 12.2. STEPWISE BACKWARD
# ------------------------------------------------------------
cat("\n--- STEPWISE BACKWARD ---\n")
##
## --- STEPWISE BACKWARD ---
model_backward <- step(
model,
direction = "backward",
trace = 1)
## Start: AIC=2963.59
## net_profit ~ revenue_growth + rd_investment + operational_cost +
## average_salary + market_volatility
##
## Df Sum of Sq RSS AIC
## - revenue_growth 1 2.8949e+07 7.9452e+09 2962.2
## - operational_cost 1 5.2675e+07 7.9689e+09 2962.7
## <none> 7.9163e+09 2963.6
## - market_volatility 1 1.0118e+08 8.0174e+09 2963.7
## - rd_investment 1 2.3681e+09 1.0284e+10 3005.3
## - average_salary 1 1.8730e+10 2.6646e+10 3164.3
##
## Step: AIC=2962.2
## net_profit ~ rd_investment + operational_cost + average_salary +
## market_volatility
##
## Df Sum of Sq RSS AIC
## - operational_cost 1 2.5307e+07 7.9705e+09 2960.7
## - market_volatility 1 8.6198e+07 8.0314e+09 2962.0
## <none> 7.9452e+09 2962.2
## - rd_investment 1 2.6256e+09 1.0571e+10 3007.9
## - average_salary 1 3.9934e+10 4.7879e+10 3260.2
##
## Step: AIC=2960.73
## net_profit ~ rd_investment + average_salary + market_volatility
##
## Df Sum of Sq RSS AIC
## - market_volatility 1 7.0227e+07 8.0407e+09 2960.2
## <none> 7.9705e+09 2960.7
## - rd_investment 1 2.6138e+09 1.0584e+10 3006.1
## - average_salary 1 4.0061e+10 4.8032e+10 3258.7
##
## Step: AIC=2960.2
## net_profit ~ rd_investment + average_salary
##
## Df Sum of Sq RSS AIC
## <none> 8.0407e+09 2960.2
## - rd_investment 1 2.9983e+09 1.1039e+10 3011.1
## - average_salary 1 4.1051e+10 4.9091e+10 3260.3
cat("\nFORMULA MODEL BACKWARD:\n")
##
## FORMULA MODEL BACKWARD:
print(formula(model_backward))
## net_profit ~ rd_investment + average_salary
cat("\nRINGKASAN MODEL BACKWARD:\n")
##
## RINGKASAN MODEL BACKWARD:
print(summary(model_backward))
##
## Call:
## lm(formula = net_profit ~ rd_investment + average_salary, data = data_model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23663 -3953 421 3234 33531
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.178e+04 1.486e+03 -7.923 3.35e-13 ***
## rd_investment 1.560e+03 1.995e+02 7.820 6.08e-13 ***
## average_salary 8.227e-01 2.843e-02 28.936 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7002 on 164 degrees of freedom
## Multiple R-squared: 0.8558, Adjusted R-squared: 0.8541
## F-statistic: 486.7 on 2 and 164 DF, p-value: < 2.2e-16
# ------------------------------------------------------------
# 12.3. STEPWISE FORWARD
# ------------------------------------------------------------
cat("\n--- STEPWISE FORWARD ---\n")
##
## --- STEPWISE FORWARD ---
model_null <- lm(
net_profit ~ 1,
data = data_model)
model_forward <- step(
model_null,
scope = list(
lower = formula(model_null),
upper = formula(model)),
direction = "forward",
trace = 1)
## Start: AIC=3279.62
## net_profit ~ 1
##
## Df Sum of Sq RSS AIC
## + average_salary 1 4.4727e+10 1.1039e+10 3011.1
## + revenue_growth 1 9.7775e+09 4.5989e+10 3249.4
## + rd_investment 1 6.6748e+09 4.9091e+10 3260.3
## + market_volatility 1 2.7393e+09 5.3027e+10 3273.2
## + operational_cost 1 7.4391e+08 5.5022e+10 3279.4
## <none> 5.5766e+10 3279.6
##
## Step: AIC=3011.12
## net_profit ~ average_salary
##
## Df Sum of Sq RSS AIC
## + rd_investment 1 2998292133 8.0407e+09 2960.2
## + market_volatility 1 454684090 1.0584e+10 3006.1
## + revenue_growth 1 147918877 1.0891e+10 3010.9
## <none> 1.1039e+10 3011.1
## + operational_cost 1 1953841 1.1037e+10 3013.1
##
## Step: AIC=2960.2
## net_profit ~ average_salary + rd_investment
##
## Df Sum of Sq RSS AIC
## <none> 8040742983 2960.2
## + market_volatility 1 70227284 7970515699 2960.7
## + operational_cost 1 9335996 8031406987 2962.0
## + revenue_growth 1 224703 8040518280 2962.2
cat("\nFORMULA MODEL FORWARD:\n")
##
## FORMULA MODEL FORWARD:
print(formula(model_forward))
## net_profit ~ average_salary + rd_investment
cat("\nRINGKASAN MODEL FORWARD:\n")
##
## RINGKASAN MODEL FORWARD:
print(summary(model_forward))
##
## Call:
## lm(formula = net_profit ~ average_salary + rd_investment, data = data_model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23663 -3953 421 3234 33531
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.178e+04 1.486e+03 -7.923 3.35e-13 ***
## average_salary 8.227e-01 2.843e-02 28.936 < 2e-16 ***
## rd_investment 1.560e+03 1.995e+02 7.820 6.08e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7002 on 164 degrees of freedom
## Multiple R-squared: 0.8558, Adjusted R-squared: 0.8541
## F-statistic: 486.7 on 2 and 164 DF, p-value: < 2.2e-16
# ============================================================
# PERBANDINGAN KRITERIA MODEL
# ============================================================
hitung_kriteria <- function(model_obj) {
n_model <- nobs(model_obj)
k <- length(coef(model_obj))
aic_value <- AIC(model_obj)
bic_value <- BIC(model_obj)
aicc_value <- aic_value +
(2 * k * (k + 1)) /
(n_model - k - 1)
hasil <- data.frame(
AIC = aic_value,
BIC = bic_value,
AICc = aicc_value
)
return(round(hasil, 4))
}
cat("\nMODEL AWAL:\n")
##
## MODEL AWAL:
print(hitung_kriteria(model))
## AIC BIC AICc
## 1 3439.515 3461.341 3440.04
cat("\nMODEL STEPWISE BOTH:\n")
##
## MODEL STEPWISE BOTH:
print(hitung_kriteria(step_model))
## AIC BIC AICc
## 1 3436.121 3448.593 3436.268
cat("\nMODEL STEPWISE BACKWARD:\n")
##
## MODEL STEPWISE BACKWARD:
print(hitung_kriteria(model_backward))
## AIC BIC AICc
## 1 3436.121 3448.593 3436.268
cat("\nMODEL STEPWISE FORWARD:\n")
##
## MODEL STEPWISE FORWARD:
print(hitung_kriteria(model_forward))
## AIC BIC AICc
## 1 3436.121 3448.593 3436.268
# ============================================================
# 13. RINGKASAN HASIL ANALISIS
# ============================================================
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("12. RINGKASAN HASIL ANALISIS\n")
## 12. RINGKASAN HASIL ANALISIS
cat(strrep("=", 60), "\n")
## ============================================================
# ------------------------------------------------------------
# MODEL AWAL
# ------------------------------------------------------------
cat("\n--- MODEL AWAL ---\n")
##
## --- MODEL AWAL ---
print(formula(model))
## net_profit ~ revenue_growth + rd_investment + operational_cost +
## average_salary + market_volatility
cat("\nPERSAMAAN REGRESI:\n")
##
## PERSAMAAN REGRESI:
koef <- coef(model)
cat("net_profit =",
round(koef[1], 4))
## net_profit = -10168.6
for (i in 2:length(koef)) {
tanda <- ifelse(koef[i] >= 0, " + ", " - ")
cat(tanda, abs(round(koef[i], 4)), "*", names(koef)[i])
}
## + 33.3675 * revenue_growth + 1579.282 * rd_investment - 43.9892 * operational_cost + 0.7981 * average_salary - 80.0348 * market_volatility
cat("\n")
# ------------------------------------------------------------
# SIGNIFIKANSI
# ------------------------------------------------------------
cat("\n--- UJI SIGNIFIKANSI ---\n")
##
## --- UJI SIGNIFIKANSI ---
cat("Uji F p-value:", format.pval(f_pvalue, digits = 5), "\n")
## Uji F p-value: < 2.22e-16
cat("\nUji t masing-masing variabel:\n")
##
## Uji t masing-masing variabel:
for (i in 2:nrow(coef_df)) {
cat(coef_df$Variabel[i], ": p-value =", format.pval(coef_df$p_value[i], digits = 5), "\n")
}
## revenue_growth : p-value = 0.444
## rd_investment : p-value = < 2.22e-16
## operational_cost : p-value = 0.3022
## average_salary : p-value = < 2.22e-16
## market_volatility : p-value = 0.1534
# ------------------------------------------------------------
# ASUMSI
# ------------------------------------------------------------
cat("\n--- UJI ASUMSI ---\n")
##
## --- UJI ASUMSI ---
cat("Shapiro-Wilk p-value:",
format.pval(shapiro_test$p.value, digits = 5),"\n")
## Shapiro-Wilk p-value: 2.5201e-09
cat("Anderson-Darling p-value:",
format.pval(ad_test$p.value, digits = 5),"\n")
## Anderson-Darling p-value: 1.4629e-07
cat("Breusch-Pagan p-value:",
format.pval(bp_test$p.value, digits = 5),"\n")
## Breusch-Pagan p-value: 3.7327e-14
cat("Durbin-Watson p-value:",
format.pval(dw_test$p.value, digits = 5),"\n")
## Durbin-Watson p-value: 0.21739
cat("VIF maksimum:",
round(max(vif_values), 4),"\n")
## VIF maksimum: 4.7974
cat("Ramsey RESET p-value:",
format.pval(reset_test$p.value, digits = 5),"\n")
## Ramsey RESET p-value: 2.2545e-13
# ------------------------------------------------------------
# KOEFISIEN DETERMINASI
# ------------------------------------------------------------
cat("\n--- KOEFISIEN DETERMINASI ---\n")
##
## --- KOEFISIEN DETERMINASI ---
cat("R-squared:", round(r_squared, 4), "\n")
## R-squared: 0.858
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.8536
cat("\n")
cat(strrep("=", 60), "\n")
## ============================================================
cat("ANALISIS SELESAI\n")
## ANALISIS SELESAI
cat(strrep("=", 60), "\n")
## ============================================================