podatki <- read.table("/cloud/project/Poglavje 4/Zgled/Stečaj.csv", header=TRUE, sep=";", dec=",", )

head(podatki)
##   ID Delovanje Stečaj Likvidnost Oprijemljiva     LnTA Panoga Tip
## 1  1       144      0   1.298389     40.42222 15.78642      0   0
## 2  2       144      0   1.204970     44.02762 15.45020      1   0
## 3  3        81      1   1.680213     13.58690 14.65593      0   0
## 4  4       144      0   1.497829     25.67354 16.27653      0   0
## 5  5        12      0   1.119343     15.62781 14.26124      2   0
## 6  6        12      0   1.563760     31.07781 17.83592      1   0
##   Subvencija
## 1          1
## 2          1
## 3          0
## 4          1
## 5          0
## 6          0

Opis spremenljivk:

podatki$StečajF <- factor(podatki$Stečaj, 
                          levels = c(0, 1), 
                          labels = c("Ne", "Da"))

podatki$PanogaF <- factor(podatki$Panoga, 
                          levels = c(0, 1, 2), 
                          labels = c("Storitve", "Proizvod.", "Gradb."))

podatki$TipF <- factor(podatki$Tip, 
                       levels = c(0, 1), 
                       labels = c("d.o.o.", "d.d."))

podatki$SubvencijaF <- factor(podatki$Subvencija, 
                              levels = c(0, 1), 
                              labels = c("NE", "DA"))
summary(podatki[c(-1, -3, -7, -8, -9)]) #Opisna statistika za izbranbe spremenljivke
##    Delovanje       Likvidnost       Oprijemljiva        LnTA      
##  Min.   :  1.0   Min.   : 0.0000   Min.   : 0.00   Min.   :11.76  
##  1st Qu.: 81.0   1st Qu.: 0.9168   1st Qu.:17.54   1st Qu.:14.49  
##  Median :144.0   Median : 1.2837   Median :33.95   Median :15.38  
##  Mean   :114.8   Mean   : 1.7813   Mean   :35.62   Mean   :15.43  
##  3rd Qu.:144.0   3rd Qu.: 2.0322   3rd Qu.:51.69   3rd Qu.:16.29  
##  Max.   :144.0   Max.   :13.1511   Max.   :96.53   Max.   :19.49  
##  StečajF        PanogaF         TipF      SubvencijaF
##  Ne:1784   Storitve :1067   d.o.o.:1979   NE: 855    
##  Da: 570   Proizvod.:1000   d.d.  : 375   DA:1499    
##            Gradb.   : 287                            
##                                                      
##                                                      
## 

Tabele preživetja

#Ustvarimo preživetveni objekt, ki je sestavljen iz spremenljivke, ki zajema časovno komponento, in indikatorja z vrednostma 0 in 1 (v konkretnem primeru je to Stečaj). Indikator ne sme biti spremenjen v faktor.

library(survival) #Aktiviramo knjižnico
prež_obj <- Surv(podatki$Delovanje, podatki$Stečaj)  
head(prež_obj, 10) #Prikaz preživetnega objekta za prvih 10 enot
##  [1] 144+ 144+  81  144+  12+  12+ 144+ 144+ 144+ 144+
library(survival) #Aktiviramo knjižnico
KM_fit <- survfit(prež_obj ~ 1) #Izračun tabele preživetja po KM metodi za vse enote skupaj

summary(KM_fit, times = c(0, 12*(1:12))) #Prikaz tabele preživetja po izbranih časovnih intervalih med 0 in 12 let
## Call: survfit(formula = prež_obj ~ 1)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0   2354       0    1.000 0.00000        1.000        1.000
##    12   2317      37    0.984 0.00256        0.979        0.989
##    24   2241      51    0.962 0.00393        0.955        0.970
##    36   2161      49    0.941 0.00489        0.932        0.951
##    48   2031      96    0.899 0.00630        0.886        0.911
##    60   1934      66    0.869 0.00707        0.855        0.883
##    72   1835      75    0.835 0.00782        0.820        0.850
##    84   1751      63    0.806 0.00835        0.790        0.822
##    96   1692      39    0.788 0.00866        0.771        0.805
##   108   1620      53    0.763 0.00904        0.745        0.781
##   120   1579      19    0.754 0.00916        0.736        0.772
##   132   1548      20    0.744 0.00930        0.726        0.762
##   144   1518       2    0.743 0.00931        0.725        0.762
KM_fit #Opisna statistika 
## Call: survfit(formula = prež_obj ~ 1)
## 
##         n events median 0.95LCL 0.95UCL
## [1,] 2354    570     NA      NA      NA
library(ggplot2)
library(ggfortify)
autoplot(KM_fit,
         ylab = "Stopnja preživetja", 
         xlab = "Čas v mesecih", 
         main = "Funkcija preživetja") +
         theme_linedraw()

KM_panoga_fit <- survfit(prež_obj ~ podatki$PanogaF) #Analiza preživetja glede na panoge

KM_panoga_fit #Izpis rezultatov
## Call: survfit(formula = prež_obj ~ podatki$PanogaF)
## 
##                              n events median 0.95LCL 0.95UCL
## podatki$PanogaF=Storitve  1067    185     NA      NA      NA
## podatki$PanogaF=Proizvod. 1000    247     NA      NA      NA
## podatki$PanogaF=Gradb.     287    138     NA      97      NA
autoplot(KM_panoga_fit, #Grafičen prikaz funkcije preživetja po skupinah
         ylab = "Stopnja preživetja", 
         xlab = "Čas v mesecih", 
         main = "Funkcija preživetja") +
         theme_linedraw() +
         theme(legend.title = element_blank())

survdiff(prež_obj ~ podatki$PanogaF) #Statistična preverba razlik v stopnji preživetja med panogami
## Call:
## survdiff(formula = prež_obj ~ podatki$PanogaF)
## 
##                              N Observed Expected (O-E)^2/E (O-E)^2/V
## podatki$PanogaF=Storitve  1067      185    269.4    26.422    50.275
## podatki$PanogaF=Proizvod. 1000      247    241.0     0.148     0.258
## podatki$PanogaF=Gradb.     287      138     59.6   103.058   115.699
## 
##  Chisq= 130  on 2 degrees of freedom, p= <2e-16
summary(KM_panoga_fit, times = c(48, 96, 144)) #Stopnje preživetja za izbrane časovne enote
## Call: survfit(formula = prež_obj ~ podatki$PanogaF)
## 
##                 podatki$PanogaF=Storitve 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    48    953      63    0.939 0.00743        0.925        0.954
##    96    821      85    0.852 0.01123        0.830        0.875
##   144    740      37    0.813 0.01243        0.789        0.838
## 
##                 podatki$PanogaF=Proizvod. 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    48    854     108    0.890  0.0100        0.870        0.910
##    96    724      96    0.787  0.0132        0.762        0.814
##   144    650      43    0.740  0.0143        0.712        0.768
## 
##                 podatki$PanogaF=Gradb. 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    48    224      62    0.783  0.0244        0.737        0.832
##    96    147      62    0.558  0.0298        0.502        0.619
##   144    128      14    0.503  0.0302        0.448        0.566
KM_tip_fit <- survfit(prež_obj ~ podatki$TipF) #Analiza preživetja glede na tip podjetja

KM_tip_fit
## Call: survfit(formula = prež_obj ~ podatki$TipF)
## 
##                        n events median 0.95LCL 0.95UCL
## podatki$TipF=d.o.o. 1979    452     NA      NA      NA
## podatki$TipF=d.d.    375    118     NA      NA      NA
autoplot(KM_tip_fit, 
         ylab = "Stopnja preživetja", 
         xlab = "Čas v mesecih", 
         main  ="Funkcija preživetja") +
         theme_linedraw() +
         theme(legend.title = element_blank())

survdiff(prež_obj ~ podatki$TipF)
## Call:
## survdiff(formula = prež_obj ~ podatki$TipF)
## 
##                        N Observed Expected (O-E)^2/E (O-E)^2/V
## podatki$TipF=d.o.o. 1979      452    485.5      2.31      15.6
## podatki$TipF=d.d.    375      118     84.5     13.27      15.6
## 
##  Chisq= 15.6  on 1 degrees of freedom, p= 8e-05

Coxova regresija

#Izvedba Coxove regresije. Odvisna spremenljivka mora biti preživetveni objekt, ustvarjen s funkcijo Surv

library(survival) #Aktiviramo knjižnico
fit0 <- coxph(prež_obj ~ 1, #Model brez pojasnjevalnih spremenljivk
             data = podatki)

summary(fit0)
## Call:  coxph(formula = prež_obj ~ 1, data = podatki)
## 
## Null model
##   log likelihood= -4313.223 
##   n= 2354
fit1 <- coxph(prež_obj ~ Likvidnost + Oprijemljiva + LnTA + PanogaF + TipF + SubvencijaF, 
              data = podatki)
library(survival) #Aktiviramo knjižnico
test_SO <- cox.zph(fit1) #Preizkus sorazmerne ogroženosti
test_SO
##               chisq df       p
## Likvidnost    0.499  1 0.48008
## Oprijemljiva  0.189  1 0.66361
## LnTA          0.349  1 0.55478
## PanogaF       4.315  2 0.11561
## TipF          3.348  1 0.06727
## SubvencijaF  13.141  1 0.00029
## GLOBAL       30.127  7   9e-05
fit1 <- coxph(prež_obj ~ Likvidnost + Oprijemljiva + LnTA + PanogaF + TipF, 
              data = podatki)
test_SO <- cox.zph(fit1)
test_SO
##                 chisq df     p
## Likvidnost   6.17e-04  1 0.980
## Oprijemljiva 1.23e-01  1 0.726
## LnTA         1.30e+00  1 0.254
## PanogaF      5.30e+00  2 0.070
## TipF         2.03e+00  1 0.154
## GLOBAL       1.20e+01  6 0.061
library(survminer) #Aktiviramo knjižnico
ggcoxzph(test_SO[1]) #Grafičen prikaz preizkusa sorazmerne ogroženosti

ggcoxzph(test_SO[2])

ggcoxzph(test_SO[3])

ggcoxzph(test_SO[4])

ggcoxzph(test_SO[5])

library(rms) #Aktiviramo knjižnico
## Loading required package: Hmisc
## Loading required package: lattice
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## The following objects are masked from 'package:DescTools':
## 
##     %nin%, Label, Mean, Quantile
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: SparseM
## 
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
## 
##     backsolve
rms::vif(fit1) #Izračun VIF statistike pri Coxovi regresiji
##       Likvidnost     Oprijemljiva             LnTA PanogaFProizvod. 
##         1.074286         1.190572         1.271263         1.353608 
##    PanogaFGradb.         TipFd.d. 
##         1.354335         1.238082
mean(rms::vif(fit1))
## [1] 1.247024
#Izračun standardiziranih ostankov
podatki$StdOst <- residuals(fit1, 
                            type = "deviance")
hist(podatki$StdOst,
     main = "Histogram standardiziranih ostankov",
     ylab = "Frekvenca",
     xlab = "Standardizirani ostanki")

head(podatki[order(-podatki$StdOst), c("ID", "StdOst")])
##        ID   StdOst
## 1473 1473 4.250923
## 1594 1594 4.241874
## 143   143 4.011675
## 785   785 3.965002
## 2084 2084 3.818250
## 377   377 3.725538
head(podatki[order(podatki$StdOst), c("ID", "StdOst")])
##        ID    StdOst
## 490   490 -2.253672
## 428   428 -1.977296
## 1842 1842 -1.787661
## 264   264 -1.769310
## 515   515 -1.696005
## 291   291 -1.683889
which(podatki$StdOst > 3) #Katere enote imajo standardizirane ostanke nad 3?
##  [1]  143  377  413  645  685  762  785 1018 1473 1594 1609 2084 2323
podatki <- podatki[c(-143, -377, -413, -645, -685, -762, -785, -1018, -1473, -1594, -1609, -2084, -2323), ]
prež_obj <- Surv(podatki$Delovanje, podatki$Stečaj)  #Ponovno ocenimo preživetveni objekt

fit0 <- coxph(prež_obj ~ 1, #Model brez pojasnjevalnih spremenljivk
              data = podatki)

fit1 <- coxph(prež_obj ~ Likvidnost + Oprijemljiva + LnTA + PanogaF + TipF, #Končni model
              data = podatki)
anova(fit0, fit1, test="Chi") #Primerjava modelov s pomočjo razmerja verjetij
## Analysis of Deviance Table
##  Cox model: response is  prež_obj
##  Model 1: ~ 1
##  Model 2: ~ Likvidnost + Oprijemljiva + LnTA + PanogaF + TipF
##    loglik  Chisq Df P(>|Chi|)    
## 1 -4212.7                        
## 2 -3928.2 568.93  6 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(fit1) #Izpis rezultatov končnega modela
## Call:
## coxph(formula = prež_obj ~ Likvidnost + Oprijemljiva + LnTA + 
##     PanogaF + TipF, data = podatki)
## 
##   n= 2341, number of events= 557 
## 
##                       coef exp(coef)  se(coef)       z Pr(>|z|)    
## Likvidnost       -1.294148  0.274131  0.089159 -14.515  < 2e-16 ***
## Oprijemljiva     -0.015238  0.984877  0.002044  -7.455 9.01e-14 ***
## LnTA             -0.254329  0.775437  0.031485  -8.078 6.60e-16 ***
## PanogaFProizvod.  0.588927  1.802053  0.099970   5.891 3.84e-09 ***
## PanogaFGradb.     1.011429  2.749528  0.115243   8.776  < 2e-16 ***
## TipFd.d.          0.756072  2.129893  0.117014   6.461 1.04e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## Likvidnost          0.2741     3.6479    0.2302    0.3265
## Oprijemljiva        0.9849     1.0154    0.9809    0.9888
## LnTA                0.7754     1.2896    0.7290    0.8248
## PanogaFProizvod.    1.8021     0.5549    1.4814    2.1921
## PanogaFGradb.       2.7495     0.3637    2.1936    3.4463
## TipFd.d.            2.1299     0.4695    1.6934    2.6789
## 
## Concordance= 0.783  (se = 0.009 )
## Likelihood ratio test= 568.9  on 6 df,   p=<2e-16
## Wald test            = 415.7  on 6 df,   p=<2e-16
## Score (logrank) test = 390.8  on 6 df,   p=<2e-16