# ==== 0) Pustaka ====
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.2.3
## Warning: package 'ggplot2' was built under R version 4.2.3
## Warning: package 'tidyr' was built under R version 4.2.3
## Warning: package 'purrr' was built under R version 4.2.3
## Warning: package 'dplyr' was built under R version 4.2.3
## Warning: package 'stringr' was built under R version 4.2.3
## Warning: package 'lubridate' was built under R version 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.0
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(janitor)
##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(survival)
library(broom)
## Warning: package 'broom' was built under R version 4.2.3
library(lubridate)
setwd("D:/KULIAH/Semester 7/Data Hauler")
df <- read.csv("dengan dummy.csv", sep=";", dec=",", stringsAsFactors = FALSE, na.strings = c("","NA"))
# ==== 1) Bersihkan nama kolom & konversi tipe ====
df <- df %>% janitor::clean_names()
num_cols <- c("d1","d2","d3","d4","d5","p1","p2","p3",
"start_time","stop_time","gap","status",
"hm_break_down","event_order")
df <- df %>%
mutate(across(all_of(num_cols), ~ suppressWarnings(as.numeric(.))))
# ==== Pastikan urutan event valid ====
df <- df %>%
filter(!is.na(stop_time), !is.na(start_time), !is.na(gap),
stop_time >= start_time, gap >= 0) %>%
group_by(equipment = if ("equipment" %in% names(.)) equipment else "UNIT") %>%
arrange(start_time, stop_time, .by_group = TRUE) %>%
mutate(event_order = if_else(is.na(event_order), row_number(), event_order)) %>%
ungroup()
# ==== 2) MAPPING D1..D5 & P1..P3 ====
comp_map_detect <- df %>%
filter(toupper(object_name) != "ENGINE") %>%
rowwise() %>%
mutate(which_d = c("d1","d2","d3","d4","d5")[which.max(c_across(d1:d5))]) %>%
ungroup() %>%
count(object_name, which_d, sort = TRUE) %>%
group_by(object_name) %>%
slice_max(n, n = 1, with_ties = FALSE) %>%
ungroup()
comp_map <- setNames(comp_map_detect$object_name, comp_map_detect$which_d)
comp_baseline <- "ENGINE"
plant_map_detect <- df %>%
filter(toupper(maintenance_plant) != "ADT") %>%
rowwise() %>%
mutate(which_p = c("p1","p2","p3")[which.max(c_across(p1:p3))]) %>%
ungroup() %>%
count(maintenance_plant, which_p, sort = TRUE) %>%
group_by(maintenance_plant) %>%
slice_max(n, n = 1, with_ties = FALSE) %>%
ungroup()
plant_map <- setNames(plant_map_detect$maintenance_plant, plant_map_detect$which_p)
plant_baseline <- "ADT"
message("Mapping komponen terdeteksi:")
## Mapping komponen terdeteksi:
print(comp_map)
## d1 d2 d3 d4
## "DIFFERENTIAL" "FINAL DRIVE" "FRAME & CHASSIS" "SUSPENSION"
## d5
## "TRANSMISSION"
message("Mapping plant terdeteksi:")
## Mapping plant terdeteksi:
print(plant_map)
## p1 p2 p3
## "BIN" "IPR" "SDJ"
# ==== 3) VALIDASI BASELINE ====
engine_ok <- df %>% filter(toupper(object_name)==comp_baseline) %>%
summarise(ok = all(d1==0 & d2==0 & d3==0 & d4==0 & d5==0)) %>% pull(ok)
adt_ok <- df %>% filter(toupper(maintenance_plant)==plant_baseline) %>%
summarise(ok = all(p1==0 & p2==0 & p3==0)) %>% pull(ok)
message("Baseline ENGINE all-zero? ", engine_ok)
## Baseline ENGINE all-zero? TRUE
message("Baseline ADT all-zero? ", adt_ok)
## Baseline ADT all-zero? TRUE
# ==== 4) PARETO ====
pareto_df <- df %>%
filter(status == 1) %>%
count(object_name, name = "freq") %>%
arrange(desc(freq)) %>%
mutate(prop = freq/sum(freq),
cum_prop = cumsum(prop),
object_name = forcats::fct_reorder(object_name, freq))
gg_pareto <- ggplot(pareto_df, aes(x = object_name, y = freq)) +
geom_col(width = 0.7) +
geom_point(aes(y = cum_prop * max(freq)), size = 2) +
geom_line(aes(y = cum_prop * max(freq), group = 1)) +
geom_hline(yintercept = 0.8 * max(pareto_df$freq), linetype = "dashed") +
scale_y_continuous(name = "Frekuensi Kerusakan",
sec.axis = sec_axis(~ . / max(pareto_df$freq),
labels = scales::percent, name = "Kumulatif")) +
labs(x = "Komponen (Object Name)",
title = "Pareto Kerusakan per Komponen") +
theme_minimal(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(gg_pareto)
# ==== 4B) ANALISIS DESKRIPTIF DATA KERUSAKAN ====
# 1️⃣ Frekuensi kerusakan per unit dan per komponen
freq_unit <- df %>%
filter(status == 1) %>%
count(equipment, name = "freq") %>%
arrange(desc(freq))
freq_comp <- df %>%
filter(status == 1) %>%
count(object_name, name = "freq") %>%
arrange(desc(freq))
message("Top 10 Unit dengan frekuensi kerusakan tertinggi:")
## Top 10 Unit dengan frekuensi kerusakan tertinggi:
print(head(freq_unit, 10))
## # A tibble: 10 × 2
## equipment freq
## <chr> <int>
## 1 HDKM78085 51
## 2 HDKM78193 51
## 3 HDKM78196 51
## 4 HDKM78205 51
## 5 HDKM78213 51
## 6 HDKM78215 51
## 7 HDKM78216 51
## 8 HDKM78220 51
## 9 HDKM78222 51
## 10 HDKM78229 51
message("Frekuensi kerusakan per komponen:")
## Frekuensi kerusakan per komponen:
print(freq_comp)
## # A tibble: 6 × 2
## object_name freq
## <chr> <int>
## 1 ENGINE 3259
## 2 SUSPENSION 1307
## 3 FRAME & CHASSIS 1284
## 4 TRANSMISSION 489
## 5 FINAL DRIVE 62
## 6 DIFFERENTIAL 39
freq_comp <- df %>%
filter(status == 1) %>%
count(object_name, name = "freq") %>%
arrange(desc(freq)) %>%
mutate(persentase = round(freq / sum(freq) * 100,2))
gg_freq_comp <- ggplot(freq_comp,
aes(x = fct_reorder(object_name, freq), y = freq)) +
geom_col(fill = "steelblue", width = 0.7) +
geom_text(
aes(label = paste0(freq, "\n", round(persentase, 1), "%")),
hjust = -0.15,
size = 4,
fontface = "bold",
lineheight = 0.9
) +
coord_flip() +
scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
labs(
x = "Komponen",
y = "Frekuensi Kerusakan",
title = "Frekuensi dan Persentase Kerusakan per Komponen"
) +
theme_minimal(base_size = 12) +
theme(
axis.title = element_text(face = "bold", size = 14),
axis.text = element_text(size = 12),
plot.title = element_text(face = "bold", size = 14, hjust = 0.5)
)
print(gg_freq_comp)
freq_site <- df %>%
filter(status == 1) %>%
count(maintenance_plant, name = "freq") %>%
arrange(desc(freq)) %>%
mutate(persentase = round(freq / sum(freq) * 100,2))
gg_freq_site <- ggplot(freq_site,
aes(x = fct_reorder(maintenance_plant, freq), y = freq)) +
geom_col(fill = "darkgreen", width = 0.7) +
geom_text(
aes(label = paste0(freq, "\n", persentase, "%")),
hjust = -0.15,
size = 4,
fontface = "bold",
lineheight = 0.9
) +
coord_flip() +
scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
labs(
x = "Site",
y = "Frekuensi Kerusakan",
title = "Frekuensi dan Persentase Kerusakan per Site"
) +
theme_minimal(base_size = 12) +
theme(
axis.title = element_text(face = "bold", size = 14), # <<< INI
axis.text = element_text(size = 12),
plot.title = element_text(face = "bold", size = 14, hjust = 0.5)
)
print(gg_freq_site)
# 2️⃣ Analisis gap time antar kerusakan
# Jika belum ada kolom gap_gt, hitung kembali untuk eksplorasi
df_gap <- df %>%
group_by(equipment) %>%
arrange(stop_time, .by_group = TRUE) %>%
mutate(
prev_stop = lag(stop_time),
gap_time = if_else(!is.na(prev_stop),
stop_time - prev_stop,
pmax(stop_time - start_time, 0))
) %>%
ungroup() %>%
filter(!is.na(gap_time), gap_time > 0)
# Statistik deskriptif gap time
gap_summary <- df_gap %>%
summarise(
n = n(),
mean_gap = mean(gap_time, na.rm = TRUE),
median_gap = median(gap_time, na.rm = TRUE),
var_gap = var(gap_time, na.rm = TRUE),
sd_gap = sd(gap_time, na.rm = TRUE),
min_gap = min(gap_time, na.rm = TRUE),
max_gap = max(gap_time, na.rm = TRUE)
)
message("Ringkasan statistik gap time:")
## Ringkasan statistik gap time:
print(gap_summary)
## # A tibble: 1 × 7
## n mean_gap median_gap var_gap sd_gap min_gap max_gap
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 6433 123. 58.5 37011. 192. 0.0100 3352
# 3️⃣ Visualisasi distribusi gap time
gg_gap_hist <- ggplot(df_gap, aes(x = gap_time)) +
geom_histogram(binwidth = diff(range(df_gap$gap_time, na.rm = TRUE))/30,
fill = "orange", color = "black", alpha = 0.7) +
labs(x = "Gap Time (jam atau HM)", y = "Frekuensi",
title = "Distribusi Interval Waktu Antar Kerusakan (Gap Time)") +
theme_minimal(base_size = 12)
print(gg_gap_hist)
# ==== 5) UJI LAPLACE ====
laplace_test <- function(times_vec){
times_vec <- sort(times_vec)
n <- length(times_vec)
if(n < 2) return(tibble(Z = NA_real_, p_value = NA_real_, n = n, T = NA_real_))
Ttot <- max(times_vec)
S <- sum((1:n) * times_vec)
Z <- (S/n - 0.5*Ttot) / sqrt(Ttot^3/(12*n))
tibble(Z = Z, p_value = 2*(1 - pnorm(abs(Z))), n = n, T = Ttot)
}
laplace_overall <- laplace_test(df %>% filter(status == 1) %>% pull(stop_time))
print(laplace_overall)
## # A tibble: 1 × 4
## Z p_value n T
## <dbl> <dbl> <int> <dbl>
## 1 3444. 0 6440 23957
# === 7) PWP–GT (baseline) ===
# Hitung ulang gap_gt sesuai definisi PWP-GT
df_gt <- df %>%
group_by(equipment) %>%
arrange(stop_time, .by_group = TRUE) %>%
mutate(
prev_stop = lag(stop_time),
gap_gt = if_else(
!is.na(prev_stop),
stop_time - prev_stop, # gap antar kerusakan
pmax(stop_time - start_time, 0) # event pertama
),
gap_gt = pmax(gap_gt, 0) # tidak boleh negatif
) %>%
ungroup() %>%
select(-prev_stop) %>%
mutate(
start_gt = 0,
stop_gt = gap_gt
)
# Hanya ambil event dengan durasi valid
df_gt <- df_gt %>% filter(stop_gt > 0 | status == 0)
# === Model PWP–GT ===
form_gt <- as.formula(
Surv(start_gt, stop_gt, status) ~
d1 + d2 + d3 + d4 + d5 +
p1 + p2 + p3 +
strata(event_order) +
cluster(equipment)
)
fit_pwp_gt <- coxph(form_gt, data = df_gt, ties = "efron")
summary(fit_pwp_gt)
## Call:
## coxph(formula = Surv(start_gt, stop_gt, status) ~ d1 + d2 + d3 +
## d4 + d5 + p1 + p2 + p3 + strata(event_order), data = df_gt,
## ties = "efron", cluster = equipment)
##
## n= 6433, number of events= 6433
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## d1 -0.37603 0.68658 0.16691 0.16437 -2.288 0.0222 *
## d2 -0.00532 0.99469 0.13341 0.21462 -0.025 0.9802
## d3 -0.22467 0.79878 0.03372 0.03696 -6.078 1.22e-09 ***
## d4 -0.16488 0.84800 0.03378 0.04006 -4.116 3.85e-05 ***
## d5 -0.08907 0.91478 0.04988 0.06181 -1.441 0.1496
## p1 -0.02952 0.97092 0.04531 0.04416 -0.668 0.5039
## p2 -0.05548 0.94603 0.02759 0.04464 -1.243 0.2139
## p3 0.08693 1.09082 0.07214 0.05829 1.491 0.1359
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## d1 0.6866 1.4565 0.4975 0.9476
## d2 0.9947 1.0053 0.6531 1.5149
## d3 0.7988 1.2519 0.7430 0.8588
## d4 0.8480 1.1792 0.7840 0.9173
## d5 0.9148 1.0932 0.8104 1.0326
## p1 0.9709 1.0300 0.8904 1.0587
## p2 0.9460 1.0570 0.8668 1.0325
## p3 1.0908 0.9167 0.9731 1.2228
##
## Concordance= 0.533 (se = 0.005 )
## Likelihood ratio test= 66.47 on 8 df, p=2e-11
## Wald test = 51.42 on 8 df, p=2e-08
## Score (logrank) test = 65.96 on 8 df, p=3e-11, Robust = 37.54 p=9e-06
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).
# Uji asumsi PH
ph_gt <- cox.zph(fit_pwp_gt, transform = "km")
print(ph_gt)
## chisq df p
## d1 0.3435 1 0.5578
## d2 4.9519 1 0.0261
## d3 10.5054 1 0.0012
## d4 0.3337 1 0.5635
## d5 1.6923 1 0.1933
## p1 0.0298 1 0.8629
## p2 4.0338 1 0.0446
## p3 0.6319 1 0.4266
## GLOBAL 22.0286 8 0.0049
# ===============================================================
# 8) PWP–GT + TVC BERDASARKAN QUARTER KALENDER (Q3_2024 → Q2_2025)
# ===============================================================
df <- df %>%
mutate(
breakdown_date = suppressWarnings(parse_date_time(
breakdown_date,
orders = c("d/m/Y","d-m-Y","Y-m-d","m/d/Y","m-d-Y")
)),
month = lubridate::month(breakdown_date),
year = lubridate::year(breakdown_date),
# Mapping quarter aktual:
cal_interval4 = case_when(
year == 2024 & month %in% 7:9 ~ "Q3_2024", # baseline
year == 2024 & month %in% 10:12 ~ "Q4_2024",
year == 2025 & month %in% 1:3 ~ "Q1_2025",
year == 2025 & month %in% 4:6 ~ "Q2_2025",
TRUE ~ NA_character_
),
cal_interval4 = factor(cal_interval4,
levels = c("Q3_2024","Q4_2024","Q1_2025","Q2_2025"))
)
# ==== Setup ulang untuk PWP–GT + TVC ====
df_tv <- df %>%
filter(!is.na(cal_interval4)) %>%
group_by(equipment) %>%
arrange(start_time, stop_time, .by_group = TRUE) %>%
mutate(event_order = row_number()) %>%
ungroup()
df_tv <- df_tv %>%
group_by(equipment) %>%
arrange(stop_time, .by_group = TRUE) %>%
mutate(
prev_stop = lag(stop_time),
gap_gt = if_else(
!is.na(prev_stop),
stop_time - prev_stop, # gap antar kejadian
pmax(stop_time - start_time, 0) # event pertama
),
gap_gt = pmax(gap_gt, 0)
) %>%
ungroup() %>%
select(-prev_stop)
df_tv <- df_tv %>%
mutate(
start_gt = 0,
stop_gt = gap_gt
) %>%
filter(stop_gt > 0 | status == 0)
# ==== FINAL MODEL PWP–GT + TVC (Quarter) ====
form_gt_tv <- as.formula(
Surv(start_gt, stop_gt, status) ~
(d1 + d2 + d3 + d4 + d5 + p1 + p2 + p3) * cal_interval4 +
strata(event_order) +
cluster(equipment)
)
fit_pwp_gt_tv <- coxph(form_gt_tv, data = df_tv, ties = "efron")
summary(fit_pwp_gt_tv)
## Call:
## coxph(formula = Surv(start_gt, stop_gt, status) ~ d1 + d2 + d3 +
## d4 + d5 + p1 + p2 + p3 + cal_interval4 + strata(event_order) +
## d1:cal_interval4 + d2:cal_interval4 + d3:cal_interval4 +
## d4:cal_interval4 + d5:cal_interval4 + p1:cal_interval4 +
## p2:cal_interval4 + p3:cal_interval4, data = df_tv, ties = "efron",
## cluster = equipment)
##
## n= 6412, number of events= 6412
##
## coef exp(coef) se(coef) robust se z
## d1 -0.2141420 0.8072337 0.3392866 0.1808961 -1.184
## d2 0.1679933 1.1829287 0.2841880 0.3143783 0.534
## d3 -0.1588208 0.8531493 0.0698437 0.0682713 -2.326
## d4 -0.1679335 0.8454101 0.0732890 0.0766060 -2.192
## d5 -0.0611987 0.9406363 0.1097452 0.0979682 -0.625
## p1 0.1362330 1.1459489 0.0788674 0.0766248 1.778
## p2 -0.0001564 0.9998436 0.0610926 0.0513441 -0.003
## p3 0.1835470 1.2014714 0.1684674 0.0973863 1.885
## cal_interval4Q4_2024 -0.2853474 0.7517530 0.0655687 0.0792357 -3.601
## cal_interval4Q1_2025 -0.5040617 0.6040721 0.0685484 0.0894079 -5.638
## cal_interval4Q2_2025 -0.6078132 0.5445404 0.0705488 0.0907712 -6.696
## d1:cal_interval4Q4_2024 -0.0163238 0.9838087 0.4510959 0.3056123 -0.053
## d1:cal_interval4Q1_2025 -0.0837809 0.9196327 0.4373875 0.3198410 -0.262
## d1:cal_interval4Q2_2025 -1.2584322 0.2840991 0.6799865 0.4657233 -2.702
## d2:cal_interval4Q4_2024 -0.6099597 0.5433727 0.4231985 0.5249778 -1.162
## d2:cal_interval4Q1_2025 0.0836979 1.0873003 0.3583106 0.5433764 0.154
## d2:cal_interval4Q2_2025 -0.2294096 0.7950028 0.4015424 0.3846385 -0.596
## d3:cal_interval4Q4_2024 -0.0508935 0.9503799 0.0963145 0.0917627 -0.555
## d3:cal_interval4Q1_2025 -0.0724035 0.9301555 0.0988453 0.1088910 -0.665
## d3:cal_interval4Q2_2025 -0.0927890 0.9113858 0.0960576 0.0896248 -1.035
## d4:cal_interval4Q4_2024 -0.0200480 0.9801517 0.1019346 0.1064305 -0.188
## d4:cal_interval4Q1_2025 0.0693768 1.0718400 0.0970716 0.1082046 0.641
## d4:cal_interval4Q2_2025 0.0303311 1.0307958 0.0986969 0.0993338 0.305
## d5:cal_interval4Q4_2024 -0.0696815 0.9326908 0.1493201 0.1532793 -0.455
## d5:cal_interval4Q1_2025 0.0446781 1.0456912 0.1491046 0.1395457 0.320
## d5:cal_interval4Q2_2025 -0.0266760 0.9736766 0.1438545 0.1429170 -0.187
## p1:cal_interval4Q4_2024 -0.2827991 0.7536712 0.1166694 0.1398628 -2.022
## p1:cal_interval4Q1_2025 -0.3404903 0.7114214 0.1346882 0.1968493 -1.730
## p1:cal_interval4Q2_2025 -0.2986517 0.7418177 0.1299465 0.1572072 -1.900
## p2:cal_interval4Q4_2024 -0.1219246 0.8852151 0.0836602 0.0888352 -1.372
## p2:cal_interval4Q1_2025 -0.0833683 0.9200122 0.0810994 0.0892318 -0.934
## p2:cal_interval4Q2_2025 0.0391287 1.0399043 0.0811575 0.0895767 0.437
## p3:cal_interval4Q4_2024 -0.0184044 0.9817639 0.2176371 0.2182787 -0.084
## p3:cal_interval4Q1_2025 0.0066461 1.0066683 0.2078511 0.1491779 0.045
## p3:cal_interval4Q2_2025 -0.5671852 0.5671195 0.2394397 0.2020950 -2.807
## Pr(>|z|)
## d1 0.236498
## d2 0.593088
## d3 0.020002 *
## d4 0.028367 *
## d5 0.532182
## p1 0.075416 .
## p2 0.997569
## p3 0.059466 .
## cal_interval4Q4_2024 0.000317 ***
## cal_interval4Q1_2025 1.72e-08 ***
## cal_interval4Q2_2025 2.14e-11 ***
## d1:cal_interval4Q4_2024 0.957403
## d1:cal_interval4Q1_2025 0.793363
## d1:cal_interval4Q2_2025 0.006890 **
## d2:cal_interval4Q4_2024 0.245285
## d2:cal_interval4Q1_2025 0.877584
## d2:cal_interval4Q2_2025 0.550889
## d3:cal_interval4Q4_2024 0.579154
## d3:cal_interval4Q1_2025 0.506104
## d3:cal_interval4Q2_2025 0.300527
## d4:cal_interval4Q4_2024 0.850589
## d4:cal_interval4Q1_2025 0.521417
## d4:cal_interval4Q2_2025 0.760103
## d5:cal_interval4Q4_2024 0.649394
## d5:cal_interval4Q1_2025 0.748841
## d5:cal_interval4Q2_2025 0.851932
## p1:cal_interval4Q4_2024 0.043179 *
## p1:cal_interval4Q1_2025 0.083684 .
## p1:cal_interval4Q2_2025 0.057468 .
## p2:cal_interval4Q4_2024 0.169914
## p2:cal_interval4Q1_2025 0.350155
## p2:cal_interval4Q2_2025 0.662244
## p3:cal_interval4Q4_2024 0.932805
## p3:cal_interval4Q1_2025 0.964465
## p3:cal_interval4Q2_2025 0.005008 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## d1 0.8072 1.2388 0.5663 1.1507
## d2 1.1829 0.8454 0.6388 2.1906
## d3 0.8531 1.1721 0.7463 0.9753
## d4 0.8454 1.1829 0.7275 0.9824
## d5 0.9406 1.0631 0.7763 1.1398
## p1 1.1459 0.8726 0.9861 1.3316
## p2 0.9998 1.0002 0.9041 1.1057
## p3 1.2015 0.8323 0.9927 1.4541
## cal_interval4Q4_2024 0.7518 1.3302 0.6436 0.8781
## cal_interval4Q1_2025 0.6041 1.6554 0.5070 0.7198
## cal_interval4Q2_2025 0.5445 1.8364 0.4558 0.6506
## d1:cal_interval4Q4_2024 0.9838 1.0165 0.5405 1.7908
## d1:cal_interval4Q1_2025 0.9196 1.0874 0.4913 1.7213
## d1:cal_interval4Q2_2025 0.2841 3.5199 0.1140 0.7078
## d2:cal_interval4Q4_2024 0.5434 1.8404 0.1942 1.5204
## d2:cal_interval4Q1_2025 1.0873 0.9197 0.3748 3.1541
## d2:cal_interval4Q2_2025 0.7950 1.2579 0.3741 1.6896
## d3:cal_interval4Q4_2024 0.9504 1.0522 0.7939 1.1376
## d3:cal_interval4Q1_2025 0.9302 1.0751 0.7514 1.1514
## d3:cal_interval4Q2_2025 0.9114 1.0972 0.7646 1.0864
## d4:cal_interval4Q4_2024 0.9802 1.0203 0.7956 1.2075
## d4:cal_interval4Q1_2025 1.0718 0.9330 0.8670 1.3251
## d4:cal_interval4Q2_2025 1.0308 0.9701 0.8484 1.2523
## d5:cal_interval4Q4_2024 0.9327 1.0722 0.6907 1.2595
## d5:cal_interval4Q1_2025 1.0457 0.9563 0.7955 1.3746
## d5:cal_interval4Q2_2025 0.9737 1.0270 0.7358 1.2884
## p1:cal_interval4Q4_2024 0.7537 1.3268 0.5730 0.9914
## p1:cal_interval4Q1_2025 0.7114 1.4056 0.4837 1.0464
## p1:cal_interval4Q2_2025 0.7418 1.3480 0.5451 1.0095
## p2:cal_interval4Q4_2024 0.8852 1.1297 0.7438 1.0536
## p2:cal_interval4Q1_2025 0.9200 1.0869 0.7724 1.0958
## p2:cal_interval4Q2_2025 1.0399 0.9616 0.8725 1.2395
## p3:cal_interval4Q4_2024 0.9818 1.0186 0.6400 1.5059
## p3:cal_interval4Q1_2025 1.0067 0.9934 0.7515 1.3485
## p3:cal_interval4Q2_2025 0.5671 1.7633 0.3816 0.8427
##
## Concordance= 0.558 (se = 0.006 )
## Likelihood ratio test= 298.9 on 35 df, p=<2e-16
## Wald test = 263.1 on 35 df, p=<2e-16
## Score (logrank) test = 299.3 on 35 df, p=<2e-16, Robust = 92.6 p=4e-07
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).
# ===============================================================
# 9 UJI SIGNIFIKANSI PARSIAL (WALD TEST + HR)
# ===============================================================
coef_final <- tidy(fit_pwp_gt_tv, exponentiate = TRUE, conf.int = TRUE)
print(coef_final)
## # A tibble: 35 × 8
## term estimate std.error robust.se statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 d1 0.807 0.339 0.181 -1.18 2.36e-1 0.566 1.15
## 2 d2 1.18 0.284 0.314 0.534 5.93e-1 0.639 2.19
## 3 d3 0.853 0.0698 0.0683 -2.33 2.00e-2 0.746 0.975
## 4 d4 0.845 0.0733 0.0766 -2.19 2.84e-2 0.728 0.982
## 5 d5 0.941 0.110 0.0980 -0.625 5.32e-1 0.776 1.14
## 6 p1 1.15 0.0789 0.0766 1.78 7.54e-2 0.986 1.33
## 7 p2 1.00 0.0611 0.0513 -0.00305 9.98e-1 0.904 1.11
## 8 p3 1.20 0.168 0.0974 1.88 5.95e-2 0.993 1.45
## 9 cal_interv… 0.752 0.0656 0.0792 -3.60 3.17e-4 0.644 0.878
## 10 cal_interv… 0.604 0.0685 0.0894 -5.64 1.72e-8 0.507 0.720
## # ℹ 25 more rows
# ===============================================================
# 10 BASELINE HAZARD MODEL AKHIR
# ===============================================================
base_haz <- basehaz(fit_pwp_gt_tv)
## Warning in survfit.coxph(fit, se.fit = FALSE): the model contains interactions;
## the default curve based on columm means of the X matrix is almost certainly not
## useful. Consider adding a newdata argument.
head(base_haz)
## hazard time strata
## 1 0.005328507 1.2 event_order=1
## 2 0.010681344 1.3 event_order=1
## 3 0.016058512 4.8 event_order=1
## 4 0.021447949 6.3 event_order=1
## 5 0.026851756 9.0 event_order=1
## 6 0.032295439 10.8 event_order=1
haz_rate <- diff(base_haz$hazard) / diff(base_haz$time)
head(haz_rate)
## [1] 0.053528374 0.001536333 0.003592958 0.002001410 0.003024268 0.003127702