################# LIBRARY #################
library(VIM)
## Warning: package 'VIM' was built under R version 4.3.3
## Loading required package: colorspace
## Loading required package: grid
## VIM is ready to use.
## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues
##
## Attaching package: 'VIM'
## The following object is masked from 'package:datasets':
##
## sleep
library(readxl)
library(survival)
## Warning: package 'survival' was built under R version 4.3.3
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(pastecs)
## Warning: package 'pastecs' was built under R version 4.3.3
##
## Attaching package: 'pastecs'
## The following objects are masked from 'package:dplyr':
##
## first, last
library(ggplot2)
library(vioplot)
## Warning: package 'vioplot' was built under R version 4.3.3
## Loading required package: sm
## Warning: package 'sm' was built under R version 4.3.3
## Package 'sm', version 2.2-6.0: type help(sm) for summary information
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(Hmisc)
## Warning: package 'Hmisc' was built under R version 4.3.3
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
##
## src, summarize
## The following objects are masked from 'package:base':
##
## format.pval, units
library(pec)
## Warning: package 'pec' was built under R version 4.3.3
## Loading required package: prodlim
## Warning: package 'prodlim' was built under R version 4.3.3
library(mice)
## Warning: package 'mice' was built under R version 4.3.3
##
## Attaching package: 'mice'
## The following object is masked from 'package:stats':
##
## filter
## The following objects are masked from 'package:base':
##
## cbind, rbind
library(survival)
library(survminer)
## Warning: package 'survminer' was built under R version 4.3.3
## Loading required package: ggpubr
## Warning: package 'ggpubr' was built under R version 4.3.3
##
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
##
## myeloma
library(flexsurv)
## Warning: package 'flexsurv' was built under R version 4.3.3
library(parfm)
## Warning: package 'parfm' was built under R version 4.3.3
## Loading required package: optimx
## Warning: package 'optimx' was built under R version 4.3.3
library(vioplot)
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:sm':
##
## muscle
## The following object is masked from 'package:dplyr':
##
## select
################# LOAD DATA #################
setwd("C:/Users/Lenovo/Downloads/SKRIPSI/Syntax")
data <- read_excel("DataFix.xlsx")
data
## # A tibble: 212 × 11
## i t event usia sist dias jk komp dm imt obat
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 303 0 66 89 161 1 1 1 26.6 1
## 2 2 353 0 74 90 158 1 1 0 24.3 1
## 3 3 37 1 77 70 140 1 0 0 31.2 0
## 4 4 327 0 54 89 171 1 0 0 NA 0
## 5 5 170 0 85 72 142 0 1 1 NA 0
## 6 6 348 0 85 91 136 0 1 0 27.1 1
## 7 7 53 0 89 91 148 1 1 0 20.4 1
## 8 8 172 0 50 77 130 0 1 0 25.4 1
## 9 9 101 0 66 65 101 0 0 0 NA 1
## 10 10 156 0 58 87 115 0 1 0 20.0 1
## # ℹ 202 more rows
################# CHECKING MISSING VALUE #################
sum(is.na(data))
## [1] 74
data_miss = aggr(data, numbers=TRUE, prop = c(TRUE, FALSE),
sortVars=TRUE, labels=names(data), cex.axis=.7, gap=3,
ylab=c("Proporsi Data Hilang Tiap Kovariat","Pola Data Hilang Dalam Data"))
##
## Variables sorted by number of missings:
## Variable Count
## imt 0.3490566
## i 0.0000000
## t 0.0000000
## event 0.0000000
## usia 0.0000000
## sist 0.0000000
## dias 0.0000000
## jk 0.0000000
## komp 0.0000000
## dm 0.0000000
## obat 0.0000000
################# DATA IMPUTATION #################
i = data$i
waktu =data$t
event = data$event
usia = data$usia
sist = data$sist
dias = data$dias
imt = data$imt
jk = factor(data$jk)
komp = factor(data$komp)
dm = factor(data$dm)
obat = factor(data$obat)
datamiss <- data.frame(i, waktu, event, usia, sist, dias, jk, komp, dm, imt, obat)
colSums(is.na(datamiss))
## i waktu event usia sist dias jk komp dm imt obat
## 0 0 0 0 0 0 0 0 0 74 0
# Run MICE with Predictive Mean Matching (PMM)
set.seed(12345)
imp <- mice(datamiss, m = 5, method = "pmm", maxit = 10, seed = 12345)
##
## iter imp variable
## 1 1 imt
## 1 2 imt
## 1 3 imt
## 1 4 imt
## 1 5 imt
## 2 1 imt
## 2 2 imt
## 2 3 imt
## 2 4 imt
## 2 5 imt
## 3 1 imt
## 3 2 imt
## 3 3 imt
## 3 4 imt
## 3 5 imt
## 4 1 imt
## 4 2 imt
## 4 3 imt
## 4 4 imt
## 4 5 imt
## 5 1 imt
## 5 2 imt
## 5 3 imt
## 5 4 imt
## 5 5 imt
## 6 1 imt
## 6 2 imt
## 6 3 imt
## 6 4 imt
## 6 5 imt
## 7 1 imt
## 7 2 imt
## 7 3 imt
## 7 4 imt
## 7 5 imt
## 8 1 imt
## 8 2 imt
## 8 3 imt
## 8 4 imt
## 8 5 imt
## 9 1 imt
## 9 2 imt
## 9 3 imt
## 9 4 imt
## 9 5 imt
## 10 1 imt
## 10 2 imt
## 10 3 imt
## 10 4 imt
## 10 5 imt
summary(imp)
## Class: mids
## Number of multiple imputations: 5
## Imputation methods:
## i waktu event usia sist dias jk komp dm imt obat
## "" "" "" "" "" "" "" "" "" "pmm" ""
## PredictorMatrix:
## i waktu event usia sist dias jk komp dm imt obat
## i 0 1 1 1 1 1 1 1 1 1 1
## waktu 1 0 1 1 1 1 1 1 1 1 1
## event 1 1 0 1 1 1 1 1 1 1 1
## usia 1 1 1 0 1 1 1 1 1 1 1
## sist 1 1 1 1 0 1 1 1 1 1 1
## dias 1 1 1 1 1 0 1 1 1 1 1
# Show the imputation result of variable 'imt'
imp$imp$imt
## 1 2 3 4 5
## 4 23.42209 33.29865 28.00000 27.47138 24.22145
## 5 27.05515 16.01562 16.49396 24.16327 16.01562
## 9 20.02884 16.49396 21.08281 27.53123 20.44444
## 11 25.39062 23.30905 23.30905 27.34375 23.30905
## 14 23.87511 20.06920 28.51562 20.06920 30.62925
## 18 29.71429 23.30905 21.08281 23.30905 25.15315
## 20 18.82711 45.26935 26.70940 35.67182 24.22145
## 29 44.29630 24.22145 28.51562 26.70940 33.29865
## 30 23.30905 19.90997 35.20430 23.30905 23.30905
## 35 28.51562 22.89282 26.70940 29.40227 23.30905
## 38 26.66667 24.97399 26.03749 25.39062 21.41094
## 48 22.86237 24.91990 27.34375 21.41094 23.49524
## 49 24.22145 33.16327 20.06920 23.30905 28.00000
## 50 45.26935 19.77238 27.53123 22.89282 30.22222
## 53 28.51562 33.29865 28.51562 28.51562 18.82711
## 65 26.66667 20.06920 28.51562 18.82711 20.06920
## 66 22.60026 12.48699 23.30905 24.97399 25.33333
## 68 22.54658 26.70940 26.03749 23.43750 23.30905
## 70 34.13111 23.49524 23.30905 23.30905 23.42209
## 82 16.49396 23.30905 16.49396 27.05515 16.01562
## 84 28.35306 24.80159 20.02884 21.08281 23.30905
## 86 33.29865 23.49524 23.30905 24.55775 25.91068
## 89 18.82711 25.39062 23.30905 33.20312 22.60026
## 92 27.26801 16.01562 26.83865 27.05515 21.08281
## 97 24.80159 26.03749 20.06920 24.80159 22.89282
## 100 30.22222 26.03749 26.66667 45.26935 26.43807
## 104 16.49396 23.42209 24.34176 27.39226 23.30905
## 107 25.20398 26.03749 33.20312 24.97399 26.66667
## 109 23.30905 23.30905 33.16327 35.67182 25.06575
## 110 23.30905 23.30905 22.49135 25.22409 23.30905
## 116 24.34176 20.02884 26.56250 27.34375 25.39062
## 118 23.87511 22.89282 23.30905 23.93606 17.08744
## 121 34.96358 25.80645 21.08281 22.89282 25.39062
## 126 22.89282 23.43750 26.70940 33.29865 23.30905
## 130 23.42209 27.53123 21.08281 21.08281 30.49353
## 132 19.77238 41.09139 44.29630 30.62925 45.26935
## 133 24.80159 18.13202 23.30905 22.89282 23.30905
## 134 30.49353 27.77427 28.00000 26.70940 24.55775
## 135 27.39226 23.30905 23.30905 22.89282 24.60973
## 136 33.69494 21.64412 28.00000 23.42209 22.60026
## 140 23.30905 23.30905 34.13111 18.82711 26.43807
## 142 23.30905 23.30905 23.30905 23.30905 26.03749
## 145 45.26935 25.06575 28.00000 24.80159 20.06920
## 147 23.87511 23.30905 27.05515 25.33333 22.49135
## 151 22.89282 23.30905 18.82711 18.82711 44.29630
## 154 23.30905 25.29938 25.15315 29.40227 43.70447
## 158 24.34176 23.30905 24.60973 16.84747 24.88889
## 161 25.22409 26.43807 22.49135 45.26935 35.67182
## 162 30.62925 24.97399 23.30905 24.03441 23.30905
## 163 26.50212 33.16327 29.17488 23.43750 20.06920
## 164 18.13202 23.49524 21.08281 22.89282 24.55775
## 172 16.01562 22.30815 26.66667 18.82711 43.70447
## 175 23.42209 25.29938 23.87511 25.22409 23.30905
## 176 23.30905 30.49353 21.08281 33.20312 23.30905
## 179 33.29865 28.51562 28.51562 34.13111 24.22145
## 180 23.30905 26.56250 16.01562 23.30905 23.87511
## 181 26.70940 23.30905 25.39062 35.67182 26.03749
## 183 20.31250 29.71429 27.05515 22.60026 34.96358
## 184 45.26935 22.60026 25.20398 23.30905 17.08744
## 188 26.50212 28.51562 33.29865 33.29865 24.60973
## 189 22.54658 23.93606 24.60973 24.03441 41.09139
## 190 23.30905 24.60973 27.53123 23.30905 23.30905
## 191 41.09139 24.14152 25.80645 23.93606 22.82996
## 192 23.30905 27.47138 24.22145 18.82711 24.22145
## 193 22.89282 25.71166 16.49396 23.43750 25.39062
## 196 26.56250 23.30905 23.30905 24.22145 26.66667
## 200 27.77427 27.26801 25.71166 22.49135 27.34375
## 202 23.30905 24.60973 23.30905 24.22145 41.09139
## 203 23.30905 24.22145 20.06920 29.40227 33.29865
## 204 24.22145 24.55775 18.13202 23.30905 23.30905
## 205 27.05515 25.20398 27.05515 22.54658 23.30905
## 207 23.30905 22.49135 23.30905 16.49396 26.70940
## 208 18.82711 21.41094 23.30905 43.70447 26.66667
## 209 18.82711 26.66667 30.62925 22.54658 44.29630
# Dataset : imputation result
data_imputed <- complete(imp)
# Checking missing value after imputation
colSums(is.na(data_imputed))
## i waktu event usia sist dias jk komp dm imt obat
## 0 0 0 0 0 0 0 0 0 0 0
data_imputed
## i waktu event usia sist dias jk komp dm imt obat
## 1 1 303 0 66 89 161 1 1 1 26.56250 1
## 2 2 353 0 74 90 158 1 1 0 24.34176 1
## 3 3 37 1 77 70 140 1 0 0 31.21748 0
## 4 4 327 0 54 89 171 1 0 0 23.42209 0
## 5 5 170 0 85 72 142 0 1 1 27.05515 0
## 6 6 348 0 85 91 136 0 1 0 27.05515 1
## 7 7 53 0 89 91 148 1 1 0 20.44444 1
## 8 8 172 0 50 77 130 0 1 0 25.39062 1
## 9 9 101 0 66 65 101 0 0 0 20.02884 1
## 10 10 156 0 58 87 115 0 1 0 20.02884 1
## 11 11 215 1 74 65 147 1 1 1 25.39062 0
## 12 12 45 0 58 77 128 1 1 1 27.05515 1
## 13 13 192 0 51 80 135 0 0 0 27.34375 0
## 14 14 30 1 56 109 179 0 0 1 23.87511 0
## 15 15 279 0 68 68 148 1 1 1 24.91990 0
## 16 16 100 0 81 80 150 1 1 0 16.49396 1
## 17 17 67 1 49 80 130 1 1 1 41.09139 0
## 18 18 146 0 78 74 146 1 1 0 29.71429 0
## 19 19 94 1 77 59 105 1 1 1 19.90997 0
## 20 20 181 0 42 80 156 1 0 0 18.82711 0
## 21 21 171 0 61 71 168 1 1 1 24.88889 1
## 22 22 30 1 73 80 130 1 1 1 21.64412 0
## 23 23 223 0 78 56 125 1 1 0 23.30905 1
## 24 24 172 0 46 70 130 1 0 0 26.66667 0
## 25 25 115 0 58 90 149 1 0 0 24.14152 0
## 26 26 91 1 74 60 130 0 1 0 16.84747 1
## 27 27 166 0 59 94 178 1 0 0 25.29938 0
## 28 28 149 0 87 76 171 1 1 0 21.08281 1
## 29 29 299 0 39 94 153 1 0 0 44.29630 0
## 30 30 83 0 71 96 180 1 0 0 23.30905 0
## 31 31 207 1 80 90 140 1 1 1 26.66667 1
## 32 32 146 0 69 69 144 1 1 1 20.31250 1
## 33 33 30 1 70 77 130 1 1 1 24.97399 1
## 34 34 320 0 69 90 130 1 1 1 25.15315 1
## 35 35 69 0 49 99 170 1 0 0 28.51562 0
## 36 36 191 0 64 87 184 0 1 1 25.39062 1
## 37 37 84 1 71 80 151 0 1 1 23.30905 0
## 38 38 143 0 64 90 160 0 1 0 26.66667 0
## 39 39 167 0 60 93 153 1 0 0 26.66667 0
## 40 40 200 0 51 60 140 1 0 1 28.35306 0
## 41 41 35 1 37 89 153 1 0 0 33.29865 0
## 42 42 359 0 71 53 134 0 1 1 34.89440 1
## 43 43 314 0 51 83 152 1 1 0 26.03749 0
## 44 44 135 0 58 65 145 0 1 1 23.87511 1
## 45 45 7 1 44 80 140 1 1 1 22.54658 0
## 46 46 200 0 56 89 130 1 0 0 23.30905 0
## 47 47 32 1 70 60 150 1 1 1 24.03441 0
## 48 48 216 0 68 96 177 0 1 0 22.86237 0
## 49 49 138 0 39 81 163 1 1 0 24.22145 0
## 50 50 72 0 66 80 130 1 1 1 45.26935 1
## 51 51 38 1 75 92 176 0 1 1 33.69494 1
## 52 52 19 1 69 90 150 0 1 1 21.41094 1
## 53 53 13 1 37 93 160 1 0 1 28.51562 0
## 54 54 171 0 53 79 140 1 1 1 24.88889 1
## 55 55 32 1 58 90 140 0 1 0 34.96358 0
## 56 56 278 0 40 78 179 1 1 0 23.30905 1
## 57 57 310 0 77 70 170 1 1 0 12.48699 1
## 58 58 311 0 38 74 138 1 0 0 33.29865 0
## 59 59 83 1 65 80 140 1 0 0 23.82812 1
## 60 60 130 1 49 109 150 1 1 0 23.30905 0
## 61 61 331 0 28 77 146 1 0 0 20.06920 0
## 62 62 38 1 54 79 128 1 1 0 27.34375 0
## 63 63 362 0 76 69 124 0 1 1 22.89282 1
## 64 64 305 0 77 78 130 1 1 1 25.71166 1
## 65 65 317 0 44 117 180 0 0 0 26.66667 0
## 66 66 178 0 64 67 108 1 0 0 22.60026 0
## 67 67 156 0 66 91 190 1 0 0 28.72008 0
## 68 68 101 0 57 89 149 1 0 0 22.54658 0
## 69 69 321 0 37 80 150 0 0 0 24.22145 0
## 70 70 199 0 55 90 140 1 1 1 34.13111 0
## 71 71 32 1 70 70 120 1 1 0 25.20398 0
## 72 72 175 1 70 80 150 1 1 0 30.22222 0
## 73 73 166 0 62 103 158 0 0 0 20.76125 0
## 74 74 282 0 50 70 141 1 1 0 23.30905 1
## 75 75 296 0 70 60 134 1 1 0 30.22222 1
## 76 76 18 1 31 90 137 0 1 0 24.22145 0
## 77 77 167 1 66 70 140 1 1 1 29.17488 1
## 78 78 228 0 35 100 140 1 1 0 23.93606 0
## 79 79 338 0 56 90 168 0 1 1 23.87511 1
## 80 80 314 0 75 69 163 1 1 1 27.53123 1
## 81 81 355 0 80 65 125 1 1 1 23.30905 1
## 82 82 311 0 86 67 157 0 1 0 16.49396 1
## 83 83 340 0 74 66 145 1 1 1 23.30905 0
## 84 84 327 0 54 62 133 1 1 0 28.35306 1
## 85 85 114 0 49 101 140 0 0 0 23.87511 0
## 86 86 115 0 60 89 149 1 0 0 33.29865 0
## 87 87 63 1 63 75 119 0 1 0 22.30815 0
## 88 88 60 1 64 85 125 1 1 1 23.30905 0
## 89 89 45 0 55 76 157 1 0 0 18.82711 0
## 90 90 160 1 46 89 121 1 1 1 34.13111 1
## 91 91 341 0 39 98 146 1 1 0 33.20312 1
## 92 92 131 0 62 92 130 0 1 0 27.26801 1
## 93 93 99 1 80 90 170 0 1 1 26.83865 1
## 94 94 306 0 57 91 157 0 1 0 23.30905 1
## 95 95 18 0 54 91 138 1 1 0 25.80645 1
## 96 96 23 1 68 77 146 0 1 1 23.30905 0
## 97 97 132 0 43 109 151 0 1 0 24.80159 0
## 98 98 37 0 25 78 127 1 0 0 22.76944 0
## 99 99 354 0 64 73 151 0 1 1 23.30905 1
## 100 100 314 0 59 94 154 0 0 0 30.22222 1
## 101 101 66 0 51 94 152 1 1 0 30.49353 1
## 102 102 201 0 57 97 161 1 1 1 23.30905 1
## 103 103 355 0 71 91 143 0 1 0 22.49135 1
## 104 104 190 1 62 85 136 0 1 0 16.49396 0
## 105 105 24 0 64 86 150 1 0 0 23.49524 0
## 106 106 131 1 75 71 169 1 1 0 23.30905 0
## 107 107 139 0 52 90 150 0 1 0 25.20398 0
## 108 108 75 1 52 78 127 1 1 1 23.30905 1
## 109 109 311 0 62 104 178 0 1 1 23.30905 0
## 110 110 18 0 77 79 165 0 0 0 23.30905 1
## 111 111 187 1 47 85 130 1 1 0 20.70312 1
## 112 112 86 0 59 84 145 1 1 0 27.26801 1
## 113 113 160 0 63 63 141 0 1 1 23.30905 1
## 114 114 27 1 31 100 140 0 1 1 29.40227 1
## 115 115 30 1 86 57 132 0 1 1 22.86237 0
## 116 116 151 0 66 90 157 0 0 0 24.34176 1
## 117 117 321 0 52 93 192 1 0 0 33.16327 0
## 118 118 332 0 61 102 155 1 0 0 23.87511 0
## 119 119 59 0 62 70 140 1 1 0 26.70940 1
## 120 120 348 0 43 66 116 0 1 1 18.13202 1
## 121 121 36 1 78 74 148 1 1 0 34.96358 0
## 122 122 258 0 40 121 158 1 0 0 32.45568 0
## 123 123 137 0 52 91 168 1 0 0 26.50212 0
## 124 124 182 1 43 113 190 1 1 1 26.70940 1
## 125 125 142 0 39 108 170 1 0 0 23.30905 0
## 126 126 342 0 55 80 150 0 0 0 22.89282 0
## 127 127 118 0 52 87 187 0 1 0 26.03749 0
## 128 128 208 0 65 113 190 1 0 0 26.70940 1
## 129 129 174 1 63 65 125 1 1 1 23.30905 0
## 130 130 151 0 66 90 140 1 1 0 23.42209 0
## 131 131 170 0 51 92 149 0 0 0 23.30905 0
## 132 132 362 0 47 97 154 0 1 1 19.77238 1
## 133 133 356 0 64 81 148 1 1 0 24.80159 0
## 134 134 200 0 35 90 130 0 1 0 30.49353 0
## 135 135 327 0 70 85 156 1 1 0 27.39226 0
## 136 136 83 0 60 97 186 0 1 1 33.69494 1
## 137 137 156 0 61 86 133 0 0 0 24.97399 0
## 138 138 7 1 48 84 159 1 1 0 45.26935 0
## 139 139 93 1 42 111 160 1 1 1 18.82711 0
## 140 140 44 0 47 100 147 0 0 0 23.30905 0
## 141 141 363 0 63 53 122 1 1 1 25.33333 1
## 142 142 339 0 59 92 153 1 1 0 23.30905 1
## 143 143 31 1 60 70 140 1 0 1 26.43807 1
## 144 144 282 0 69 90 160 0 1 0 31.11111 1
## 145 145 335 0 41 93 171 1 1 0 45.26935 0
## 146 146 11 0 23 90 140 0 1 0 27.77427 1
## 147 147 3 0 90 63 147 1 0 0 23.87511 0
## 148 148 263 0 59 58 139 1 0 0 27.39226 0
## 149 149 361 0 61 100 180 1 0 0 22.60026 0
## 150 150 116 0 62 73 149 1 1 1 23.30905 1
## 151 151 103 0 29 96 153 0 0 0 22.89282 0
## 152 152 319 0 41 87 145 0 1 0 43.70447 0
## 153 153 202 0 37 103 154 0 1 0 23.30905 0
## 154 154 345 0 54 83 136 0 0 0 23.30905 0
## 155 155 293 0 54 73 236 1 1 1 23.30905 1
## 156 156 333 0 43 136 180 0 1 0 24.60973 1
## 157 157 325 0 34 88 143 1 0 0 17.08744 0
## 158 158 171 0 45 89 147 0 1 0 24.34176 1
## 159 159 178 0 70 88 156 0 0 0 25.39062 0
## 160 160 31 1 40 99 160 1 1 0 25.06575 0
## 161 161 293 0 66 88 171 1 0 0 25.22409 0
## 162 162 60 0 69 85 181 1 1 1 30.62925 0
## 163 163 291 0 52 81 158 1 1 1 26.50212 0
## 164 164 292 0 66 80 145 1 1 0 18.13202 0
## 165 165 206 0 57 68 163 0 1 0 24.16327 1
## 166 166 133 1 59 60 148 0 1 1 29.71429 0
## 167 167 153 0 32 100 136 1 0 0 23.42209 1
## 168 168 221 0 24 140 181 1 0 0 28.51562 0
## 169 169 12 0 63 93 167 1 1 0 21.28743 1
## 170 170 209 0 52 90 160 0 1 0 25.91068 0
## 171 171 171 1 52 80 150 1 1 1 22.82996 0
## 172 172 202 0 64 90 168 0 0 0 16.01562 0
## 173 173 202 0 35 108 156 1 0 0 44.29630 0
## 174 174 1 1 52 86 130 0 1 0 26.03749 1
## 175 175 209 0 52 90 140 1 1 1 23.42209 0
## 176 176 234 0 64 65 156 1 0 0 23.30905 0
## 177 177 112 1 60 95 171 1 0 1 24.55775 1
## 178 178 156 0 48 68 129 0 1 0 24.80159 0
## 179 179 126 1 22 109 143 1 1 0 33.29865 0
## 180 180 207 0 67 87 159 0 1 0 23.30905 0
## 181 181 200 0 52 80 140 0 0 0 26.70940 0
## 182 182 101 1 55 85 145 1 1 0 23.43750 0
## 183 183 199 0 43 50 120 0 1 1 20.31250 0
## 184 184 178 0 52 90 160 1 0 0 45.26935 0
## 185 185 171 0 54 88 167 1 0 0 23.49524 0
## 186 186 138 0 49 110 234 1 1 1 35.67182 0
## 187 187 34 1 71 83 151 1 1 0 24.34176 0
## 188 188 4 1 25 64 134 0 0 1 26.50212 1
## 189 189 33 1 55 107 173 1 1 0 22.54658 0
## 190 190 92 1 71 92 168 1 1 1 23.30905 0
## 191 191 142 0 54 100 156 1 0 0 41.09139 0
## 192 192 142 0 36 100 180 0 0 0 23.30905 0
## 193 193 122 0 68 90 130 0 0 0 22.89282 0
## 194 194 96 0 35 108 165 1 0 0 28.00000 0
## 195 195 129 0 27 86 145 1 0 0 27.47138 0
## 196 196 128 0 59 107 170 0 0 0 26.56250 1
## 197 197 30 0 41 82 142 1 0 0 25.22409 1
## 198 198 115 0 50 104 162 1 0 0 22.89282 0
## 199 199 118 0 76 82 146 1 0 0 19.77238 0
## 200 200 115 0 67 100 170 1 1 0 27.77427 0
## 201 201 109 0 81 90 182 0 1 0 16.01562 0
## 202 202 100 0 42 90 165 0 0 0 23.30905 0
## 203 203 103 0 26 124 199 0 1 0 23.30905 0
## 204 204 69 0 33 82 128 1 0 0 24.22145 0
## 205 205 67 0 56 92 154 0 0 0 27.05515 0
## 206 206 3 0 59 69 139 0 0 0 26.56250 0
## 207 207 20 0 60 80 130 0 1 0 23.30905 0
## 208 208 20 0 51 93 146 1 0 0 18.82711 0
## 209 209 10 0 37 106 155 1 0 0 18.82711 0
## 210 210 12 0 47 85 135 1 0 0 35.20430 0
## 211 211 17 0 28 83 123 1 0 0 30.62925 0
## 212 212 4 0 55 79 145 1 0 0 23.30905 1
# Save imputation result
#saveRDS(data_imputed, file = "C:/Users/Lenovo/Downloads/SKRIPSI/Syntax/data_imputed.rds")
#options(max.print = 10000)
#print(data_imputed)
# Testing the mean of imt before and after imputation
# Data before imputation (contains missing values)
imt_before <- datamiss$imt
# Data after imputation (without missing value)
imt_after <- data_imputed$imt
# Check the distribution of IMT values before and after imputation
summary(imt_before)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 12.49 23.31 24.71 25.66 27.21 45.27 74
summary(imt_after)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 12.49 23.31 24.34 25.77 27.29 45.27
# Boxplot
df <- data.frame(
Group = rep(c("Sebelum Imputasi", "Setelah Imputasi"),
c(length(imt_before), length(imt_after))),
IMT = c(imt_before, imt_after))
ggplot(df, aes(x = Group, y = IMT, fill = Group)) +
geom_boxplot() +
theme_minimal() +
labs(title = "Distribusi IMT Sebelum dan Setelah Imputasi",
x = "Kelompok",
y = "IMT")
## Warning: Removed 74 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
# Normality test
shapiro.test(imt_before)
##
## Shapiro-Wilk normality test
##
## data: imt_before
## W = 0.89016, p-value = 1.135e-08
shapiro.test(imt_after)
##
## Shapiro-Wilk normality test
##
## data: imt_after
## W = 0.86853, p-value = 1.422e-12
# Use wilcoxon because data is not normal
wilcox.test(imt_before, imt_after, paired = FALSE)
##
## Wilcoxon rank sum test with continuity correction
##
## data: imt_before and imt_after
## W = 14790, p-value = 0.8606
## alternative hypothesis: true location shift is not equal to 0
## Karena p-val=0.8606>0.05, tidak ada perbedaan signifikan antara imt_before dan imt_after.
################# DESCRIPTIVE STATISTICS #################
data = data_imputed
par(mfrow=c(1,1))
#Categorized Value
data$jk = factor(data$jk)
data$komp = factor(data$komp)
data$dm = factor(data$dm)
data$obat = factor(data$obat)
data$event = factor(data$event)
summary(data$t)
## Length Class Mode
## 0 NULL NULL
summary(data$event)
## 0 1
## 160 52
summary(data$usia)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 22.00 47.75 58.00 56.85 67.00 90.00
summary(data$sist)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 50.00 76.75 87.00 85.24 92.25 140.00
summary(data$dias)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 101.0 138.8 149.0 150.5 160.2 236.0
summary(data$imt)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 12.49 23.31 24.34 25.77 27.29 45.27
summary_by_event <- aggregate(usia ~ event, data = data, FUN = function(x) {
c(
Mean = mean(x),
Median = median(x),
SD = sd(x),
Min = min(x),
Max = max(x),
Count = length(x)
)
})
summary_by_event
## event usia.Mean usia.Median usia.SD usia.Min usia.Max usia.Count
## 1 0 56.16875 57.00000 14.45758 23.00000 90.00000 160.00000
## 2 1 58.96154 61.00000 15.38265 22.00000 86.00000 52.00000
################# DATA CHECKING #################
boxplot(data[, c("usia", "sist", "dias", "imt")], main="Boxplot Variabel Numerik")
# check outlier
check_outliers <- function(x) {
Q1 <- quantile(x, 0.25, na.rm = TRUE)
Q3 <- quantile(x, 0.75, na.rm = TRUE)
IQR_value <- Q3 - Q1
lower_bound <- Q1 - 1.5 * IQR_value
upper_bound <- Q3 + 1.5 * IQR_value
sum(x < lower_bound | x > upper_bound, na.rm = TRUE)
}
sapply(data[, c("usia", "sist", "dias", "imt")], check_outliers)
## usia sist dias imt
## 0 8 5 30
# Data 'IMT'
Q1 <- quantile(data$imt, 0.25)
Q3 <- quantile(data$imt, 0.75)
IQR <- Q3 - Q1
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
which(data$imt < lower_bound | data$imt > upper_bound)
## [1] 16 17 26 29 41 42 50 51 55 57 58 70 82 86 90 104 121 136 138
## [20] 145 152 157 172 173 179 184 186 191 201 210
data[data$imt < lower_bound | data$imt > upper_bound, ]
## i waktu event usia sist dias jk komp dm imt obat
## 16 16 100 0 81 80 150 1 1 0 16.49396 1
## 17 17 67 1 49 80 130 1 1 1 41.09139 0
## 26 26 91 1 74 60 130 0 1 0 16.84747 1
## 29 29 299 0 39 94 153 1 0 0 44.29630 0
## 41 41 35 1 37 89 153 1 0 0 33.29865 0
## 42 42 359 0 71 53 134 0 1 1 34.89440 1
## 50 50 72 0 66 80 130 1 1 1 45.26935 1
## 51 51 38 1 75 92 176 0 1 1 33.69494 1
## 55 55 32 1 58 90 140 0 1 0 34.96358 0
## 57 57 310 0 77 70 170 1 1 0 12.48699 1
## 58 58 311 0 38 74 138 1 0 0 33.29865 0
## 70 70 199 0 55 90 140 1 1 1 34.13111 0
## 82 82 311 0 86 67 157 0 1 0 16.49396 1
## 86 86 115 0 60 89 149 1 0 0 33.29865 0
## 90 90 160 1 46 89 121 1 1 1 34.13111 1
## 104 104 190 1 62 85 136 0 1 0 16.49396 0
## 121 121 36 1 78 74 148 1 1 0 34.96358 0
## 136 136 83 0 60 97 186 0 1 1 33.69494 1
## 138 138 7 1 48 84 159 1 1 0 45.26935 0
## 145 145 335 0 41 93 171 1 1 0 45.26935 0
## 152 152 319 0 41 87 145 0 1 0 43.70447 0
## 157 157 325 0 34 88 143 1 0 0 17.08744 0
## 172 172 202 0 64 90 168 0 0 0 16.01562 0
## 173 173 202 0 35 108 156 1 0 0 44.29630 0
## 179 179 126 1 22 109 143 1 1 0 33.29865 0
## 184 184 178 0 52 90 160 1 0 0 45.26935 0
## 186 186 138 0 49 110 234 1 1 1 35.67182 0
## 191 191 142 0 54 100 156 1 0 0 41.09139 0
## 201 201 109 0 81 90 182 0 1 0 16.01562 0
## 210 210 12 0 47 85 135 1 0 0 35.20430 0
which(data$imt < 19 | data$imt > 31)
## [1] 3 16 17 20 26 29 41 42 50 51 55 57 58 70 82 86 89 90 91
## [20] 104 117 120 121 122 136 138 139 144 145 152 157 164 172 173 179 184 186 191
## [39] 201 208 209 210
# Data 'SIST'
Q1 <- quantile(data$sist, 0.25)
Q3 <- quantile(data$sist, 0.75)
IQR <- Q3 - Q1
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
which(data$sist < lower_bound | data$sist > upper_bound)
## [1] 42 65 122 141 156 168 183 203
data[data$sist < lower_bound | data$sist > upper_bound, ]
## i waktu event usia sist dias jk komp dm imt obat
## 42 42 359 0 71 53 134 0 1 1 34.89440 1
## 65 65 317 0 44 117 180 0 0 0 26.66667 0
## 122 122 258 0 40 121 158 1 0 0 32.45568 0
## 141 141 363 0 63 53 122 1 1 1 25.33333 1
## 156 156 333 0 43 136 180 0 1 0 24.60973 1
## 168 168 221 0 24 140 181 1 0 0 28.51562 0
## 183 183 199 0 43 50 120 0 1 1 20.31250 0
## 203 203 103 0 26 124 199 0 1 0 23.30905 0
# Data 'DIAS'
Q1 <- quantile(data$dias, 0.25)
Q3 <- quantile(data$dias, 0.75)
IQR <- Q3 - Q1
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
which(data$dias < lower_bound | data$dias > upper_bound)
## [1] 9 19 155 186 203
data[data$dias < lower_bound | data$dias > upper_bound, ]
## i waktu event usia sist dias jk komp dm imt obat
## 9 9 101 0 66 65 101 0 0 0 20.02884 1
## 19 19 94 1 77 59 105 1 1 1 19.90997 0
## 155 155 293 0 54 73 236 1 1 1 23.30905 1
## 186 186 138 0 49 110 234 1 1 1 35.67182 0
## 203 203 103 0 26 124 199 0 1 0 23.30905 0
#Hapus outlier
#rows_to_remove <- c(3,9,16,17,19,26,29,41,42,50,51,55,57,58,65,70,82,86,89,90,91,104,117,120,121,122,136,138,139,144,141,145,152,155,91,156,157,164,168,172,173,179,183,184,186,191,201,203,210,208,209)
#data_new <- data[-unique(rows_to_remove), ]
#data_new
# check corelation
#cor(data[, c("usia", "sist", "dias", "imt")])
# check event proportion
#table(data$event)
#prop.table(table(data$event))
################# DATA VISUALISATION #################
# Boxplot untuk usia
ggplot(data, aes(x = factor(event), y = usia, fill = factor(event))) +
geom_boxplot() +
ylim(0, 100) +
labs(title = "USIA", x = "Status", y = "Usia") +
scale_fill_manual(values = c("lightpink", "violet")) +
theme_minimal() +
theme(legend.position = "none") +
scale_x_discrete(labels = c("0" = "Sensor", "1" = "Event"))
# Boxplot untuk sistolik
ggplot(data, aes(x = factor(event), y = sist, fill = factor(event))) +
geom_boxplot() +
ylim(0, 150) +
labs(title = "SISTOLIK", x = "Status", y = "Sistolik") +
scale_fill_manual(values = c("lightpink", "violet")) +
theme_minimal() +
theme(legend.position = "none") +
scale_x_discrete(labels = c("0" = "Sensor", "1" = "Event"))
# Boxplot untuk diastolik
ggplot(data, aes(x = factor(event), y = dias, fill = factor(event))) +
geom_boxplot() +
ylim(0, 250) +
labs(title = "DIASTOLIK", x = "Status", y = "Diastolik") +
scale_fill_manual(values = c("lightpink", "violet")) +
theme_minimal() +
theme(legend.position = "none") +
scale_x_discrete(labels = c("0" = "Sensor", "1" = "Event"))
# Boxplot untuk IMT
ggplot(data, aes(x = factor(event), y = imt, fill = factor(event))) +
geom_boxplot() +
ylim(0, 50) +
labs(title = "IMT", x = "Status", y = "IMT") +
scale_fill_manual(values = c("lightpink", "violet")) +
theme_minimal() +
theme(legend.position = "none") + # Menghilangkan legend
scale_x_discrete(labels = c("0" = "Sensor", "1" = "Event"))
# Stacked Bar Chart : Jenis Kelamin
ggplot(data, aes(x=jk, fill=event))+
geom_bar(stat="count", position="stack", width=0.7)+
geom_text(aes(label = after_stat(count)), stat="count", vjust=-0.5, color="grey17", size=4,
fontface="bold", position=position_stack(vjust = 0.5))+
labs(fill= "Status", y="Jumlah Pasien")+
scale_fill_manual(values=c("0"="salmon1", "1"="turquoise3"),
labels=c("0"="Sensor", "1"="Event"))+
scale_x_discrete(labels=c("0"="Laki-Laki", "1"="Perempuan"))+
guides(fill=guide_legend(title="Jenis Kelamin"))+
theme_minimal() +
theme(panel.grid=element_blank(),legend.position="top",
legend.justification="left",
legend.title=element_text(size=11, face="bold"),
axis.text.x=element_blank(),axis.title.y=element_blank(),
plot.title=element_text(face="bold", size=12),
legend.box.margin=margin(b=-5), legend.margin=margin(l=10))+coord_flip()
# Stacked Bar Chart : Penyakit Komplikasi
ggplot(data, aes(x=komp, fill=event))+
geom_bar(stat="count", position="stack", width=0.7)+
geom_text(aes(label = after_stat(count)), stat="count", vjust=-0.5, color="grey17", size=4,
fontface="bold", position=position_stack(vjust = 0.5))+
labs(fill= "Status", y="Jumlah Pasien")+
scale_fill_manual(values=c("0"="salmon1", "1"="turquoise3"),
labels=c("0"="Sensor", "1"="Event"))+
scale_x_discrete(labels=c("0"="Tidak Ada", "1"="Ada"))+
guides(fill=guide_legend(title="Penyakit Komplikasi"))+
theme_minimal() +
theme(panel.grid=element_blank(),legend.position="top",
legend.justification="left",
legend.title=element_text(size=11, face="bold"),
axis.text.x=element_blank(),axis.title.y=element_blank(),
plot.title=element_text(face="bold", size=12),
legend.box.margin=margin(b=-5), legend.margin=margin(l=10))+coord_flip()
# Stacked Bar Chart : Riwayat Diabetes Melitus
ggplot(data, aes(x=dm, fill=event))+
geom_bar(stat="count", position="stack", width=0.7)+
geom_text(aes(label = after_stat(count)), stat="count", vjust=-0.5, color="grey17", size=4,
fontface="bold", position=position_stack(vjust = 0.5))+
labs(fill= "Status", y="Jumlah Pasien")+
scale_fill_manual(values=c("0"="salmon1", "1"="turquoise3"),
labels=c("0"="Sensor", "1"="Event"))+
scale_x_discrete(labels=c("0"="Tidak Ada", "1"="Ada"))+
guides(fill=guide_legend(title="Riwayat Diabetes Melitus"))+
theme_minimal() +
theme(panel.grid=element_blank(),legend.position="top",
legend.justification="left",
legend.title=element_text(size=11, face="bold"),
axis.text.x=element_blank(),axis.title.y=element_blank(),
plot.title=element_text(face="bold", size=12),
legend.box.margin=margin(b=-5), legend.margin=margin(l=10))+coord_flip()
# Stacked Bar Chart : Kepatuhan terhadap Pengobatan
ggplot(data, aes(x=obat, fill=event))+
geom_bar(stat="count", position="stack", width=0.7)+
geom_text(aes(label = after_stat(count)), stat="count", vjust=-0.5, color="grey17", size=4,
fontface="bold", position=position_stack(vjust = 0.5))+
labs(fill= "Status", y="Jumlah Pasien")+
scale_fill_manual(values=c("0"="salmon1", "1"="turquoise3"),
labels=c("0"="Sensor", "1"="Event"))+
scale_x_discrete(labels=c("0"="Tidak Patuh", "1"="Patuh"))+
guides(fill=guide_legend(title="Kepatuhan terhadap Pengobatan"))+
theme_minimal() +
theme(panel.grid=element_blank(),legend.position="top",
legend.justification="left",
legend.title=element_text(size=11, face="bold"),
axis.text.x=element_blank(),axis.title.y=element_blank(),
plot.title=element_text(face="bold", size=12),
legend.box.margin=margin(b=-5), legend.margin=margin(l=10))+coord_flip()
# Bar Chart Event
data_event <- subset(data, event == 1)
ggplot(data_event, aes(x = factor(`event`))) +
geom_bar(aes(y = after_stat(count)), fill = "skyblue", color = "black") +
geom_text(aes(label = after_stat(count), y = after_stat(count)),
stat = "count", vjust = -0.5) +
labs(x = "Pasien dengan Komplikasi Hipertensi", y = "Jumlah") + theme_minimal()
# Crosstab
tabelcross <- function(var, event) {
tab <- table(var, event)
prop <- prop.table(tab, margin = 1) * 100
hasil <- matrix(paste0(sprintf("%.2f", prop), "% (", tab, ")"),
nrow = nrow(tab), dimnames = dimnames(tab))
total_row <- rowSums(tab)
total_prop <- total_row / sum(tab) * 100
hasil <- cbind(hasil, Total = paste0(sprintf("%.2f", total_prop), "% (", total_row, ")"))
total_col <- colSums(tab)
col_prop <- total_col / sum(tab) * 100
total_total <- paste0("100.00% (", sum(tab), ")")
hasil <- rbind(hasil, Total = c(paste0(sprintf("%.2f", col_prop), "% (", total_col, ")"), total_total))
return(hasil)
}
tabelcross_jk <- tabelcross(data$jk, data$event);tabelcross_jk
## 0 1 Total
## 0 "79.75% (63)" "20.25% (16)" "37.26% (79)"
## 1 "72.93% (97)" "27.07% (36)" "62.74% (133)"
## Total "75.47% (160)" "24.53% (52)" "100.00% (212)"
tabelcross_komp <- tabelcross(data$komp, data$event);tabelcross_komp
## 0 1 Total
## 0 "90.12% (73)" "9.88% (8)" "38.21% (81)"
## 1 "66.41% (87)" "33.59% (44)" "61.79% (131)"
## Total "75.47% (160)" "24.53% (52)" "100.00% (212)"
tabelcross_dm <- tabelcross(data$dm, data$event);tabelcross_dm
## 0 1 Total
## 0 "85.03% (125)" "14.97% (22)" "69.34% (147)"
## 1 "53.85% (35)" "46.15% (30)" "30.66% (65)"
## Total "75.47% (160)" "24.53% (52)" "100.00% (212)"
tabelcross_obat <- tabelcross(data$obat, data$event);tabelcross_obat
## 0 1 Total
## 0 "73.68% (98)" "26.32% (35)" "62.74% (133)"
## 1 "78.48% (62)" "21.52% (17)" "37.26% (79)"
## Total "75.47% (160)" "24.53% (52)" "100.00% (212)"
################# KAPLAN MEIER AND LOG-RANK TEST #################
Y = Surv(data$waktu, data$event == 1)
kmfit <- survfit(Y~1)
km = ggsurvplot(kmfit, data=data, xlab="Waktu hingga pasien mengalami komplikasi HHD (hari)",
ylab=expression(bold(hat(S)(t))), legend="none", surv.median.line = "hv",
palette="turquoise4",break.x.by=11)
## Warning in .add_surv_median(p, fit, type = surv.median.line, fun = fun, :
## Median survival not reached.
ggpar(km, font.main=c(12,"bold"), font.x=c(10,"bold"), font.y=c(10,"bold"),
font.tickslab=c(10))
#Jenis Kelamin
kmfit1 <- survfit(Y~data$jk)
km1 = ggsurvplot(kmfit1, data=data, pval=TRUE, pval.method=TRUE, pval.size=3.5,
pval.coord=c(70,0.47), pval.method.coord=c(70,0.53),
legend.title="Jenis Kelamin",legend.labs=c("Laki-Laki","Perempuan"),
palette=c("turquoise3","tomato2"),break.x.by=11,
xlab="Waktu hingga pasien mengalami komplikasi HHD (hari)", ylab=expression(bold(hat(S)(t))))
ggpar(km1, font.main=c(12,"bold"), font.x=c(10,"bold"), font.y=c(10,"bold"),
font.tickslab=c(10), font.legend=c(9,"bold"), legend=c(0.8,0.2))
#Penyakit Komplikasi
kmfit2 <- survfit(Y~data$komp)
km2 = ggsurvplot(kmfit2, data=data, pval=TRUE, pval.method=TRUE, pval.size=3.5,
pval.coord=c(70,0.47), pval.method.coord=c(70,0.53),
legend.title="Penyakit Komplikasi",legend.labs=c("Tidak Ada","Ada"),
palette=c("turquoise3","tomato2"),break.x.by=11,
xlab="Waktu hingga pasien mengalami komplikasi HHD (hari)", ylab=expression(bold(hat(S)(t))))
ggpar(km2, font.main=c(12,"bold"), font.x=c(10,"bold"), font.y=c(10,"bold"),
font.tickslab=c(10), font.legend=c(9,"bold"), legend=c(0.8,0.2))
#Riwayat Diabetes Melitus
kmfit3 <- survfit(Y~data$dm)
km3 = ggsurvplot(kmfit3, data=data, pval=TRUE, pval.method=TRUE, pval.size=3.5,
pval.coord=c(70,0.47), pval.method.coord=c(70,0.53),
legend.title="Riwayat Diabetes Melitus",legend.labs=c("Tidak Ada","Ada"),
palette=c("turquoise3","tomato2"),break.x.by=11,
xlab="Waktu hingga pasien mengalami komplikasi HHD (hari)", ylab=expression(bold(hat(S)(t))))
ggpar(km3, font.main=c(12,"bold"), font.x=c(10,"bold"), font.y=c(10,"bold"),
font.tickslab=c(10), font.legend=c(9,"bold"), legend=c(0.8,0.2))
#Kepatuhan terhadap Pengobatan
kmfit4 <- survfit(Y~data$obat)
km4 = ggsurvplot(kmfit4, data=data, pval=TRUE, pval.method=TRUE, pval.size=3.5,
pval.coord=c(70,0.47), pval.method.coord=c(70,0.53),
legend.title="Kepatuhan terhadap Pengobatan",legend.labs=c("Tidak Ada","Ada"),
palette=c("turquoise3","tomato2"),break.x.by=11,
xlab="Waktu hingga pasien mengalami komplikasi HHD (hari)", ylab=expression(bold(hat(S)(t))))
ggpar(km4, font.main=c(12,"bold"), font.x=c(10,"bold"), font.y=c(10,"bold"),
font.tickslab=c(10), font.legend=c(9,"bold"), legend=c(0.8,0.2))
## LOG RANK TEST
LR1 <- survdiff(Y~data$jk); LR1; LR1$pvalue; LR1$chisq
## Call:
## survdiff(formula = Y ~ data$jk)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## data$jk=0 79 16 19.6 0.663 1.07
## data$jk=1 133 36 32.4 0.401 1.07
##
## Chisq= 1.1 on 1 degrees of freedom, p= 0.3
## [1] 0.301493
## [1] 1.067578
LR2 <- survdiff(Y~data$komp); LR2; LR2$pvalue; LR2$chisq
## Call:
## survdiff(formula = Y ~ data$komp)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## data$komp=0 81 8 18.8 6.21 9.77
## data$komp=1 131 44 33.2 3.52 9.77
##
## Chisq= 9.8 on 1 degrees of freedom, p= 0.002
## [1] 0.001770674
## [1] 9.773268
LR3 <- survdiff(Y~data$dm); LR3; LR3$pvalue; LR3$chisq
## Call:
## survdiff(formula = Y ~ data$dm)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## data$dm=0 147 22 36.2 5.56 18.3
## data$dm=1 65 30 15.8 12.71 18.3
##
## Chisq= 18.3 on 1 degrees of freedom, p= 2e-05
## [1] 1.858429e-05
## [1] 18.32914
LR4 <- survdiff(Y~data$obat); LR4; LR4$pvalue; LR4$chisq
## Call:
## survdiff(formula = Y ~ data$obat)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## data$obat=0 133 35 31.8 0.330 0.853
## data$obat=1 79 17 20.2 0.518 0.853
##
## Chisq= 0.9 on 1 degrees of freedom, p= 0.4
## [1] 0.3555966
## [1] 0.8533846
################# LN(LN(S(t))) PLOT #################
data$event <- as.numeric(as.character(data$event))
#JENIS KELAMIN
data$jk <- factor(data$jk, levels = c("0", "1"), labels = c("Laki", "Perempuan"))
surv_obj <- Surv(time = data$waktu, event = data$event)
fit <- survfit(surv_obj ~ jk, data = data)
ggsurvplot(
fit,
data = data,
fun = "cloglog",
xscale = "d_m",
conf.int = FALSE,
legend.title = "Jenis Kelamin",
legend.labs = c("Laki", "Perempuan"),
palette = c("blue", "red"),
xlab = "Time in Months",
ylab = "ln(-ln(S(t)))",
ggtheme = theme_minimal())
#PENYAKIT KOMPLIKASI
data$komp <- factor(data$komp, levels = c("0", "1"), labels = c("Tidak", "Ya"))
surv_obj <- Surv(time = data$waktu, event = data$event)
fit <- survfit(surv_obj ~ komp, data = data)
ggsurvplot(
fit,
data = data,
fun = "cloglog",
xscale = "d_m",
conf.int = FALSE,
legend.title = "Komplikasi Penyakit",
legend.labs = c("Tidak", "Ya"),
palette = c("blue", "red"),
xlab = "Time in Months",
ylab = "ln(-ln(S(t)))",
ggtheme = theme_minimal())
#RIWAYAT DIABETES MELITUS
data$dm <- factor(data$dm, levels = c("0", "1"), labels = c("Tidak", "Ya"))
surv_obj <- Surv(time = data$waktu, event = data$event)
fit <- survfit(surv_obj ~ dm, data = data)
ggsurvplot(
fit,
data = data,
fun = "cloglog",
xscale = "d_m",
conf.int = FALSE,
legend.title = "Diabetes Melitus",
legend.labs = c("Tidak", "Ya"),
palette = c("blue", "red"),
xlab = "Time in Months",
ylab = "ln(-ln(S(t)))",
ggtheme = theme_minimal())
#KEPATUHAN TERHADAP PENGOBATAN
data$obat <- factor(data$obat, levels = c("0", "1"), labels = c("Tidak", "Ya"))
surv_obj <- Surv(time = data$waktu, event = data$event)
fit <- survfit(surv_obj ~ obat, data = data)
ggsurvplot(
fit,
data = data,
fun = "cloglog",
xscale = "d_m",
conf.int = FALSE,
legend.title = "Kepatuhan terhadap Pengobatan",
legend.labs = c("Tidak", "Ya"),
palette = c("blue", "red"),
xlab = "Time in Months",
ylab = "ln(-ln(S(t)))",
ggtheme = theme_minimal())
# Kembalikan label menjadi “0” dan “1” lagi
data$obat <- factor(data$obat,
levels = c("Tidak", "Ya"),
labels = c("0", "1"))
data$dm <- factor(data$dm,
levels = c("Tidak", "Ya"),
labels = c("0", "1"))
data$komp <- factor(data$komp,
levels = c("Tidak", "Ya"),
labels = c("0", "1"))
data$jk <- factor(data$jk,
levels = c("Laki", "Perempuan"),
labels = c("0", "1"))
################# REGRESION MODEL WITH ALL DISTRIBUTION #################
reg_weibull <- survreg(Surv(waktu, as.numeric(event == 1)) ~ jk+usia+komp+dm+obat+sist+dias+imt,
data = data, dist = "weibull")
reg_exponential <- survreg(Surv(waktu, as.numeric(event == 1)) ~ jk+usia+komp+dm+obat+sist+dias+imt,
data = data, dist = "exponential")
reg_lognormal <- survreg(Surv(waktu, as.numeric(event == 1)) ~ jk+usia+komp+dm+obat+sist+dias+imt,
data = data, dist = "lognormal")
reg_loglogistic <- survreg(Surv(waktu, as.numeric(event == 1)) ~ jk+usia+komp+dm+obat+sist+dias+imt,
data = data, dist = "loglogistic")
# Tampilkan hasil ringkasan dari setiap model
summary(reg_weibull)
##
## Call:
## survreg(formula = Surv(waktu, as.numeric(event == 1)) ~ jk +
## usia + komp + dm + obat + sist + dias + imt, data = data,
## dist = "weibull")
## Value Std. Error z p
## (Intercept) 6.52682 2.10205 3.10 0.0019
## jk1 -0.23563 0.39016 -0.60 0.5459
## usia -0.00704 0.01517 -0.46 0.6427
## komp1 -1.17279 0.56950 -2.06 0.0395
## dm1 -1.32933 0.45031 -2.95 0.0032
## obat1 1.21994 0.40950 2.98 0.0029
## sist -0.02770 0.01595 -1.74 0.0824
## dias 0.03128 0.01219 2.57 0.0103
## imt -0.01766 0.03122 -0.57 0.5717
## Log(scale) 0.22439 0.12389 1.81 0.0701
##
## Scale= 1.25
##
## Weibull distribution
## Loglik(model)= -369.1 Loglik(intercept only)= -387.6
## Chisq= 37.08 on 8 degrees of freedom, p= 1.1e-05
## Number of Newton-Raphson Iterations: 8
## n= 212
summary(reg_exponential)
##
## Call:
## survreg(formula = Surv(waktu, as.numeric(event == 1)) ~ jk +
## usia + komp + dm + obat + sist + dias + imt, data = data,
## dist = "exponential")
## Value Std. Error z p
## (Intercept) 6.21031 1.68119 3.69 0.00022
## jk1 -0.19094 0.31117 -0.61 0.53946
## usia -0.00615 0.01216 -0.51 0.61291
## komp1 -0.93888 0.44146 -2.13 0.03344
## dm1 -1.08379 0.34110 -3.18 0.00149
## obat1 1.04815 0.31874 3.29 0.00101
## sist -0.02383 0.01255 -1.90 0.05766
## dias 0.02659 0.00965 2.76 0.00586
## imt -0.01452 0.02498 -0.58 0.56098
##
## Scale fixed at 1
##
## Exponential distribution
## Loglik(model)= -370.9 Loglik(intercept only)= -390.3
## Chisq= 38.8 on 8 degrees of freedom, p= 5.4e-06
## Number of Newton-Raphson Iterations: 6
## n= 212
summary(reg_lognormal)
##
## Call:
## survreg(formula = Surv(waktu, as.numeric(event == 1)) ~ jk +
## usia + komp + dm + obat + sist + dias + imt, data = data,
## dist = "lognormal")
## Value Std. Error z p
## (Intercept) 5.09737 2.24143 2.27 0.02296
## jk1 -0.08212 0.42194 -0.19 0.84569
## usia 0.00158 0.01618 0.10 0.92247
## komp1 -0.95976 0.50753 -1.89 0.05862
## dm1 -1.58977 0.46613 -3.41 0.00065
## obat1 1.01299 0.44797 2.26 0.02374
## sist -0.02030 0.01841 -1.10 0.27033
## dias 0.03387 0.01361 2.49 0.01282
## imt -0.03496 0.03577 -0.98 0.32844
## Log(scale) 0.72386 0.11118 6.51 7.5e-11
##
## Scale= 2.06
##
## Log Normal distribution
## Loglik(model)= -369.7 Loglik(intercept only)= -385.4
## Chisq= 31.43 on 8 degrees of freedom, p= 0.00012
## Number of Newton-Raphson Iterations: 4
## n= 212
summary(reg_loglogistic)
##
## Call:
## survreg(formula = Surv(waktu, as.numeric(event == 1)) ~ jk +
## usia + komp + dm + obat + sist + dias + imt, data = data,
## dist = "loglogistic")
## Value Std. Error z p
## (Intercept) 6.06493 2.21606 2.74 0.0062
## jk1 -0.17598 0.40369 -0.44 0.6629
## usia -0.00619 0.01609 -0.38 0.7004
## komp1 -1.05225 0.53849 -1.95 0.0507
## dm1 -1.45227 0.45804 -3.17 0.0015
## obat1 1.17002 0.43014 2.72 0.0065
## sist -0.02513 0.01707 -1.47 0.1411
## dias 0.03016 0.01243 2.43 0.0152
## imt -0.02287 0.03367 -0.68 0.4969
## Log(scale) 0.08526 0.12172 0.70 0.4836
##
## Scale= 1.09
##
## Log logistic distribution
## Loglik(model)= -369.2 Loglik(intercept only)= -386.6
## Chisq= 34.93 on 8 degrees of freedom, p= 2.8e-05
## Number of Newton-Raphson Iterations: 5
## n= 212
# Bandingkan model menggunakan AIC
AIC(reg_weibull, reg_exponential, reg_lognormal, reg_loglogistic)
## df AIC
## reg_weibull 10 758.1122
## reg_exponential 9 759.7169
## reg_lognormal 10 759.3433
## reg_loglogistic 10 758.3666
################# REGRESSION MODEL WITH WEIBULL DISTRIBUTION #################
reg_weibull <- survreg(Surv(waktu, as.numeric(event == 1)) ~ jk+usia+komp+dm+obat+sist+dias+imt,
data = data, dist = "weibull")
reg_weibull
## Call:
## survreg(formula = Surv(waktu, as.numeric(event == 1)) ~ jk +
## usia + komp + dm + obat + sist + dias + imt, data = data,
## dist = "weibull")
##
## Coefficients:
## (Intercept) jk1 usia komp1 dm1 obat1
## 6.526824567 -0.235628568 -0.007038683 -1.172786479 -1.329328455 1.219944889
## sist dias imt
## -0.027701749 0.031283254 -0.017656260
##
## Scale= 1.251561
##
## Loglik(model)= -369.1 Loglik(intercept only)= -387.6
## Chisq= 37.08 on 8 degrees of freedom, p= 1.11e-05
## n= 212
summary(reg_weibull)
##
## Call:
## survreg(formula = Surv(waktu, as.numeric(event == 1)) ~ jk +
## usia + komp + dm + obat + sist + dias + imt, data = data,
## dist = "weibull")
## Value Std. Error z p
## (Intercept) 6.52682 2.10205 3.10 0.0019
## jk1 -0.23563 0.39016 -0.60 0.5459
## usia -0.00704 0.01517 -0.46 0.6427
## komp1 -1.17279 0.56950 -2.06 0.0395
## dm1 -1.32933 0.45031 -2.95 0.0032
## obat1 1.21994 0.40950 2.98 0.0029
## sist -0.02770 0.01595 -1.74 0.0824
## dias 0.03128 0.01219 2.57 0.0103
## imt -0.01766 0.03122 -0.57 0.5717
## Log(scale) 0.22439 0.12389 1.81 0.0701
##
## Scale= 1.25
##
## Weibull distribution
## Loglik(model)= -369.1 Loglik(intercept only)= -387.6
## Chisq= 37.08 on 8 degrees of freedom, p= 1.1e-05
## Number of Newton-Raphson Iterations: 8
## n= 212
logLik(reg_weibull)
## 'log Lik.' -369.0561 (df=10)
#backward
reg_weibull_back <- stepAIC(reg_weibull, direction = "backward")
## Start: AIC=758.11
## Surv(waktu, as.numeric(event == 1)) ~ jk + usia + komp + dm +
## obat + sist + dias + imt
##
## Df AIC
## - usia 1 756.33
## - imt 1 756.42
## - jk 1 756.49
## <none> 758.11
## - sist 1 759.23
## - komp 1 761.02
## - dias 1 764.08
## - obat 1 765.99
## - dm 1 766.54
##
## Step: AIC=756.33
## Surv(waktu, as.numeric(event == 1)) ~ jk + komp + dm + obat +
## sist + dias + imt
##
## Df AIC
## - imt 1 754.56
## - jk 1 754.75
## <none> 756.33
## - sist 1 757.29
## - komp 1 759.64
## - dias 1 762.08
## - obat 1 763.99
## - dm 1 765.05
##
## Step: AIC=754.56
## Surv(waktu, as.numeric(event == 1)) ~ jk + komp + dm + obat +
## sist + dias
##
## Df AIC
## - jk 1 753.09
## <none> 754.56
## - sist 1 755.74
## - komp 1 757.86
## - dias 1 760.22
## - obat 1 762.30
## - dm 1 763.22
##
## Step: AIC=753.09
## Surv(waktu, as.numeric(event == 1)) ~ komp + dm + obat + sist +
## dias
##
## Df AIC
## <none> 753.09
## - sist 1 754.41
## - komp 1 756.16
## - dias 1 758.89
## - obat 1 761.50
## - dm 1 762.80
summary(reg_weibull_back)
##
## Call:
## survreg(formula = Surv(waktu, as.numeric(event == 1)) ~ komp +
## dm + obat + sist + dias, data = data, dist = "weibull")
## Value Std. Error z p
## (Intercept) 5.4825 1.5709 3.49 0.00048
## komp1 -1.1771 0.5644 -2.09 0.03701
## dm1 -1.3951 0.4478 -3.12 0.00184
## obat1 1.2283 0.4028 3.05 0.00230
## sist -0.0261 0.0144 -1.81 0.07018
## dias 0.0308 0.0120 2.56 0.01058
## Log(scale) 0.2267 0.1239 1.83 0.06739
##
## Scale= 1.25
##
## Weibull distribution
## Loglik(model)= -369.5 Loglik(intercept only)= -387.6
## Chisq= 36.1 on 5 degrees of freedom, p= 9.1e-07
## Number of Newton-Raphson Iterations: 8
## n= 212
logLik(reg_weibull_back)
## 'log Lik.' -369.5474 (df=7)
# parameter scale and shape
scale <- reg_weibull_back$scale;scale
## [1] 1.254412
shape <- 1 / scale;shape
## [1] 0.7971865
################# REGRESSION MODEL WITH WEIBULL DISTRIBUTION WITH FRAILTY MODEL GAMMA #################
data$i <- as.numeric(as.character(data$i))
data$event <- as.numeric(as.character(data$event))
data$usia <- as.numeric(as.character(data$usia))
data$sist <- as.numeric(as.character(data$sist))
data$dias <- as.numeric(as.character(data$dias))
data$imt <- as.numeric(as.character(data$imt))
data$jk <- as.numeric(as.character(data$jk))
data$komp <- as.numeric(as.character(data$komp))
data$dm <- as.numeric(as.character(data$dm))
data$obat <- as.numeric(as.character(data$obat))
str(data)
## 'data.frame': 212 obs. of 11 variables:
## $ i : num 1 2 3 4 5 6 7 8 9 10 ...
## $ waktu: num 303 353 37 327 170 348 53 172 101 156 ...
## $ event: num 0 0 1 0 0 0 0 0 0 0 ...
## $ usia : num 66 74 77 54 85 85 89 50 66 58 ...
## $ sist : num 89 90 70 89 72 91 91 77 65 87 ...
## $ dias : num 161 158 140 171 142 136 148 130 101 115 ...
## $ jk : num 1 1 1 1 0 0 1 0 0 0 ...
## $ komp : num 1 1 0 0 1 1 1 1 0 1 ...
## $ dm : num 1 0 0 0 1 0 0 0 0 0 ...
## $ imt : num 26.6 24.3 31.2 23.4 27.1 ...
## $ obat : num 1 1 0 0 0 1 1 1 1 1 ...
# WEIBULL - GAMMA FRAILTY
model_frailty_weibull <- parfm(Surv(waktu, event) ~ jk+usia+komp+dm+obat+sist+dias+imt ,
cluster = "i",
data = data,
dist = "weibull", # baseline hazard
frailty = "gamma") # frailty
model_frailty_weibull
##
## Frailty distribution: gamma
## Baseline hazard distribution: Weibull
## Loglikelihood: -369.017
##
## ESTIMATE SE p-val
## theta 0.276 0.964
## rho 0.834 0.146
## lambda 0.005 0.002
## jk 0.183 0.329 0.577
## usia 0.006 0.011 0.597
## komp 0.948 0.448 0.034 *
## dm 1.147 0.446 0.01 *
## obat -1.012 0.355 0.004 **
## sist 0.023 0.013 0.071 .
## dias -0.026 0.009 0.007 **
## imt 0.016 0.025 0.513
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
as.data.frame(model_frailty_weibull)
## ESTIMATE SE p-val
## theta 0.275575421 0.963896230 NA
## rho 0.833622673 0.146067636 NA
## lambda 0.004662349 0.001978549 NA
## jk 0.183335443 0.328726991 0.577040116
## usia 0.005810627 0.010987824 0.596927428
## komp 0.948434500 0.447861762 0.034200928
## dm 1.146956157 0.446145348 0.010146016
## obat -1.011836204 0.354712939 0.004337006
## sist 0.022730660 0.012577779 0.070729905
## dias -0.025730138 0.009490880 0.006707289
## imt 0.016108552 0.024641206 0.513289511
length(unique(data$i))
## [1] 212
#summary
summary(model_frailty_weibull)
## ESTIMATE SE p-val
## Min. :-1.011836 Min. :0.001978 Min. :0.004337
## 1st Qu.: 0.005236 1st Qu.:0.011783 1st Qu.:0.009286
## Median : 0.022731 Median :0.146068 Median :0.052465
## Mean : 0.218152 Mean :0.249735 Mean :0.226672
## 3rd Qu.: 0.554599 3rd Qu.:0.400429 3rd Qu.:0.529227
## Max. : 1.146956 Max. :0.963896 Max. :0.596927
## NA's :3
#check data structure
#str(model_frailty_weibull)
#ekstract frailty parameter
model_frailty_weibull["theta", ]
## ESTIMATE SE p-val
## 0.2755754 0.9638962 NA
#show attributes
#attributes(model_frailty_weibull)
#confidence interval
coefs <- model_frailty_weibull[, "ESTIMATE"];coefs
## theta rho lambda jk usia komp
## 0.275575421 0.833622673 0.004662349 0.183335443 0.005810627 0.948434500
## dm obat sist dias imt
## 1.146956157 -1.011836204 0.022730660 -0.025730138 0.016108552
se <- model_frailty_weibull[, "SE"]
lower_ci <- coefs - 1.96 * se
upper_ci <- coefs + 1.96 * se
conf_intervals <- data.frame(Estimate = coefs, Lower_CI = lower_ci, Upper_CI = upper_ci)
print(conf_intervals)
## Estimate Lower_CI Upper_CI
## theta 0.275575421 -1.6136611905 2.164812032
## rho 0.833622673 0.5473301052 1.119915240
## lambda 0.004662349 0.0007843937 0.008540304
## jk 0.183335443 -0.4609694595 0.827640346
## usia 0.005810627 -0.0157255079 0.027346762
## komp 0.948434500 0.0706254476 1.826243553
## dm 1.146956157 0.2725112739 2.021401040
## obat -1.011836204 -1.7070735638 -0.316598843
## sist 0.022730660 -0.0019217882 0.047383107
## dias -0.025730138 -0.0443322624 -0.007128014
## imt 0.016108552 -0.0321882108 0.064405316
pred_values <- predict(model_frailty_weibull, type = "risk")
hist(pred_values, main = "Distribution of Frailty", xlab = "Frailty Value", col = "blue")
#comparing AIC weibull and weibull with gamma frailty
AIC(reg_weibull_back)
## [1] 753.0947
AIC(model_frailty_weibull)
## [1] 760.0341
#comparing AIC weibull and weibull with gamma frailty
logLik(reg_weibull_back)
## 'log Lik.' -369.5474 (df=7)
logLik(model_frailty_weibull)
## [1] -369.0171
## attr(,"df")
## [1] 11
## attr(,"nobs")
## [1] 212
## coba untuk backward
#cox_model <- coxph(Surv(waktu, event) ~ jk + usia + komp + dm + obat + sist + dias + imt, data = data)
#step_model <- step(cox_model, direction = "backward")
#STEP BY STEP BACKWARD
model_frailty_weibull2 <- parfm(Surv(waktu, event) ~ jk+komp+dm+obat+sist+dias+imt ,
cluster = "i",
data = data,
dist = "weibull", # baseline hazard
frailty = "gamma") # frailty
AIC(model_frailty_weibull2)
## [1] 758.2377
model_frailty_weibull3 <- parfm(Surv(waktu, event) ~ komp+dm+obat+sist+dias ,
cluster = "i",
data = data,
dist = "weibull", # baseline hazard
frailty = "gamma") # frailty
model_frailty_weibull3
##
## Frailty distribution: gamma
## Baseline hazard distribution: Weibull
## Loglikelihood: -369.478
##
## ESTIMATE SE p-val
## theta 0.395 1.082
## rho 0.846 0.159
## lambda 0.011 0.010
## komp 0.955 0.441 0.031 *
## dm 1.225 0.450 0.006 **
## obat -1.025 0.358 0.004 **
## sist 0.021 0.012 0.08 .
## dias -0.025 0.009 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_frailty_weibull3)
## [1] 754.9565
pred_values <- predict(model_frailty_weibull3, type = "risk")
hist(pred_values, main = "Distribution of Frailty", xlab = "Frailty Value", col = "blue")
lambda <- model_frailty_weibull3[3] # misal parameter skala di posisi 1
rho <- model_frailty_weibull3[2]
# Hitung intercept b0 dari lambda dan rho
b0 <- -log(lambda) / rho
cat("Intercept (b0) =", b0, "\n")
## Intercept (b0) = 5.307939
#comparing AIC weibull and weibull with gamma frailty
AIC(reg_weibull)
## [1] 758.1122
AIC(reg_weibull_back)
## [1] 753.0947
AIC(model_frailty_weibull)
## [1] 760.0341
AIC(model_frailty_weibull3)
## [1] 754.9565