Nothing
#######################
# observable survival
########################
my.survnp <- function(tt,st, D,ratetable,konc,conf.int=0.95){
md.init(ratetable)
#estimates event-free (EFS) and net survival
#tt ... cas
#st ... censoring indicator (1=dogodek, nic sicer)
#R ... spremenljivke za populacijsko tabelo
#ratetable ... uporabljena populacijska tabela
#konc ... maximum potential observation time for this individual
if(missing(konc)) konc <- rep(max(tt),length(tt)) #ce ni podan koncni datum krnjenja - das vsem na zadnjega opazenega
if(any((tt>konc)&(st==1)))stop("Events are observed after administrative censoring \n") #ce imamo za nekoga podan konc, ne more imeti dogodka po tem
if(any((tt!=konc)&(st==0)))stop("Only administrative censoring is allowed - max times should be reported for all individuals \n") #za nekoga imamo cas krnjenja, ki ni enak koncu
data <- data.frame(time=round(tt,8),cens=st,konc=konc)
data <- data[order(tt),] #uredim po casih
D <- D[order(tt),] #uredim po casih
#ttun <- sort(unique(c(data$time,data$konc))) #unique casi dogodkov ali censoringa
ttun <- sort(unique(c(data$time,data$konc,1:max(data$time,data$konc)))) #unique casi dogodkov ali censoringa - za net je vazno, da jih je cimvec, za ostale je vseeno
fk = (attr(ratetable, 'type') != 1)
#fk <- (attributes(ratetable)$factor != 1) #spremenljivke iz popul. tabel, ki se vecajo v casu
nfk <- length(fk) #stevilo spr. v popul. tabelah
difft <- c(ttun[1],diff(ttun)) #sprememba casa na posameznem intervalu
le.net <- vari.net <- rep(NA,length(ttun)) #hazard, stevilo ljudi ob vsakem casu in prirastek variance
le.efs <- vari.efs <- rep(NA,length(ttun)) #hazard, stevilo ljudi ob vsakem casu in prirastek variance
Yosebni <- rep(1,nrow(data)) #Yosebni je verjetnost, da je se ziv v populaciji
Yosebnin <- Spicum <- rep(1,nrow(data)) #Spicum potrebujemo za net survival, na zacetku 1 za vse, Yosebnin je za net
Y.efs = Y.net = rep(0, length(ttun))
test = rep(NA,length(ttun))
for(it in 1:length(ttun)){
#population survival probabilities
Spi <- md.expprep(D, difft[it]) #populacijska verjetnost, da zivijo se naprej
Spicum <- Spicum*Spi #unconditional (od zacetka do konca), potrebujem za net survival
dLpi <- -log(Spi)
test[it]=Spi[1]
#at risk indicator
Yosebni <- Yosebni*Spi #Yosebni se zmanjsa - nekateri umrejo
#number of events
Nei <- data$time==ttun[it]&data$cens==1 #individual event indicator
who.event <- which(Nei==1) #kdo ima dogodek ob tem casu
Ne <- length(who.event) #stevilo dogodkov ob tem casu
Npi <- Yosebni*dLpi
Np <- sum(Npi)
#net number of events
Neisp <- Nei/Spicum
Nesp <- sum(Neisp) #stevilo dogodkov ob tem casu v hyp world
#hazard - net
Y.net[it] = Yoni <- sum(Yosebni/Spicum) #number at risk in hyp world
le.net[it] <- Nesp/Yoni #trenutno tveganje
vari.net[it] <- sum(Nei/Spicum^2)/Yoni^2 #prirastek k varianci
#hazard - event free
Y.efs[it] = sum(Yosebni)
le.efs[it] <- (Ne+Np)/sum(Yosebni) #trenutno tveganje
vari.efs[it] <- sum((Nei + Npi)^2)/sum(Yosebni)^2 #prirastek k varianci
#preparing for the next time interval
who.cens <- which(data$konc==ttun[it]) #kdo je krnjen oz. konca opazovanje ob tem casu
Yosebni[who.event] <- Yosebni[who.event]-1 #odstejem 1 pri tistemu, ki ima dogodek
Yosebni[who.cens] <- 0 #dokoncno vrzem ven vse, ki imajo konc (tudi ce je imel prej event)
D$age = D$age + difft[it]
}
se.fac <- sqrt(qchisq(conf.int, 1)) #factor needed for confidence interval
#NET
vari.net <- c(0,cumsum(vari.net))
surv.net <- exp(-c(0,cumsum(le.net)))
#log scale conf. interval
lower.net <- exp(as.vector(log(surv.net) - sqrt(vari.net)*se.fac))
upper.net <- exp(as.vector(log(surv.net) + sqrt(vari.net)*se.fac))
#EFS
vari.efs <- c(0,cumsum(vari.efs))
surv.efs <- exp(-c(0,cumsum(le.efs)))
#log scale conf. interval
lower.efs <- exp(as.vector(log(surv.efs) - sqrt(vari.efs)*se.fac))
upper.efs <- exp(as.vector(log(surv.efs) + sqrt(vari.efs)*se.fac))
out <- list(time=c(0,ttun),
surv.net=surv.net,std.err.net=sqrt(vari.net),lower.net=lower.net,upper.net=upper.net,
surv.efs=surv.efs,std.err.efs=sqrt(vari.efs),lower.efs=lower.efs,upper.efs=upper.efs,
Y.net, Y.efs, test=test
)
out
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.