Ringkasan Data

library(dplyr)
library(tidyr)

# Fungsi untuk membuat tabel frekuensi + persentase
make_wide_table <- function(varname, label_name) {
  mydata %>%
    count(!!sym(varname)) %>%
    mutate(
      n_pct = paste0(n, " (", round(100 * n / sum(n), 1), "%)")
    ) %>%
    select(Category = !!sym(varname), n_pct) %>%
    pivot_wider(names_from = Category, values_from = n_pct) %>%
    mutate(Variable = label_name) %>%
    select(Variable, everything())
}

# Daftar nama variabel
variabel_list <- c("STATUS","KRT", "WILAYAH", "KAWIN", "DIDIK", "JK", 
                   "SERTIF_LATIH", "SERTIF_MAGANG", "MRISEN", "UMK")

# Gunakan nama label sama dengan nama variabel (bisa disesuaikan kalau perlu)
label_list <- variabel_list

# Buat tabel ringkas
tabel_ringkas_lain <- purrr::map2_dfr(variabel_list, label_list, make_wide_table)

# Tampilkan tabel
print(tabel_ringkas_lain)

## # A tibble: 10 × 21
##    Variable     `0`   `1`   `bukan KRT` KRT   Perdesaan Perkotaan `Pernah Kawin`
##    <chr>        <chr> <chr> <chr>       <chr> <chr>     <chr>     <chr>         
##  1 STATUS       290 … 88 (… <NA>        <NA>  <NA>      <NA>      <NA>          
##  2 KRT          <NA>  <NA>  375 (99.2%) 3 (0… <NA>      <NA>      <NA>          
##  3 WILAYAH      <NA>  <NA>  <NA>        <NA>  63 (16.7… 315 (83.… <NA>          
##  4 KAWIN        <NA>  <NA>  <NA>        <NA>  <NA>      <NA>      8 (2.1%)      
##  5 DIDIK        <NA>  <NA>  <NA>        <NA>  <NA>      <NA>      <NA>          
##  6 JK           <NA>  <NA>  <NA>        <NA>  <NA>      <NA>      <NA>          
##  7 SERTIF_LATIH <NA>  <NA>  <NA>        <NA>  <NA>      <NA>      <NA>          
##  8 SERTIF_MAGA… <NA>  <NA>  <NA>        <NA>  <NA>      <NA>      <NA>          
##  9 MRISEN       <NA>  <NA>  <NA>        <NA>  <NA>      <NA>      <NA>          
## 10 UMK          <NA>  <NA>  <NA>        <NA>  <NA>      <NA>      <NA>          
## # ℹ 13 more variables: `Belum Pernah Kawin` <chr>, `SMA/Sederajat` <chr>,
## #   `D1/D2/D3` <chr>, `D4/S1` <chr>, Perempuan <chr>, `Laki-laki` <chr>,
## #   `Tidak Memiliki` <chr>, Memiliki <chr>, `Bukan Migran` <chr>, Migran <chr>,
## #   rendah <chr>, sedang <chr>, tinggi <chr>

summary(mydata$TIME)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.000   6.000   9.619  14.000  72.000

summary(mydata$UMUR)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   17.00   19.00   20.00   20.31   22.00   24.00

Kaplan Maiyer

# Library
library(survival)
library(survminer)

## Loading required package: ggpubr

## 
## Attaching package: 'survminer'

## The following object is masked from 'package:survival':
## 
##     myeloma

library(flexsurv)
library(eha)

## 
## Attaching package: 'eha'

## The following objects are masked from 'package:flexsurv':
## 
##     dgompertz, dllogis, hgompertz, Hgompertz, hllogis, Hllogis, hlnorm,
##     Hlnorm, hweibull, Hweibull, pgompertz, pllogis, qgompertz, qllogis,
##     rgompertz, rllogis

## The following object is masked from 'package:Hmisc':
## 
##     logrank

library(haven) # jika impor dari .sav

# Asumsi dataset sudah dibersihkan dan bernama `data`
# Buat objek Surv
mydata$s_obj <- Surv(time = mydata$TIME, event = mydata$STATUS)

# Library yang diperlukan
library(survival)
library(survminer)


# Buat model Kaplan-Meier
km_fit <- survfit(Surv(time = mydata$TIME, event = mydata$STATUS) ~ 1)

# Plot KM yang telah ditingkatkan
# Plot KM yang telah ditingkatkan - tanpa p-value karena model tidak memiliki grup
ggsurvplot(
  fit = km_fit,
  data = mydata,
  conf.int = TRUE,                   # Interval kepercayaan
  pval = FALSE,                      # Tidak perlu nilai p karena tidak ada grup dibandingkan
  risk.table = TRUE,                # Tabel risiko
  risk.table.title = "Jumlah individu yang berisiko per waktu",
  risk.table.y.text.col = TRUE,     
  risk.table.y.text = FALSE,
  surv.median.line = "hv",          
  xlab = "Waktu (bulan)",           
  ylab = "Probabilitas Survival",   
  ggtheme = theme_bw(base_size = 14),
  palette = c("#0073C2FF"),         
  censor.shape = 124,               
  censor.size = 3,
  legend = "bottom",                
  legend.labs = c("Keseluruhan Sample"), 
  break.time.by = 6
) +
  labs(
    title = "Kurva Kaplan-Meier Keseluruhan",
    subtitle = "Dengan interval kepercayaan dan garis median",
    caption = "Sumber: Sakernas 2023"
  )

Uji Log rank

library(survival)
library(survminer)
variabel_list <- c("STATUS","KRT", "WILAYAH", "KAWIN", "DIDIK", "JK", 
                   "SERTIF_LATIH", "SERTIF_MAGANG", "MRISEN", "UMK")
# Loop eksplorasi Kaplan-Meier dan log-rank
for (v in variabel_list) {
  cat("\n\n==============================\n")
  cat("Variable:", v, "\n")
  
  # Coba-catch agar error tidak hentikan loop
  tryCatch({
    # Survival object

    
    # Pastikan kolom sebagai faktor (untuk KM)
    mydata[[v]] <- as.factor(mydata[[v]])
    
    # Model survival
    fit <- survfit(mydata$s_obj ~ mydata[[v]])
    
    # Plot KM
    print(ggsurvplot(
      fit, 
      data = mydata,           # <- ganti dari data.frame(group_var)
      pval = TRUE,
      conf.int = TRUE,
      risk.table = TRUE,
      ggtheme = theme_minimal(),
      title = paste("Kaplan-Meier -", v),
      legend.title = v,
      legend.labs = levels(mydata[[v]])
    ))
    
    # Uji log-rank
    print(survdiff(mydata$s_obj ~ mydata[[v]], data = mydata))
    
    # Plot log-log survival
    print(ggsurvplot(
      fit, 
      data = mydata,
      fun = "cloglog", 
      conf.int = FALSE,
      ggtheme = theme_minimal(),
      title = paste("Log-Log Survival -", v),
      legend.title = v,
      legend.labs = levels(mydata[[v]])
    ))
  }, error = function(e) {
    cat("Error processing variable", v, ":", conditionMessage(e), "\n")
  })
}

## 
## 
## ==============================
## Variable: STATUS

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=0 290        0     71.1      71.1       414
## mydata[[v]]=1  88       88     16.9     299.0       414
## 
##  Chisq= 414  on 1 degrees of freedom, p= <2e-16

## 
## 
## ==============================
## Variable: KRT

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                         N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=bukan KRT 375       86   87.516    0.0263         5
## mydata[[v]]=KRT         3        2    0.484    4.7498         5
## 
##  Chisq= 5  on 1 degrees of freedom, p= 0.03

## 
## 
## ==============================
## Variable: WILAYAH

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                         N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=Perdesaan  63       11     15.7     1.395      1.77
## mydata[[v]]=Perkotaan 315       77     72.3     0.302      1.77
## 
##  Chisq= 1.8  on 1 degrees of freedom, p= 0.2

## 
## 
## ==============================
## Variable: KAWIN

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                                  N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=Pernah Kawin         8        3     1.88    0.6738     0.719
## mydata[[v]]=Belum Pernah Kawin 370       85    86.12    0.0147     0.719
## 
##  Chisq= 0.7  on 1 degrees of freedom, p= 0.4

## 
## 
## ==============================
## Variable: DIDIK

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                             N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=SMA/Sederajat 332       71    78.78     0.769      7.69
## mydata[[v]]=D1/D2/D3       15        5     2.98     1.377      1.49
## mydata[[v]]=D4/S1          31       12     6.24     5.313      5.98
## 
##  Chisq= 7.8  on 2 degrees of freedom, p= 0.02

## 
## 
## ==============================
## Variable: JK

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                         N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=Perempuan 171       40     38.5    0.0623     0.116
## mydata[[v]]=Laki-laki 207       48     49.5    0.0484     0.116
## 
##  Chisq= 0.1  on 1 degrees of freedom, p= 0.7

## 
## 
## ==============================
## Variable: SERTIF_LATIH

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                              N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=Tidak Memiliki 298       63     69.1     0.536      2.73
## mydata[[v]]=Memiliki        80       25     18.9     1.957      2.73
## 
##  Chisq= 2.7  on 1 degrees of freedom, p= 0.1

## 
## 
## ==============================
## Variable: SERTIF_MAGANG

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                              N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=Tidak Memiliki 345       81    81.84   0.00858     0.152
## mydata[[v]]=Memiliki        33        7     6.16   0.11395     0.152
## 
##  Chisq= 0.2  on 1 degrees of freedom, p= 0.7

## 
## 
## ==============================
## Variable: MRISEN

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                            N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=Bukan Migran 374       85   87.201    0.0555       6.4
## mydata[[v]]=Migran         4        3    0.799    6.0615       6.4
## 
##  Chisq= 6.4  on 1 degrees of freedom, p= 0.01

## 
## 
## ==============================
## Variable: UMK

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                      N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=rendah  44       11     10.8   0.00487   0.00579
## mydata[[v]]=sedang 117       25     23.9   0.04682   0.06759
## mydata[[v]]=tinggi 217       52     53.3   0.03112   0.08319
## 
##  Chisq= 0.1  on 2 degrees of freedom, p= 1

#library(survival)
library(survminer)

# Tambahkan UMUR dan JUMUR ke list variabel
variabel_list <- c("UMUR", "JUMUR", "KRT", "DIDIK", "MRISEN")

# Ubah jadi faktor sebelum loop
mydata$UMUR <- as.factor(mydata$UMUR)

# Loop eksplorasi Kaplan-Meier dan log-rank
for (v in variabel_list) {
  cat("\n\n==============================\n")
  cat("Variable:", v, "\n")
  
  tryCatch({
    # Pastikan faktor
    mydata[[v]] <- as.factor(mydata[[v]])
    
    # Survival fit
    fit <- survfit(mydata$s_obj ~ mydata[[v]])
    
    # Plot KM
    print(ggsurvplot(
      fit, 
      data = mydata,
      pval = TRUE,
      conf.int = TRUE,
      risk.table = TRUE,
      ggtheme = theme_minimal(),
      title = paste("Kaplan-Meier -", v),
      legend.title = v,
      legend.labs = levels(mydata[[v]])
    ))
    
    # Log-rank test
    print(survdiff(mydata$s_obj ~ mydata[[v]], data = mydata))
    
    # Plot log-log survival
    print(ggsurvplot(
      fit, 
      data = mydata,
      fun = "cloglog", 
      conf.int = FALSE,
      ggtheme = theme_minimal(),
      title = paste("Log-Log Survival -", v),
      legend.title = v,
      legend.labs = levels(mydata[[v]])
    ))
  }, error = function(e) {
    cat("Error processing variable", v, ":", conditionMessage(e), "\n")
  })
}

## 
## 
## ==============================
## Variable: UMUR

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=17 15        2     1.35   0.31839   0.34902
## mydata[[v]]=18 72       13     9.93   0.95001   1.15221
## mydata[[v]]=19 73       21    17.43   0.73210   0.95748
## mydata[[v]]=20 57       16    15.25   0.03679   0.04645
## mydata[[v]]=21 46       12    12.55   0.02375   0.02888
## mydata[[v]]=22 42        6    11.99   2.99508   3.62625
## mydata[[v]]=23 41        9    10.70   0.26928   0.31960
## mydata[[v]]=24 32        9     8.81   0.00409   0.00527
## 
##  Chisq= 5.7  on 7 degrees of freedom, p= 0.6

## 
## 
## ==============================
## Variable: JUMUR 
## Error processing variable JUMUR : Assigned data `as.factor(mydata[[v]])` must be compatible with existing
## data.
## ✖ Existing data has 378 rows.
## ✖ Assigned data has 0 rows.
## ℹ Only vectors of size 1 are recycled.
## Caused by error in `vectbl_recycle_rhs_rows()`:
## ! Can't recycle input of size 0 to size 378. 
## 
## 
## ==============================
## Variable: KRT

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                         N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=bukan KRT 375       86   87.516    0.0263         5
## mydata[[v]]=KRT         3        2    0.484    4.7498         5
## 
##  Chisq= 5  on 1 degrees of freedom, p= 0.03

## 
## 
## ==============================
## Variable: DIDIK

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                             N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=SMA/Sederajat 332       71    78.78     0.769      7.69
## mydata[[v]]=D1/D2/D3       15        5     2.98     1.377      1.49
## mydata[[v]]=D4/S1          31       12     6.24     5.313      5.98
## 
##  Chisq= 7.8  on 2 degrees of freedom, p= 0.02

## 
## 
## ==============================
## Variable: MRISEN

## Call:
## survdiff(formula = mydata$s_obj ~ mydata[[v]], data = mydata)
## 
##                            N Observed Expected (O-E)^2/E (O-E)^2/V
## mydata[[v]]=Bukan Migran 374       85   87.201    0.0555       6.4
## mydata[[v]]=Migran         4        3    0.799    6.0615       6.4
## 
##  Chisq= 6.4  on 1 degrees of freedom, p= 0.01

Model Cox PH

mydata$UMUR<-as.numeric(mydata$UMUR)
cox_model <-coxph(s_obj ~ KRT + DIDIK + MRISEN, data = mydata)

summary(cox_model)

## Call:
## coxph(formula = s_obj ~ KRT + DIDIK + MRISEN, data = mydata)
## 
##   n= 378, number of events= 88 
## 
##                 coef exp(coef) se(coef)     z Pr(>|z|)  
## KRTKRT        1.6574    5.2459   0.7223 2.295   0.0217 *
## DIDIKD1/D2/D3 0.4638    1.5901   0.4930 0.941   0.3468  
## DIDIKD4/S1    0.7751    2.1708   0.3152 2.459   0.0139 *
## MRISENMigran  1.1948    3.3028   0.6263 1.908   0.0564 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## KRTKRT            5.246     0.1906    1.2735    21.609
## DIDIKD1/D2/D3     1.590     0.6289    0.6050     4.179
## DIDIKD4/S1        2.171     0.4607    1.1704     4.026
## MRISENMigran      3.303     0.3028    0.9677    11.272
## 
## Concordance= 0.552  (se = 0.022 )
## Likelihood ratio test= 12.39  on 4 df,   p=0.01
## Wald test            = 16.52  on 4 df,   p=0.002
## Score (logrank) test = 18.95  on 4 df,   p=8e-04

# Cek asumsi PH dengan cox.zph
ph_test <- cox.zph(cox_model)
print(ph_test)

##        chisq df     p
## KRT     1.89  1 0.170
## DIDIK   5.64  2 0.060
## MRISEN  2.14  1 0.143
## GLOBAL  9.81  4 0.044

# Plot grafik asumsi PH tiap variabel
plot(ph_test)

Model Stratified Cox

str_model <- coxph(s_obj ~ KRT + DIDIK + MRISEN+strata(UMUR), data = mydata)
summary(str_model)

## Call:
## coxph(formula = s_obj ~ KRT + DIDIK + MRISEN + strata(UMUR), 
##     data = mydata)
## 
##   n= 378, number of events= 88 
## 
##                 coef exp(coef) se(coef)     z Pr(>|z|)    
## KRTKRT        1.8044    6.0761   0.7489 2.409   0.0160 *  
## DIDIKD1/D2/D3 1.3817    3.9816   0.5695 2.426   0.0153 *  
## DIDIKD4/S1    1.8783    6.5426   0.4689 4.006 6.18e-05 ***
## MRISENMigran  1.3120    3.7136   0.6760 1.941   0.0523 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## KRTKRT            6.076     0.1646    1.4001     26.37
## DIDIKD1/D2/D3     3.982     0.2512    1.3041     12.16
## DIDIKD4/S1        6.543     0.1528    2.6100     16.40
## MRISENMigran      3.714     0.2693    0.9871     13.97
## 
## Concordance= 0.565  (se = 0.018 )
## Likelihood ratio test= 25.03  on 4 df,   p=5e-05
## Wald test            = 26.23  on 4 df,   p=3e-05
## Score (logrank) test = 33.32  on 4 df,   p=1e-06

ph_test_str <- cox.zph(str_model)
print(ph_test_str)

##        chisq df    p
## KRT    0.423  1 0.52
## DIDIK  1.595  2 0.45
## MRISEN 0.459  1 0.50
## GLOBAL 2.841  4 0.58

plot(ph_test_str)

AIC(str_model)

## [1] 599.8247

# Cek asumsi PH dengan cox.zph
ph_test <- cox.zph(cox_model)
print(ph_test)

##        chisq df     p
## KRT     1.89  1 0.170
## DIDIK   5.64  2 0.060
## MRISEN  2.14  1 0.143
## GLOBAL  9.81  4 0.044

# Plot grafik asumsi PH tiap variabel
plot(ph_test)

Parametrik

library(survival)
library(flexsurv)
library(eha)

# Model Exponential
exp_model <- survreg(s_obj ~ 1, data = mydata, dist = "exponential")

# Model Weibull
weib_model <- survreg(s_obj ~ 1, data = mydata, dist = "weibull")

# Model Log-normal
lnorm_model <- survreg(s_obj ~ 1, data = mydata, dist = "lognormal")

# Model Log-logistic
llogis_model <- survreg(s_obj ~ 1, data = mydata, dist = "loglogistic")

# Buat tabel AIC awal dari model parametrik
aic_table <- data.frame(
  Distribusi = c("exponential", "weibull", "lognormal", "loglogistic"),
  AIC = c(
    AIC(exp_model),
    AIC(weib_model),
    AIC(lnorm_model),
    AIC(llogis_model)
  )
)

# Tambahkan AIC dari model Cox
cox_row <- data.frame(
  Distribusi = "coxph",
  AIC = AIC(cox_model)
)

# Gabungkan semua ke dalam satu tabel
aic_table <- rbind(aic_table, cox_row)

# Urutkan berdasarkan AIC terkecil
aic_table <- aic_table[order(aic_table$AIC), ]

# Tampilkan tabel akhir
print(aic_table)

##    Distribusi      AIC
## 3   lognormal 805.1600
## 4 loglogistic 816.3117
## 2     weibull 823.2771
## 1 exponential 832.9493
## 5       coxph 938.4964

# Model survival
fit_llogis <- flexsurvreg(s_obj ~ 1, data = mydata, dist = "llogis")
fit_lnorm  <- flexsurvreg(s_obj ~ 1, data = mydata, dist = "lognormal")

# Plot survival functions
plot(fit_llogis, 
     col = "red", 
     lwd = 2, 
     lwd.ci = 0,
     xlab = "Time",
     ylab = "Survival Probability",
     main = "Survival Function: Loglogistic vs Lognormal")

# Tambahkan model lognormal ke plot
plot(fit_lnorm, 
     col = "blue", 
     lwd = 2, 
     lwd.ci = 0, 
     add = TRUE)

# Tambahkan legenda
legend("bottomleft", 
       legend = c("Loglogistic", "Lognormal"), 
       col = c("red", "blue"), 
       lty = 1, 
       lwd = 1)

Model LogNormal

library(flexsurv)

model_terbaik1 <- flexsurvreg(s_obj ~UMUR+KRT + DIDIK + MRISEN, data = mydata, dist = "lognormal")

(model_terbaik1)

## Call:
## flexsurvreg(formula = s_obj ~ UMUR + KRT + DIDIK + MRISEN, data = mydata, 
##     dist = "lognormal")
## 
## Estimates: 
##                data mean  est       L95%      U95%      se        exp(est)
## meanlog              NA    2.27312   1.58250   2.96375   0.35237        NA
## sdlog                NA    1.77995   1.51929   2.08534   0.14380        NA
## UMUR            4.31481    0.33197   0.18019   0.48374   0.07744   1.39371
## KRTKRT          0.00794   -2.18333  -4.45024   0.08358   1.15661   0.11267
## DIDIKD1/D2/D3   0.03968   -1.44553  -2.68346  -0.20759   0.63161   0.23562
## DIDIKD4/S1      0.08201   -1.75047  -2.70480  -0.79614   0.48691   0.17369
## MRISENMigran    0.01058   -1.63104  -3.51434   0.25226   0.96089   0.19573
##                L95%      U95%    
## meanlog              NA        NA
## sdlog                NA        NA
## UMUR            1.19745   1.62213
## KRTKRT          0.01168   1.08717
## DIDIKD1/D2/D3   0.06833   0.81254
## DIDIKD4/S1      0.06688   0.45106
## MRISENMigran    0.02977   1.28693
## 
## N = 378,  Events: 88,  Censored: 290
## Total time at risk: 3636
## Log-likelihood = -387.1856, df = 7
## AIC = 788.3712

AIC(model_terbaik1)

## [1] 788.3712

Model Loglogistic

library(flexsurv)

model_terbaik2 <- flexsurvreg(s_obj ~UMUR+KRT + DIDIK + MRISEN, data = mydata, dist = "llogis")

(model_terbaik2)

## Call:
## flexsurvreg(formula = s_obj ~ UMUR + KRT + DIDIK + MRISEN, data = mydata, 
##     dist = "llogis")
## 
## Estimates: 
##                data mean  est       L95%      U95%      se        exp(est)
## shape                NA    0.99422   0.83987   1.17693   0.08558        NA
## scale                NA    8.63516   4.25519  17.52354   3.11798        NA
## UMUR            4.31481    0.34966   0.18854   0.51078   0.08221   1.41858
## KRTKRT          0.00794   -2.24436  -4.21444  -0.27429   1.00516   0.10600
## DIDIKD1/D2/D3   0.03968   -1.43872  -2.69014  -0.18730   0.63849   0.23723
## DIDIKD4/S1      0.08201   -1.94019  -2.85955  -1.02082   0.46907   0.14368
## MRISENMigran    0.01058   -1.61929  -3.29242   0.05383   0.85365   0.19804
##                L95%      U95%    
## shape                NA        NA
## scale                NA        NA
## UMUR            1.20748   1.66658
## KRTKRT          0.01478   0.76011
## DIDIKD1/D2/D3   0.06787   0.82919
## DIDIKD4/S1      0.05729   0.36030
## MRISENMigran    0.03716   1.05531
## 
## N = 378,  Events: 88,  Censored: 290
## Total time at risk: 3636
## Log-likelihood = -391.2678, df = 7
## AIC = 796.5355

AIC(model_terbaik2)

## [1] 796.5355

# Buat tabel AIC dari kedua model
aic_perbandingan <- data.frame(
  Model = c("Lognormal", "Loglogistic","Cox PH (Tanpa Umur)","Stratified Cox"),
  AIC = c(AIC(model_terbaik1), AIC(model_terbaik2),AIC(cox_model),AIC(str_model))
)

# Urutkan dari AIC terkecil
aic_perbandingan <- aic_perbandingan[order(aic_perbandingan$AIC), ]

# Tampilkan tabel
print(aic_perbandingan)

##                 Model      AIC
## 4      Stratified Cox 599.8247
## 1           Lognormal 788.3712
## 2         Loglogistic 796.5355
## 3 Cox PH (Tanpa Umur) 938.4964

plot(model_terbaik2,
     type = "survival",
     col = "blue",
     lwd = 2,
     xlab = "Time",
     ylab = "Survival Probability",
     main = "Survival Curve - Loglogistik AFT Model")

plot(model_terbaik1,
     type = "survival",
     col = "blue",
     lwd = 2,
     xlab = "Time",
     ylab = "Survival Probability",
     main = "Survival Curve - Lognormal AFT Model")

summary(str_model)

## Call:
## coxph(formula = s_obj ~ KRT + DIDIK + MRISEN + strata(UMUR), 
##     data = mydata)
## 
##   n= 378, number of events= 88 
## 
##                 coef exp(coef) se(coef)     z Pr(>|z|)    
## KRTKRT        1.8044    6.0761   0.7489 2.409   0.0160 *  
## DIDIKD1/D2/D3 1.3817    3.9816   0.5695 2.426   0.0153 *  
## DIDIKD4/S1    1.8783    6.5426   0.4689 4.006 6.18e-05 ***
## MRISENMigran  1.3120    3.7136   0.6760 1.941   0.0523 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## KRTKRT            6.076     0.1646    1.4001     26.37
## DIDIKD1/D2/D3     3.982     0.2512    1.3041     12.16
## DIDIKD4/S1        6.543     0.1528    2.6100     16.40
## MRISENMigran      3.714     0.2693    0.9871     13.97
## 
## Concordance= 0.565  (se = 0.018 )
## Likelihood ratio test= 25.03  on 4 df,   p=5e-05
## Wald test            = 26.23  on 4 df,   p=3e-05
## Score (logrank) test = 33.32  on 4 df,   p=1e-06