podatki <- read.table("/cloud/project/Poglavje 4/Naloga 2/Aretacija.csv", header=TRUE, sep=";", dec=",", )
str(podatki)
## 'data.frame': 1000 obs. of 10 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Aretacija : int 1 1 0 0 0 1 1 0 0 1 ...
## $ Čas : int 74 239 208 730 730 287 528 730 730 431 ...
## $ Starost : int 22 19 28 36 27 29 24 28 36 25 ...
## $ Starost_kateg: int 2 1 2 3 2 2 2 2 3 2 ...
## $ Poročen : int 0 0 0 1 1 1 0 1 1 0 ...
## $ Izobrazba : int 1 2 1 2 2 2 2 1 1 1 ...
## $ Zaposlen : int 1 0 0 1 0 1 1 1 1 1 ...
## $ Spol : int 0 1 0 0 0 0 0 0 0 0 ...
## $ Nasilje : int 0 0 0 0 0 0 0 0 0 0 ...
Opis spremenljivk:
podatki$StarostF <- factor(podatki$Starost_kateg,
levels = c(1, 2, 3, 4),
labels = c("18-21", "22-29", "30-39", "40+"))
podatki$PoročenF <- factor(podatki$Poročen,
levels = c(0, 1),
labels = c("Ne", "Da"))
podatki$IzobrazbaF <- factor(podatki$Izobrazba,
levels = c(1, 2, 3, 4),
labels = c("Brez SŠ", "SŠ", "Višja", "Univerza"))
podatki$ZaposlenF <- factor(podatki$Zaposlen,
levels = c(0, 1),
labels = c("Ne", "Da"))
podatki$SpolF <- factor(podatki$Spol,
levels = c(0, 1),
labels = c("Moški", "Ženska"))
podatki$NasiljeF <- factor(podatki$Nasilje,
levels = c(0, 1),
labels = c("Ne", "Da"))
str(podatki)
## 'data.frame': 1000 obs. of 16 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Aretacija : int 1 1 0 0 0 1 1 0 0 1 ...
## $ Čas : int 74 239 208 730 730 287 528 730 730 431 ...
## $ Starost : int 22 19 28 36 27 29 24 28 36 25 ...
## $ Starost_kateg: int 2 1 2 3 2 2 2 2 3 2 ...
## $ Poročen : int 0 0 0 1 1 1 0 1 1 0 ...
## $ Izobrazba : int 1 2 1 2 2 2 2 1 1 1 ...
## $ Zaposlen : int 1 0 0 1 0 1 1 1 1 1 ...
## $ Spol : int 0 1 0 0 0 0 0 0 0 0 ...
## $ Nasilje : int 0 0 0 0 0 0 0 0 0 0 ...
## $ StarostF : Factor w/ 4 levels "18-21","22-29",..: 2 1 2 3 2 2 2 2 3 2 ...
## $ PoročenF : Factor w/ 2 levels "Ne","Da": 1 1 1 2 2 2 1 2 2 1 ...
## $ IzobrazbaF : Factor w/ 4 levels "Brez SŠ","SŠ",..: 1 2 1 2 2 2 2 1 1 1 ...
## $ ZaposlenF : Factor w/ 2 levels "Ne","Da": 2 1 1 2 1 2 2 2 2 2 ...
## $ SpolF : Factor w/ 2 levels "Moški","Ženska": 1 2 1 1 1 1 1 1 1 1 ...
## $ NasiljeF : Factor w/ 2 levels "Ne","Da": 1 1 1 1 1 1 1 1 1 1 ...
library(survival)
prež_obj <- Surv(podatki$Čas, podatki$Aretacija)
head(prež_obj, 10)
## [1] 74 239 208+ 730+ 730+ 287 528 730+ 730+ 431
KM_fit <- survfit(prež_obj ~ 1)
summary(KM_fit, times=c(0, 60*(1:12)))
## Call: survfit(formula = prež_obj ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 1000 0 1.000 0.00000 1.000 1.000
## 60 944 56 0.944 0.00727 0.930 0.958
## 120 872 72 0.872 0.01056 0.852 0.893
## 180 798 70 0.802 0.01261 0.777 0.827
## 240 719 75 0.726 0.01413 0.699 0.754
## 300 656 61 0.665 0.01497 0.636 0.695
## 360 600 55 0.609 0.01549 0.579 0.640
## 420 561 37 0.571 0.01572 0.541 0.603
## 480 509 50 0.520 0.01588 0.490 0.552
## 540 461 43 0.476 0.01590 0.446 0.508
## 600 420 41 0.434 0.01580 0.404 0.466
## 660 389 30 0.403 0.01565 0.373 0.435
## 720 358 25 0.377 0.01548 0.347 0.408
library(ggplot2)
library(ggfortify)
KM_fit
## Call: survfit(formula = prež_obj ~ 1)
##
## n events median 0.95LCL 0.95UCL
## [1,] 1000 620 501 465 552
autoplot(KM_fit,
ylab="Stopnja preživetja",
xlab="Čas v dnevih",
main = "Funkcija preživetja") +
theme_linedraw()
KM_starost_fit <- survfit(prež_obj ~ StarostF,
data = podatki)
KM_starost_fit
## Call: survfit(formula = prež_obj ~ StarostF, data = podatki)
##
## n events median 0.95LCL 0.95UCL
## StarostF=18-21 256 206 376 336 456
## StarostF=22-29 374 258 426 359 489
## StarostF=30-39 257 119 NA 683 NA
## StarostF=40+ 113 37 NA NA NA
autoplot(KM_starost_fit,
ylab="Stopnja preživetja",
xlab="Čas v dnevih",
main="Funkcija preživetja") +
theme_linedraw() +
theme(legend.title = element_blank())
survdiff(prež_obj ~ podatki$StarostF)
## Call:
## survdiff(formula = prež_obj ~ podatki$StarostF)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## podatki$StarostF=18-21 256 206 137.5 34.1 44.3
## podatki$StarostF=22-29 374 258 208.1 12.0 18.1
## podatki$StarostF=30-39 257 119 184.5 23.2 33.3
## podatki$StarostF=40+ 113 37 89.9 31.1 36.7
##
## Chisq= 102 on 3 degrees of freedom, p= <2e-16
KM_spol_fit <- survfit(prež_obj ~ SpolF,
data = podatki)
KM_spol_fit
## Call: survfit(formula = prež_obj ~ SpolF, data = podatki)
##
## n events median 0.95LCL 0.95UCL
## SpolF=Moški 845 529 481 445 538
## SpolF=Ženska 155 91 616 506 714
autoplot(KM_spol_fit,
ylab="Stopnja preživetja",
xlab="Čas v dnevih",
main="Funkcija preživetja") +
theme_linedraw() +
theme(legend.title = element_blank())
survdiff(prež_obj ~ podatki$SpolF)
## Call:
## survdiff(formula = prež_obj ~ podatki$SpolF)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## podatki$SpolF=Moški 845 529 514 0.462 2.7
## podatki$SpolF=Ženska 155 91 106 2.229 2.7
##
## Chisq= 2.7 on 1 degrees of freedom, p= 0.1
fit1 <- coxph(prež_obj ~ Starost + PoročenF + IzobrazbaF + ZaposlenF + SpolF,
data = podatki)
summary(fit1)
## Call:
## coxph(formula = prež_obj ~ Starost + PoročenF + IzobrazbaF +
## ZaposlenF + SpolF, data = podatki)
##
## n= 1000, number of events= 620
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Starost -0.05642 0.94515 0.00584 -9.661 < 2e-16 ***
## PoročenFDa -0.59269 0.55284 0.08549 -6.933 4.13e-12 ***
## IzobrazbaFSŠ -0.10897 0.89676 0.08957 -1.217 0.223746
## IzobrazbaFVišja 0.18764 1.20640 0.12008 1.563 0.118131
## IzobrazbaFUniverza -0.14070 0.86875 0.27747 -0.507 0.612104
## ZaposlenFDa -0.31118 0.73258 0.08191 -3.799 0.000145 ***
## SpolFŽenska -0.16243 0.85008 0.11393 -1.426 0.153966
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Starost 0.9451 1.0580 0.9344 0.9560
## PoročenFDa 0.5528 1.8088 0.4676 0.6537
## IzobrazbaFSŠ 0.8968 1.1151 0.7524 1.0688
## IzobrazbaFVišja 1.2064 0.8289 0.9534 1.5265
## IzobrazbaFUniverza 0.8688 1.1511 0.5043 1.4965
## ZaposlenFDa 0.7326 1.3650 0.6239 0.8602
## SpolFŽenska 0.8501 1.1764 0.6799 1.0628
##
## Concordance= 0.644 (se = 0.011 )
## Likelihood ratio test= 175.2 on 7 df, p=<2e-16
## Wald test = 159.3 on 7 df, p=<2e-16
## Score (logrank) test = 164.6 on 7 df, p=<2e-16
fit2 <- coxph(prež_obj ~ Starost + PoročenF + IzobrazbaF + ZaposlenF + SpolF + NasiljeF,
data = podatki)
anova(fit1, fit2, test="Chi")
## Analysis of Deviance Table
## Cox model: response is prež_obj
## Model 1: ~ Starost + PoročenF + IzobrazbaF + ZaposlenF + SpolF
## Model 2: ~ Starost + PoročenF + IzobrazbaF + ZaposlenF + SpolF + NasiljeF
## loglik Chisq Df P(>|Chi|)
## 1 -3932.0
## 2 -3929.4 5.2888 1 0.02146 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(fit2)
## Call:
## coxph(formula = prež_obj ~ Starost + PoročenF + IzobrazbaF +
## ZaposlenF + SpolF + NasiljeF, data = podatki)
##
## n= 1000, number of events= 620
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Starost -0.05620 0.94535 0.00582 -9.655 < 2e-16 ***
## PoročenFDa -0.59829 0.54975 0.08558 -6.991 2.72e-12 ***
## IzobrazbaFSŠ -0.10635 0.89911 0.08960 -1.187 0.235233
## IzobrazbaFVišja 0.18949 1.20864 0.12012 1.578 0.114660
## IzobrazbaFUniverza -0.13743 0.87159 0.27755 -0.495 0.620487
## ZaposlenFDa -0.31170 0.73220 0.08195 -3.804 0.000143 ***
## SpolFŽenska -0.16138 0.85097 0.11394 -1.416 0.156679
## NasiljeFDa 0.35697 1.42899 0.14780 2.415 0.015728 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Starost 0.9454 1.0578 0.9346 0.9562
## PoročenFDa 0.5497 1.8190 0.4649 0.6501
## IzobrazbaFSŠ 0.8991 1.1122 0.7543 1.0717
## IzobrazbaFVišja 1.2086 0.8274 0.9551 1.5295
## IzobrazbaFUniverza 0.8716 1.1473 0.5059 1.5016
## ZaposlenFDa 0.7322 1.3657 0.6236 0.8598
## SpolFŽenska 0.8510 1.1751 0.6807 1.0639
## NasiljeFDa 1.4290 0.6998 1.0696 1.9091
##
## Concordance= 0.647 (se = 0.011 )
## Likelihood ratio test= 180.5 on 8 df, p=<2e-16
## Wald test = 165.7 on 8 df, p=<2e-16
## Score (logrank) test = 170.3 on 8 df, p=<2e-16