library(survival)
kokoner<-read.table("kokon.txt",header=T)
attach(kokoner)
options('contrasts=contr.treatment')

################### Variabel-def start ###############
treatF<-0*as.numeric(kokoner$treat);treatF[kokoner$treat=='F']<-1
treatI<-0*as.numeric(kokoner$treat);treatI[kokoner$treat=='I']<-1
treatK<-0*as.numeric(kokoner$treat);treatK[kokoner$treat=='K']<-1
treat<-treatF+2*treatI

age.gr1<-0*age+1*(age<10);
age.gr2<-0*age+1*( (age>=10)&(age<15) );
age.gr<-age.gr1+2*age.gr2
age.gr<-2-(2*age.gr1+1*age.gr2)


###################Kaplan-Meier plots #################

surv.all<-survfit(Surv(time,status==1))
par(mfrow=c(1,1))
plot(surv.all,mark.time=F)

surv.treat<-survfit(Surv(time,status==1)~treat)
plot(surv.treat)
plot(surv.treat,conf.int=T)


################## Cox-regression ##################
fit1<-coxph(Surv(time,status==1)~treatF)

########### plotkum-funktion til model kontrol####


plotkum<-function(x){
fit<-x
l<-length(survfit(fit)$strata)#Antal strata
no<-numeric(l)
no[1:l]<-survfit(fit)$strata[1:l];
no1<-cumsum(no[1:l])
xup<-max(survfit(fit)$time)
ylow<<-min(log(-log(survfit(fit)$surv)))
yup<<-max(log(-log(survfit(fit)$surv)))
plot(c(survfit(fit)$time[1:no1[1]],xup*1.05),
     log(-log(survfit(fit)$surv[c(1:no1[1],no1[1])]
              )),type='s',ylim=c(ylow,yup),xlim=c(0,xup),xlab='Time',
                ylab='Log(Cumulative hazard)',lty=c(1),lwd=1.5,main='Log-Cumulative hazard plots')
for (j in 2:l){
     lines(c(survfit(fit)$time[(no1[j-1]+1):no1[j]],xup*1.05),
     log(-log(survfit(fit)$surv[c((no1[j-1]+1):no1[j],no1[j])])),type='s',
      lty=c(j+1),lwd=1.5)
}
legend(xup*0.55,0.5*ylow,names(survfit(fit)$strata)[1:l],lty=c(1:l),lwd=1.5)
}