Analiza zgodovine dogodkov
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 0Opis spremenljivk:
- ID: Oznaka podjetja.
- Delovanje: Število mesecev delovanja podjetja v obravnavanem obdobju.
- Stečaj: 0:NE, 1:DA.
- Likvidnost: Razmerje med kratkoročnimi sredstvi in kratkoročnimi obveznostmi.
- Oprijemljiva: Odstotek oprijemljivih sredstev v celotnih sredstvih.
- LnTA: Naravni logaritem celotnih sredstev.
- Panoga: 0:storitve, 1:proizvodnja, 2:gradbeništvo.
- Tip: 0:d.o.o., 1:d.d.
- Subvencija: Je podjetje prejelo državno subvencijo? 0:NE, 1:DA.
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.762KM_fit #Opisna statistika ## Call: survfit(formula = prež_obj ~ 1)
##
## n events median 0.95LCL 0.95UCL
## [1,] 2354 570 NA NA NAlibrary(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 NAautoplot(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-16summary(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.566KM_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 NAautoplot(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-05Coxova 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= 2354fit1 <- 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-05fit1 <- 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.061library(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':
##
## backsolverms::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.238082mean(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.725538head(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.683889which(podatki$StdOst > 3) #Katere enote imajo standardizirane ostanke nad 3?## [1] 143 377 413 645 685 762 785 1018 1473 1594 1609 2084 2323podatki <- 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 ' ' 1summary(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