Library

library(nnet)
library(MASS)
library(broom)
library(generalhoslem)

Dataset

Pada studi ini, digunakan data BPJS dengan jumlah baris lebih dari 1 juta. Data ini merupakan hasil kombinasi data pelayanan peserta BPJS dan data masing-masing peserta BPJS. Data di bawah dapat diakses melalui link berikut: https://ipb.link/data-kapitasi

data <- read.csv("C:/Users/Dicky Girsang/Desktop/Semester 6/Data Challenge/data_kapitasi combine.csv")

Mengubah Tipe Data

data$FKP10 <- factor(data$FKP10, levels=sort(unique(data$FKP10)))
data$FKP12 <- factor(data$FKP12, levels=sort(unique(data$FKP12)))
data$FKP13 <- factor(data$FKP13, levels=sort(unique(data$FKP13)))
data$FKP19 <- factor(data$FKP19, levels=sort(unique(data$FKP19)))

data_train <- data[,c("PSTV15.x","FKP10","FKP12", "FKP13", "FKP19")]
head(data_train)

##   PSTV15.x FKP10    FKP12           FKP13         FKP19
## 1 22.90139  RJTP PBI APBN KUNJUNGAN SEHAT TIDAK DIRUJUK
## 2 22.90139  RJTP PBI APBN KUNJUNGAN SEHAT TIDAK DIRUJUK
## 3 22.90139  RJTP PBI APBN KUNJUNGAN SEHAT TIDAK DIRUJUK
## 4 22.90139  RJTP PBI APBN KUNJUNGAN SEHAT TIDAK DIRUJUK
## 5 22.90139  RJTP PBI APBN KUNJUNGAN SEHAT TIDAK DIRUJUK
## 6 22.90139  RJTP PBI APBN KUNJUNGAN SEHAT TIDAK DIRUJUK

“Jenis segmen peserta saat peserta mengakses FKTP” -> mencirikan klaster peserta
“Status peserta di akhir masa perawatan” -> mengetahui track record peserta yang berkunjung ke FKTP
“Jenis fasilitas kesehatan tujuan rujukan” -> mengetahui tindak lanjut dari kunjungan peserta pasca kedatangan pertama

Pembagian Data Train

set.seed(12345)
indices <- sample(nrow(data))
train_prop <- 0.8
train_rows <- round(train_prop * nrow(data))

train <- data[indices[1:train_rows], ]
validasi <- data[indices[(train_rows + 1):nrow(data)], ]

Tujuan

Mengetahui faktor-faktor yang berpengaruh terhadap tingkat pelayanan FKTP (FKP10) Mengetahui karakteristik pada faktor-faktor yang memengaruhi tingkat pelayanan FKTP

Pemodelan Regresi Logistik Multinomial

model <- multinom(FKP10 ~ FKP12 + FKP13 + FKP19, data = train, weights = train$PSTV15.x)

## # weights:  42 (26 variable)
## initial  value 96745207.652179 
## iter  10 value 29120053.777244
## iter  20 value 22703884.403366
## iter  30 value 16515243.237017
## iter  40 value 16374564.403773
## iter  50 value 16351656.786595
## iter  60 value 16349444.037452
## iter  70 value 16349132.797620
## final  value 16349123.790090 
## converged

summary(model)

## Warning in sqrt(diag(vc)): NaNs produced

## Call:
## multinom(formula = FKP10 ~ FKP12 + FKP13 + FKP19, data = train, 
##     weights = train$PSTV15.x)
## 
## Coefficients:
##      (Intercept) FKP12PBI APBD FKP12PBI APBN FKP12PBPU FKP12PPU
## RITP  -27.797003     0.8469995     0.8161325 0.1941181 1.060421
## RJTP    2.384584     0.8025720     0.4881726 0.4324593 1.108456
##      FKP13KUNJUNGAN SEHAT FKP13LAIN-LAIN FKP13MENINGGAL FKP13PULANG PAKSA
## RITP           59.2445687      154.89638       409.9896          237.8886
## RJTP           -0.6476432       56.68069       343.7799         -982.6433
##      FKP13RUJUK LANJUT FKP13SEMBUH FKP19RUMAH SAKIT FKP19TIDAK DIRUJUK
## RITP          8.371893    78.05904       16.3660061        -36.1688961
## RJTP          2.019370   -12.80553       -0.2220405          0.3652155
## 
## Std. Errors:
##      (Intercept) FKP12PBI APBD FKP12PBI APBN  FKP12PBPU    FKP12PPU
## RITP 0.365609690   0.026496461   0.024091271 0.03069935 0.025297902
## RJTP 0.009259433   0.002520505   0.002034993 0.00254065 0.002256389
##      FKP13KUNJUNGAN SEHAT FKP13LAIN-LAIN FKP13MENINGGAL FKP13PULANG PAKSA
## RITP          0.913477351   8.177924e-10      0.4571496               NaN
## RJTP          0.001188116   8.790740e-10      0.4571496                 0
##      FKP13RUJUK LANJUT FKP13SEMBUH FKP19RUMAH SAKIT FKP19TIDAK DIRUJUK
## RITP        0.18285705 1.918225667       0.18286578        0.548138520
## RJTP        0.01838719 0.000158703       0.02795592        0.009234275
## 
## Residual Deviance: 32698248 
## AIC: 32698296

Uji Likelihood Ratio

Compact Model

# Compute the log-likelihood of the fitted model
llh <- logLik(model)
# Fit the null model
nullModel <- multinom(FKP10 ~ 1, data = train)

## # weights:  6 (2 variable)
## initial  value 921581.904486 
## iter  10 value 169225.204208
## final  value 168967.544488 
## converged

# Compute the log-likelihood of the null model
llhNull <- logLik(nullModel)
llhNull

## 'log Lik.' -168967.5 (df=2)

# Compute the G2 statistic
G2 <- 2 * (llh - llhNull)
G2

## 'log Lik.' -32360312 (df=24)

## bentuk model
test3 <- chisq.test(table(data_train$FKP10, data_train$FKP12))
test4 <- chisq.test(table(data_train$FKP10, data_train$FKP13))

## Warning in chisq.test(table(data_train$FKP10, data_train$FKP13)): Chi-squared
## approximation may be incorrect

test6 <- chisq.test(table(data_train$FKP10, data_train$FKP19))

## Menghitung degrees of freedom untuk model alternative
df_alt <- sum(test3$parameter, test4$parameter, test6$parameter)

## Menghitung degrees of freedom untuk model null
df_null <- 0

## Menghitung chi-square critical value pada alpha = 0.05 dengan df = 0:24
chi_crit <- qchisq(0.05, df = 0:24)

## Hasil
result <- data.frame('llh' = llh, 'llhNull' = llhNull, 'G2' = G2, 'χ2(0,05:11)' = chi_crit[df_alt+1])
print(result)

##         llh   llhNull        G2 χ2.0.05.11.
## 1 -16349124 -168967.5 -32360312    13.84843

Full Model

llFull <- logLik(model)
# Fit the null model
nullModel <- multinom(FKP10 ~ 1, data = data)

## # weights:  6 (2 variable)
## initial  value 1151977.380613 
## iter  10 value 212241.181455
## final  value 211322.212686 
## converged

# Compute the log-likelihood of the null model
llhNull <- logLik(nullModel)
llhNull

## 'log Lik.' -211322.2 (df=2)

# Compute the G2 statistic
G2 <- 2 * (llh - llhNull)
G2

## 'log Lik.' -32275603 (df=24)

## bentuk model
test5 <- chisq.test(table(data$FKP10, data$FKP05))
test7 <- chisq.test(table(data$FKP10, data$FKP07))

## Warning in chisq.test(table(data$FKP10, data$FKP07)): Chi-squared approximation
## may be incorrect

test8 <- chisq.test(table(data$FKP10, data$FKP08))
test10 <- chisq.test(table(data$FKP10, data$FKP10))
test12 <- chisq.test(table(data$FKP10, data$FKP12))
test13 <- chisq.test(table(data$FKP10, data$FKP13))

## Warning in chisq.test(table(data$FKP10, data$FKP13)): Chi-squared approximation
## may be incorrect

test19 <- chisq.test(table(data$FKP10, data$FKP19))

## Menghitung degrees of freedom untuk model alternative
df_altfull <- sum(test5$parameter,test7$parameter,test8$parameter,test10$parameter,test12$parameter,test13$parameter,test19$parameter)

## Menghitung degrees of freedom untuk model null
df_null <- 0

## Menghitung chi-square critical value pada alpha = 0.05 dengan df = 0:24
chi_crit <- qchisq(0.05, df = 0:120)

## Hasil
result <- data.frame('llh' = llFull, 'llhNull' = llhNull, 'G2' = G2, 'χ2(0,05:11)' = chi_crit[df_altfull+1])
print(result)

##         llh   llhNull        G2 χ2.0.05.11.
## 1 -16349124 -211322.2 -32275603    95.70464

Diperoleh kesimpulan H0 ditolak,artinya paling sedikit ada satu variabel prediktor yang berpengaruh terhadap variabel respon.

Uji Wald (Uji Parsial)

## Get the model coefficients and their standard errors
coef_data <- tidy(model)

## Warning in sqrt(diag(vc)): NaNs produced

View(coef_data)

## Calculate the Wald test statistics and p-values for each predictor
coef_data$wald_stat <- (coef_data$estimate / coef_data$std.error) ^ 2
coef_data$wald_pvalue <- 1 - pchisq(coef_data$wald_stat, df = 1)

result_df <- data.frame(
  "Variabel Prediktor" = rownames(coef_data),
  "Hasil Uji Wald" = paste0(round(coef_data$wald_stat, 2), " (p = ", round(coef_data$wald_pvalue, 4), ")"),
  "Kesimpulan" = ifelse(coef_data$wald_pvalue < 0.05, "Berpengaruh terhadap model", "Tidak berpengaruh terhadap model")
)
print(result_df)

##    Variabel.Prediktor               Hasil.Uji.Wald                 Kesimpulan
## 1                   1              5780.44 (p = 0) Berpengaruh terhadap model
## 2                   2              1021.86 (p = 0) Berpengaruh terhadap model
## 3                   3              1147.63 (p = 0) Berpengaruh terhadap model
## 4                   4                39.98 (p = 0) Berpengaruh terhadap model
## 5                   5              1757.06 (p = 0) Berpengaruh terhadap model
## 6                   6              4206.31 (p = 0) Berpengaruh terhadap model
## 7                   7 3.58753700700951e+22 (p = 0) Berpengaruh terhadap model
## 8                   8            804320.12 (p = 0) Berpengaruh terhadap model
## 9                   9                NaN (p = NaN)                       <NA>
## 10                 10              2096.16 (p = 0) Berpengaruh terhadap model
## 11                 11              1655.95 (p = 0) Berpengaruh terhadap model
## 12                 12              8009.78 (p = 0) Berpengaruh terhadap model
## 13                 13              4354.01 (p = 0) Berpengaruh terhadap model
## 14                 14             66321.82 (p = 0) Berpengaruh terhadap model
## 15                 15            101389.49 (p = 0) Berpengaruh terhadap model
## 16                 16             57546.78 (p = 0) Berpengaruh terhadap model
## 17                 17              28973.5 (p = 0) Berpengaruh terhadap model
## 18                 18            241328.52 (p = 0) Berpengaruh terhadap model
## 19                 19            297135.28 (p = 0) Berpengaruh terhadap model
## 20                 20 4.15737676082383e+21 (p = 0) Berpengaruh terhadap model
## 21                 21            565515.23 (p = 0) Berpengaruh terhadap model
## 22                 22                  Inf (p = 0) Berpengaruh terhadap model
## 23                 23             12061.49 (p = 0) Berpengaruh terhadap model
## 24                 24        6510657461.82 (p = 0) Berpengaruh terhadap model
## 25                 25                63.08 (p = 0) Berpengaruh terhadap model
## 26                 26               1564.2 (p = 0) Berpengaruh terhadap model

Output akhir: peubah respon, peubah prediktor per kategori, beta, Standard eror, wald-pvalue, kesimpulan.
Ada model logit RITP sama RJTP pake beta

Odds-Ratio

exp(coef(model)) # pake tabel vertikal sesuai skripsi ya

##       (Intercept) FKP12PBI APBD FKP12PBI APBN FKP12PBPU FKP12PPU
## RITP 8.470618e-13      2.332637      2.261736  1.214240 2.887586
## RJTP 1.085455e+01      2.231272      1.629336  1.541043 3.029676
##      FKP13KUNJUNGAN SEHAT FKP13LAIN-LAIN FKP13MENINGGAL FKP13PULANG PAKSA
## RITP         5.365242e+25   1.864846e+67  1.138185e+178     2.059319e+103
## RJTP         5.232776e-01   4.131520e+24  2.003195e+149      0.000000e+00
##      FKP13RUJUK LANJUT  FKP13SEMBUH FKP19RUMAH SAKIT FKP19TIDAK DIRUJUK
## RITP       4323.813092 7.954477e+33     1.281345e+07       1.959061e-16
## RJTP          7.533577 2.745546e-06     8.008829e-01       1.440825e+00

Goodness of Fit

predicted <- fitted(model)
observed <- train$FKP10
logitgof(observed,predicted, ord = FALSE)

## Warning in logitgof(observed, predicted, ord = FALSE): At least one cell in the
## expected frequencies table is < 1. Chi-square approximation may be incorrect.

## Warning in logitgof(observed, predicted, ord = FALSE): Not possible to compute
## 10 rows. There might be too few observations.

## 
##  Hosmer and Lemeshow test (multinomial model)
## 
## data:  observed, predicted
## X-squared = 1114.9, df = 8, p-value < 2.2e-16

Ketepatan Klasifikasi

# Melakukan klasifikasi menggunakan data uji untuk membentuk hasil prediksi
predict_validasi <- predict(model, newdata = validasi, type = "class")
head(predict_validasi)

## [1] RJTP RJTP RJTP RJTP RJTP RJTP
## Levels: PROMOTIF RITP RJTP

str(predict_validasi)

##  Factor w/ 3 levels "PROMOTIF","RITP",..: 3 3 3 3 3 3 3 3 3 3 ...

# confusion
table(predict_validasi, validasi$FKP10)

##                 
## predict_validasi PROMOTIF   RITP   RJTP
##         PROMOTIF        0      0      0
##         RITP            0   1037      0
##         RJTP         8464     76 200138

Kesimpulan: Dalam keseluruhan, model memiliki akurasi yang tinggi (95.9%), tetapi kesesuaian prediksi model dengan data aktual relatif rendah (nilai kappa 18.92%). Oleh karena itu, diperlukan pengembangan dan peningkatan model untuk meningkatkan kesesuaian prediksi model dengan data aktual.

Regresi Logistik Multinomial: Faktor-faktor yang Memengaruhi Tingkat Pelayanan FKTP pada Pasien Peserta BPJS

Dicky Girsang - Faadiyah Ramdani - Zaki Ramadhan - Angelika Batubara - Ervina Dwi

2023-04-11