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