################# 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(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(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("DataFix2.xlsx")
data
## # A tibble: 212 × 11
## i t event usia dias sist 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
## dias 0.0000000
## sist 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(1234)
imp <- mice(datamiss, m = 5, method = "pmm", maxit = 10, seed = 1234)
##
## 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 26.70940 24.88889 23.93606 30.62925 41.09139
## 5 21.08281 21.41094 26.03749 23.30905 21.08281
## 9 20.44444 25.33333 16.49396 23.30905 25.71166
## 11 30.49353 23.49524 25.06575 23.42209 23.49524
## 14 23.30905 23.30905 23.30905 35.20430 24.60973
## 18 16.49396 20.02884 23.30905 27.77427 27.05515
## 20 26.66667 23.49524 24.60973 23.30905 23.30905
## 29 20.70312 27.34375 29.71429 18.82711 26.70940
## 30 16.49396 34.89440 21.41094 28.72008 31.11111
## 35 24.88889 23.49524 26.66667 45.26935 24.14152
## 38 27.05515 25.39062 33.16327 28.35306 23.30905
## 48 27.26801 25.15315 33.29865 23.30905 16.01562
## 49 20.06920 23.30905 20.06920 24.22145 23.30905
## 50 22.89282 25.39062 19.77238 23.30905 20.44444
## 53 26.70940 23.30905 26.03749 44.29630 20.06920
## 65 23.30905 26.70940 23.42209 34.13111 18.82711
## 66 27.05515 20.76125 27.05515 43.70447 23.30905
## 68 34.89440 31.21748 16.84747 43.70447 23.30905
## 70 23.87511 26.03749 23.30905 24.22145 19.90997
## 82 23.30905 22.49135 23.30905 16.49396 21.08281
## 84 27.26801 29.17488 26.03749 25.39062 29.17488
## 86 16.01562 24.97399 31.11111 23.43750 23.49524
## 89 23.30905 23.82812 27.39226 26.03749 24.16327
## 92 20.44444 31.11111 23.30905 21.28743 22.89282
## 97 43.70447 24.34176 43.70447 29.40227 24.97399
## 100 27.05515 25.39062 16.49396 16.49396 23.42209
## 104 26.70940 33.69494 23.30905 23.87511 21.64412
## 107 23.87511 24.91990 23.30905 30.62925 20.31250
## 109 29.17488 24.14152 20.06920 23.30905 26.70940
## 110 16.49396 22.89282 21.08281 23.30905 27.05515
## 116 20.44444 31.11111 23.30905 16.49396 23.30905
## 118 30.22222 24.97399 23.30905 20.70312 24.60973
## 121 22.89282 21.28743 23.30905 23.30905 23.30905
## 126 24.34176 24.03441 25.91068 29.17488 30.62925
## 130 30.22222 23.30905 24.88889 23.93606 21.08281
## 132 26.03749 22.60026 24.22145 24.88889 23.30905
## 133 26.66667 25.39062 27.47138 33.16327 23.49524
## 134 41.09139 24.55775 23.30905 25.06575 23.30905
## 135 19.77238 23.30905 22.76944 23.42209 27.53123
## 136 25.33333 21.41094 22.54658 22.86237 19.77238
## 140 19.90997 30.49353 16.01562 23.30905 22.30815
## 142 30.22222 27.26801 27.34375 25.20398 23.30905
## 145 17.08744 23.49524 24.22145 24.22145 30.62925
## 147 20.44444 22.86237 27.05515 31.11111 16.49396
## 151 23.30905 24.55775 23.30905 23.30905 23.30905
## 154 12.48699 25.91068 19.90997 23.87511 23.30905
## 158 23.87511 26.56250 23.30905 23.30905 24.34176
## 161 23.30905 24.16327 22.89282 26.03749 23.49524
## 162 23.30905 26.03749 28.51562 34.13111 21.08281
## 163 33.20312 23.30905 24.22145 25.06575 27.34375
## 164 22.86237 33.69494 30.62925 23.30905 25.33333
## 172 21.08281 23.49524 29.17488 23.30905 29.17488
## 175 23.87511 23.42209 27.47138 18.82711 25.39062
## 176 23.30905 23.30905 34.96358 27.39226 24.14152
## 179 20.06920 34.13111 23.93606 33.29865 28.51562
## 180 22.86237 23.30905 41.09139 24.14152 25.71166
## 181 26.70940 25.91068 26.83865 23.30905 34.96358
## 183 23.30905 24.22145 28.51562 41.09139 22.30815
## 184 21.28743 24.80159 18.13202 26.03749 24.55775
## 188 41.09139 24.22145 25.80645 22.49135 18.82711
## 189 23.30905 24.97399 23.30905 43.70447 20.76125
## 190 23.30905 24.80159 17.08744 34.96358 26.83865
## 191 22.86237 23.30905 25.15315 24.55775 23.30905
## 192 26.03749 23.87511 17.08744 23.30905 17.08744
## 193 20.44444 23.30905 27.05515 23.49524 16.84747
## 196 16.49396 12.48699 21.08281 16.49396 29.17488
## 200 12.48699 27.53123 30.62925 43.70447 23.30905
## 202 23.30905 27.34375 26.50212 26.66667 23.43750
## 203 24.22145 25.22409 18.82711 22.76944 27.77427
## 204 23.87511 33.20312 21.28743 24.22145 27.77427
## 205 25.71166 26.56250 26.70940 26.56250 23.49524
## 207 31.11111 33.69494 24.97399 34.96358 23.30905
## 208 31.11111 21.64412 16.01562 24.55775 23.30905
## 209 28.72008 24.14152 23.30905 45.26935 26.66667
# 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 161 89 1 1 1 26.56250 1
## 2 2 353 0 74 158 90 1 1 0 24.34176 1
## 3 3 37 1 77 140 70 1 0 0 31.21748 0
## 4 4 327 0 54 171 89 1 0 0 26.70940 0
## 5 5 170 0 85 142 72 0 1 1 21.08281 0
## 6 6 348 0 85 136 91 0 1 0 27.05515 1
## 7 7 53 0 89 148 91 1 1 0 20.44444 1
## 8 8 172 0 50 130 77 0 1 0 25.39062 1
## 9 9 101 0 66 101 65 0 0 0 20.44444 1
## 10 10 156 0 58 115 87 0 1 0 20.02884 1
## 11 11 215 1 74 147 65 1 1 1 30.49353 0
## 12 12 45 0 58 128 77 1 1 1 27.05515 1
## 13 13 192 0 51 135 80 0 0 0 27.34375 0
## 14 14 30 1 56 179 109 0 0 1 23.30905 0
## 15 15 279 0 68 148 68 1 1 1 24.91990 0
## 16 16 100 0 81 150 80 1 1 0 16.49396 1
## 17 17 67 1 49 130 80 1 1 1 41.09139 0
## 18 18 146 0 78 146 74 1 1 0 16.49396 0
## 19 19 94 1 77 105 59 1 1 1 19.90997 0
## 20 20 181 0 42 156 80 1 0 0 26.66667 0
## 21 21 171 0 61 168 71 1 1 1 24.88889 1
## 22 22 30 1 73 130 80 1 1 1 21.64412 0
## 23 23 223 0 78 125 56 1 1 0 23.30905 1
## 24 24 172 0 46 130 70 1 0 0 26.66667 0
## 25 25 115 0 58 149 90 1 0 0 24.14152 0
## 26 26 91 1 74 130 60 0 1 0 16.84747 1
## 27 27 166 0 59 178 94 1 0 0 25.29938 0
## 28 28 149 0 87 171 76 1 1 0 21.08281 1
## 29 29 299 0 39 153 94 1 0 0 20.70312 0
## 30 30 83 0 71 180 96 1 0 0 16.49396 0
## 31 31 207 1 80 140 90 1 1 1 26.66667 1
## 32 32 146 0 69 144 69 1 1 1 20.31250 1
## 33 33 30 1 70 130 77 1 1 1 24.97399 1
## 34 34 320 0 69 130 90 1 1 1 25.15315 1
## 35 35 69 0 49 170 99 1 0 0 24.88889 0
## 36 36 191 0 64 184 87 0 1 1 25.39062 1
## 37 37 84 1 71 151 80 0 1 1 23.30905 0
## 38 38 143 0 64 160 90 0 1 0 27.05515 0
## 39 39 167 0 60 153 93 1 0 0 26.66667 0
## 40 40 200 0 51 140 60 1 0 1 28.35306 0
## 41 41 35 1 37 153 89 1 0 0 33.29865 0
## 42 42 359 0 71 134 53 0 1 1 34.89440 1
## 43 43 314 0 51 152 83 1 1 0 26.03749 0
## 44 44 135 0 58 145 65 0 1 1 23.87511 1
## 45 45 7 1 44 140 80 1 1 1 22.54658 0
## 46 46 200 0 56 130 89 1 0 0 23.30905 0
## 47 47 32 1 70 150 60 1 1 1 24.03441 0
## 48 48 216 0 68 177 96 0 1 0 27.26801 0
## 49 49 138 0 39 163 81 1 1 0 20.06920 0
## 50 50 72 0 66 130 80 1 1 1 22.89282 1
## 51 51 38 1 75 176 92 0 1 1 33.69494 1
## 52 52 19 1 69 150 90 0 1 1 21.41094 1
## 53 53 13 1 37 160 93 1 0 1 26.70940 0
## 54 54 171 0 53 140 79 1 1 1 24.88889 1
## 55 55 32 1 58 140 90 0 1 0 34.96358 0
## 56 56 278 0 40 179 78 1 1 0 23.30905 1
## 57 57 310 0 77 170 70 1 1 0 12.48699 1
## 58 58 311 0 38 138 74 1 0 0 33.29865 0
## 59 59 83 1 65 140 80 1 0 0 23.82812 1
## 60 60 130 1 49 150 109 1 1 0 23.30905 0
## 61 61 331 0 28 146 77 1 0 0 20.06920 0
## 62 62 38 1 54 128 79 1 1 0 27.34375 0
## 63 63 362 0 76 124 69 0 1 1 22.89282 1
## 64 64 305 0 77 130 78 1 1 1 25.71166 1
## 65 65 317 0 44 180 117 0 0 0 23.30905 0
## 66 66 178 0 64 108 67 1 0 0 27.05515 0
## 67 67 156 0 66 190 91 1 0 0 28.72008 0
## 68 68 101 0 57 149 89 1 0 0 34.89440 0
## 69 69 321 0 37 150 80 0 0 0 24.22145 0
## 70 70 199 0 55 140 90 1 1 1 23.87511 0
## 71 71 32 1 70 120 70 1 1 0 25.20398 0
## 72 72 175 1 70 150 80 1 1 0 30.22222 0
## 73 73 166 0 62 158 103 0 0 0 20.76125 0
## 74 74 282 0 50 141 70 1 1 0 23.30905 1
## 75 75 296 0 70 134 60 1 1 0 30.22222 1
## 76 76 18 1 31 137 90 0 1 0 24.22145 0
## 77 77 167 1 66 140 70 1 1 1 29.17488 1
## 78 78 228 0 35 140 100 1 1 0 23.93606 0
## 79 79 338 0 56 168 90 0 1 1 23.87511 1
## 80 80 314 0 75 163 69 1 1 1 27.53123 1
## 81 81 355 0 80 125 65 1 1 1 23.30905 1
## 82 82 311 0 86 157 67 0 1 0 23.30905 1
## 83 83 340 0 74 145 66 1 1 1 23.30905 0
## 84 84 327 0 54 133 62 1 1 0 27.26801 1
## 85 85 114 0 49 140 101 0 0 0 23.87511 0
## 86 86 115 0 60 149 89 1 0 0 16.01562 0
## 87 87 63 1 63 119 75 0 1 0 22.30815 0
## 88 88 60 1 64 125 85 1 1 1 23.30905 0
## 89 89 45 0 55 157 76 1 0 0 23.30905 0
## 90 90 160 1 46 121 89 1 1 1 34.13111 1
## 91 91 341 0 39 146 98 1 1 0 33.20312 1
## 92 92 131 0 62 130 92 0 1 0 20.44444 1
## 93 93 99 1 80 170 90 0 1 1 26.83865 1
## 94 94 306 0 57 157 91 0 1 0 23.30905 1
## 95 95 18 0 54 138 91 1 1 0 25.80645 1
## 96 96 23 1 68 146 77 0 1 1 23.30905 0
## 97 97 132 0 43 151 109 0 1 0 43.70447 0
## 98 98 37 0 25 127 78 1 0 0 22.76944 0
## 99 99 354 0 64 151 73 0 1 1 23.30905 1
## 100 100 314 0 59 154 94 0 0 0 27.05515 1
## 101 101 66 0 51 152 94 1 1 0 30.49353 1
## 102 102 201 0 57 161 97 1 1 1 23.30905 1
## 103 103 355 0 71 143 91 0 1 0 22.49135 1
## 104 104 190 1 62 136 85 0 1 0 26.70940 0
## 105 105 24 0 64 150 86 1 0 0 23.49524 0
## 106 106 131 1 75 169 71 1 1 0 23.30905 0
## 107 107 139 0 52 150 90 0 1 0 23.87511 0
## 108 108 75 1 52 127 78 1 1 1 23.30905 1
## 109 109 311 0 62 178 104 0 1 1 29.17488 0
## 110 110 18 0 77 165 79 0 0 0 16.49396 1
## 111 111 187 1 47 130 85 1 1 0 20.70312 1
## 112 112 86 0 59 145 84 1 1 0 27.26801 1
## 113 113 160 0 63 141 63 0 1 1 23.30905 1
## 114 114 27 1 31 140 100 0 1 1 29.40227 1
## 115 115 30 1 86 132 57 0 1 1 22.86237 0
## 116 116 151 0 66 157 90 0 0 0 20.44444 1
## 117 117 321 0 52 192 93 1 0 0 33.16327 0
## 118 118 332 0 61 155 102 1 0 0 30.22222 0
## 119 119 59 0 62 140 70 1 1 0 26.70940 1
## 120 120 348 0 43 116 66 0 1 1 18.13202 1
## 121 121 36 1 78 148 74 1 1 0 22.89282 0
## 122 122 258 0 40 158 121 1 0 0 32.45568 0
## 123 123 137 0 52 168 91 1 0 0 26.50212 0
## 124 124 182 1 43 190 113 1 1 1 26.70940 1
## 125 125 142 0 39 170 108 1 0 0 23.30905 0
## 126 126 342 0 55 150 80 0 0 0 24.34176 0
## 127 127 118 0 52 187 87 0 1 0 26.03749 0
## 128 128 208 0 65 190 113 1 0 0 26.70940 1
## 129 129 174 1 63 125 65 1 1 1 23.30905 0
## 130 130 151 0 66 140 90 1 1 0 30.22222 0
## 131 131 170 0 51 149 92 0 0 0 23.30905 0
## 132 132 362 0 47 154 97 0 1 1 26.03749 1
## 133 133 356 0 64 148 81 1 1 0 26.66667 0
## 134 134 200 0 35 130 90 0 1 0 41.09139 0
## 135 135 327 0 70 156 85 1 1 0 19.77238 0
## 136 136 83 0 60 186 97 0 1 1 25.33333 1
## 137 137 156 0 61 133 86 0 0 0 24.97399 0
## 138 138 7 1 48 159 84 1 1 0 45.26935 0
## 139 139 93 1 42 160 111 1 1 1 18.82711 0
## 140 140 44 0 47 147 100 0 0 0 19.90997 0
## 141 141 363 0 63 122 53 1 1 1 25.33333 1
## 142 142 339 0 59 153 92 1 1 0 30.22222 1
## 143 143 31 1 60 140 70 1 0 1 26.43807 1
## 144 144 282 0 69 160 90 0 1 0 31.11111 1
## 145 145 335 0 41 171 93 1 1 0 17.08744 0
## 146 146 11 0 23 140 90 0 1 0 27.77427 1
## 147 147 3 0 90 147 63 1 0 0 20.44444 0
## 148 148 263 0 59 139 58 1 0 0 27.39226 0
## 149 149 361 0 61 180 100 1 0 0 22.60026 0
## 150 150 116 0 62 149 73 1 1 1 23.30905 1
## 151 151 103 0 29 153 96 0 0 0 23.30905 0
## 152 152 319 0 41 145 87 0 1 0 43.70447 0
## 153 153 202 0 37 154 103 0 1 0 23.30905 0
## 154 154 345 0 54 136 83 0 0 0 12.48699 0
## 155 155 293 0 54 236 73 1 1 1 23.30905 1
## 156 156 333 0 43 180 136 0 1 0 24.60973 1
## 157 157 325 0 34 143 88 1 0 0 17.08744 0
## 158 158 171 0 45 147 89 0 1 0 23.87511 1
## 159 159 178 0 70 156 88 0 0 0 25.39062 0
## 160 160 31 1 40 160 99 1 1 0 25.06575 0
## 161 161 293 0 66 171 88 1 0 0 23.30905 0
## 162 162 60 0 69 181 85 1 1 1 23.30905 0
## 163 163 291 0 52 158 81 1 1 1 33.20312 0
## 164 164 292 0 66 145 80 1 1 0 22.86237 0
## 165 165 206 0 57 163 68 0 1 0 24.16327 1
## 166 166 133 1 59 148 60 0 1 1 29.71429 0
## 167 167 153 0 32 136 100 1 0 0 23.42209 1
## 168 168 221 0 24 181 140 1 0 0 28.51562 0
## 169 169 12 0 63 167 93 1 1 0 21.28743 1
## 170 170 209 0 52 160 90 0 1 0 25.91068 0
## 171 171 171 1 52 150 80 1 1 1 22.82996 0
## 172 172 202 0 64 168 90 0 0 0 21.08281 0
## 173 173 202 0 35 156 108 1 0 0 44.29630 0
## 174 174 1 1 52 130 86 0 1 0 26.03749 1
## 175 175 209 0 52 140 90 1 1 1 23.87511 0
## 176 176 234 0 64 156 65 1 0 0 23.30905 0
## 177 177 112 1 60 171 95 1 0 1 24.55775 1
## 178 178 156 0 48 129 68 0 1 0 24.80159 0
## 179 179 126 1 22 143 109 1 1 0 20.06920 0
## 180 180 207 0 67 159 87 0 1 0 22.86237 0
## 181 181 200 0 52 140 80 0 0 0 26.70940 0
## 182 182 101 1 55 145 85 1 1 0 23.43750 0
## 183 183 199 0 43 120 50 0 1 1 23.30905 0
## 184 184 178 0 52 160 90 1 0 0 21.28743 0
## 185 185 171 0 54 167 88 1 0 0 23.49524 0
## 186 186 138 0 49 234 110 1 1 1 35.67182 0
## 187 187 34 1 71 151 83 1 1 0 24.34176 0
## 188 188 4 1 25 134 64 0 0 1 41.09139 1
## 189 189 33 1 55 173 107 1 1 0 23.30905 0
## 190 190 92 1 71 168 92 1 1 1 23.30905 0
## 191 191 142 0 54 156 100 1 0 0 22.86237 0
## 192 192 142 0 36 180 100 0 0 0 26.03749 0
## 193 193 122 0 68 130 90 0 0 0 20.44444 0
## 194 194 96 0 35 165 108 1 0 0 28.00000 0
## 195 195 129 0 27 145 86 1 0 0 27.47138 0
## 196 196 128 0 59 170 107 0 0 0 16.49396 1
## 197 197 30 0 41 142 82 1 0 0 25.22409 1
## 198 198 115 0 50 162 104 1 0 0 22.89282 0
## 199 199 118 0 76 146 82 1 0 0 19.77238 0
## 200 200 115 0 67 170 100 1 1 0 12.48699 0
## 201 201 109 0 81 182 90 0 1 0 16.01562 0
## 202 202 100 0 42 165 90 0 0 0 23.30905 0
## 203 203 103 0 26 199 124 0 1 0 24.22145 0
## 204 204 69 0 33 128 82 1 0 0 23.87511 0
## 205 205 67 0 56 154 92 0 0 0 25.71166 0
## 206 206 3 0 59 139 69 0 0 0 26.56250 0
## 207 207 20 0 60 130 80 0 1 0 31.11111 0
## 208 208 20 0 51 146 93 1 0 0 31.11111 0
## 209 209 10 0 37 155 106 1 0 0 28.72008 0
## 210 210 12 0 47 135 85 1 0 0 35.20430 0
## 211 211 17 0 28 123 83 1 0 0 30.62925 0
## 212 212 4 0 55 145 79 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 22.89 24.19 25.25 27.06 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.90933, p-value = 4.455e-10
# 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 = 15418, p-value = 0.392
## 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.
## 101.0 138.8 149.0 150.5 160.2 236.0
summary(data$dias)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 50.00 76.75 87.00 85.24 92.25 140.00
summary(data$imt)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 12.49 22.89 24.19 25.25 27.06 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 VISUALISATION #################
# Violinplot usia
ggplot(data, aes(x = factor(event), y = usia, fill = factor(event))) +
geom_violin(trim = FALSE, alpha = 0.6) +
geom_boxplot(width = 0.1, outlier.shape = NA, color = "black", fill = "white") + # boxplot tipis di tengah
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"))
## Warning: Removed 46 rows containing missing values or values outside the scale range
## (`geom_violin()`).
# Violinplot sistolik
ggplot(data, aes(x = factor(event), y = sist, fill = factor(event))) +
geom_violin(trim = FALSE, alpha = 0.6) +
geom_boxplot(width = 0.1, outlier.shape = NA, color = "black", fill = "white") +
ylim(0, 300) +
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"))
# Violinplot diastolik
ggplot(data, aes(x = factor(event), y = dias, fill = factor(event))) +
geom_violin(trim = FALSE, alpha = 0.6) +
geom_boxplot(width = 0.1, outlier.shape = NA, color = "black", fill = "white") +
ylim(0, 200) +
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"))
# Violinplot untuk IMT
ggplot(data, aes(x = factor(event), y = imt, fill = factor(event))) +
geom_violin(trim = FALSE, alpha = 0.6) +
geom_boxplot(width = 0.1, outlier.shape = NA, color = "black", fill = "white") +
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"))
################# 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
## 7.163315685 -0.292981679 -0.009262852 -1.136502801 -1.325445394 1.214833345
## sist dias imt
## 0.031356180 -0.027803799 -0.036966658
##
## Scale= 1.255649
##
## Loglik(model)= -368.6 Loglik(intercept only)= -387.6
## Chisq= 38.06 on 8 degrees of freedom, p= 7.32e-06
## 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) 7.16332 2.19896 3.26 0.0011
## jk1 -0.29298 0.39004 -0.75 0.4526
## usia -0.00926 0.01563 -0.59 0.5535
## komp1 -1.13650 0.57373 -1.98 0.0476
## dm1 -1.32545 0.45283 -2.93 0.0034
## obat1 1.21483 0.41089 2.96 0.0031
## sist 0.03136 0.01218 2.57 0.0100
## dias -0.02780 0.01597 -1.74 0.0818
## imt -0.03697 0.03223 -1.15 0.2514
## Log(scale) 0.22765 0.12393 1.84 0.0662
##
## Scale= 1.26
##
## Weibull distribution
## Loglik(model)= -368.6 Loglik(intercept only)= -387.6
## Chisq= 38.06 on 8 degrees of freedom, p= 7.3e-06
## Number of Newton-Raphson Iterations: 8
## n= 212
logLik(reg_weibull)
## 'log Lik.' -368.5641 (df=10)
#backward
reg_weibull_back <- stepAIC(reg_weibull, direction = "backward")
## Start: AIC=757.13
## Surv(waktu, as.numeric(event == 1)) ~ jk + usia + komp + dm +
## obat + sist + dias + imt
##
## Df AIC
## - usia 1 755.48
## - jk 1 755.71
## - imt 1 756.42
## <none> 757.13
## - dias 1 758.27
## - komp 1 759.63
## - sist 1 763.17
## - obat 1 764.83
## - dm 1 765.38
##
## Step: AIC=755.48
## Surv(waktu, as.numeric(event == 1)) ~ jk + komp + dm + obat +
## sist + dias + imt
##
## Df AIC
## - jk 1 754.12
## - imt 1 754.56
## <none> 755.48
## - dias 1 756.28
## - komp 1 758.53
## - sist 1 761.19
## - obat 1 762.84
## - dm 1 764.09
##
## Step: AIC=754.12
## Surv(waktu, as.numeric(event == 1)) ~ komp + dm + obat + sist +
## dias + imt
##
## Df AIC
## - imt 1 753.09
## <none> 754.12
## - dias 1 755.06
## - komp 1 756.82
## - sist 1 759.97
## - obat 1 762.02
## - dm 1 764.01
##
## Step: AIC=753.09
## Surv(waktu, as.numeric(event == 1)) ~ komp + dm + obat + sist +
## dias
##
## Df AIC
## <none> 753.09
## - dias 1 754.41
## - komp 1 756.16
## - sist 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.0308 0.0120 2.56 0.01058
## dias -0.0261 0.0144 -1.81 0.07018
## 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$scale;scale
## [1] 1.255649
shape <- 1 / scale;shape
## [1] 0.796401
################# 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 161 158 140 171 142 136 148 130 101 115 ...
## $ dias : num 89 90 70 89 72 91 91 77 65 87 ...
## $ 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 26.7 21.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: -368.509
##
## ESTIMATE SE p-val
## theta 0.416 1.212
## rho 0.848 0.166
## lambda 0.003 0.000
## jk 0.229 0.334 0.493
## usia 0.007 0.012 0.523
## komp 0.925 0.459 0.044 *
## dm 1.172 0.476 0.014 *
## obat -1.018 0.369 0.006 **
## sist -0.026 0.010 0.008 **
## dias 0.023 0.013 0.08 .
## imt 0.034 0.027 0.209
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#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] 757.317
model_frailty_weibull3 <- parfm(Surv(waktu, event) ~ jk+komp+dm+obat+sist+dias ,
cluster = "i",
data = data,
dist = "weibull", # baseline hazard
frailty = "gamma") # frailty
AIC(model_frailty_weibull3)
## [1] 756.5122
model_frailty_weibull4 <- parfm(Surv(waktu, event) ~ komp+dm+obat+sist+dias ,
cluster = "i",
data = data,
dist = "weibull", # baseline hazard
frailty = "gamma") # frailty
model_frailty_weibull4
##
## 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.025 0.009 0.004 **
## dias 0.021 0.012 0.08 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model_frailty_weibull4)
## [1] 754.9565
#comparing AIC weibull and weibull with gamma frailty
AIC(reg_weibull)
## [1] 757.1282
AIC(reg_weibull_back)
## [1] 753.0947
AIC(model_frailty_weibull)
## [1] 759.0183
AIC(model_frailty_weibull4)
## [1] 754.9565