library(survival)

#***************************** Feigl and Zelen
 fz <- read.table("~/LifeData/FZ/feiglz.dat", header=T)
attach(fz)
 dead <- rep(1,33) ; lwbc <- log(wbc) ; y <- log(tim) 
 
 #  plots with continuous variables: 
 summary(wbc) ; (b <- summary(wbc)[c(1,2,3,5,6)] )
 wbcGp <- cut(wbc, breaks=b, include.lowest=T, labels=(c("q1","q2","q3","q4"))) 
  # provide a variable which indicates which quartile wbc values are in.
  i <- order(wbc) ;  cbind(wbc[i],wbcGp[i]) 
  # look at result. 
  medlwgp <- c(median(lwbc[wbcGp=="q1"]),  median(lwbc[wbcGp=="q2"]), median(lwbc[wbcGp=="q3"]), median(lwbc[wbcGp=="q4"]))
 
 km.wGp <-  survfit(Surv(tim, dead) ~ wbcGp, conf.type="log-log")
 plot(km.wGp, lty=1:4, xlab="Time (weeks)", ylab="Estimated survival", lwd=1.5, main="Survival of leukaemia patients") 
 legend(80,1, c("Lowest quartile","Second quartile","Third  quartile","Highest quartile"),lty=1:4, cex=1.5, bty="n")
# This is how to look at the general pattern of estimated survival functions with a continous variable. 
 km.wGp
 
  plot(km.wGp, lty=1:4, fun="cloglog" ,  xlab="Time (weeks)", ylab="Estimated survival", lwd=1.5, main="Survival of leukaemia patients")  
  legend(25,-1, c("Lowest quartile","Second quartile","Third  quartile","Highest quartile"),lty=1:4, cex=1.1, bty="n")
  
  # Fitted lines
       t <- 1:140
  fz.lw <- survreg(Surv(tim, dead) ~ lwbc,dist="exponential")
    t <- 1:140  
    fzlamq <- exp(-7.3755014 + 0.4012504 * medlwgp)
     # sets up four vectors of fitted lines. 
 lines(t, exp(-fzlamq[1]*t ),lty=1); lines(t, exp(-fzlamq[2]*t ),lty=2)  
  lines(t, exp(-fzlamq[3]*t ),lty=3); lines(t, exp(-fzlamq[4]*t ),lty=4)
 
 ( fz.lwgp <- survreg(Surv(tim, dead) ~ wbcGp,dist="exponential") )
  summary(fz.lwgp)
  
  (mb <- c(b[1],b[3], b[5]) )
 wbcMed <- cut(wbc, breaks=mb, include.lowest=T, labels=(c("m1","m2"))) 
  # provide a variable which indicates which quartile wbc values are in.
  c(median(wbc[wbcMed=="m1"]),  median(wbc[wbcMed=="m2"]))
  # look at result. 
 fz.lwMi <- survreg(Surv(tim, dead) ~ ag*wbcMed,dist="exponential")
   # with interaction.
 fz.lwM<- survreg(Surv(tim, dead) ~ ag+wbcMed,dist="exponential")
  # without interaction. Loglik(model)= -143.9   Loglik(intt only)= -155.5
  fz.lwag <- survreg(Surv(tim, dead) ~ ag+lwbc,dist="exponential")
  #  Loglik(model)= -146.5

 # ***************************** Look at graphs of this. 
  
 km.wMag <-  survfit(Surv(tim, dead) ~ ag+wbcMed, conf.type="log-log")
 plot(km.wMag, lty=c(3,2,4,1), xlab="Time (weeks)", ylab="Estimated survival", lwd=1.5, main="Survival of leukaemia patients") 
 legend(80,1, c("AG-, M1","AG-, M2","AG+, M1","AG+, M2"),lty=c(3,2,4,1), cex=1.1, bty="n")
# This is how to look at the general pattern of estimated survival functions with a continous variable. 
 km.wGp
   fz.lwM$coef    # coefficients of the model. 
   (fz.lwM$coef%*%c(1,0,0)) # linear predictor for ag=-ve, wbc = m1.     
  (fz.lwM$coef%*%c(1,1,1)) # linear predictor for ag=+, wbc = m2. 

  lines(t, exp (- exp(-fz.lwM$coef%*%c(1,0,0)) *t ),lty=3)
  lines(t, exp (- exp(-fz.lwM$coef%*%c(1,0,1))*t ),lty=2)  
  lines(t, exp (- exp(-fz.lwM$coef%*%c(1,1,0))*t ),lty=4)
  lines(t, exp (- exp(-fz.lwM$coef%*%c(1,1,1))*t ),lty=1)
 
  plot(km.wMag, lty=1:4, fun="cloglog" ,  xlab="Time (weeks)", ylab="Estimated survival", lwd=1.5, main="Survival of leukaemia patients")  
  legend(1,0.7,c("AG-, M1","AG-, M2","AG+, M1","AG+, M2"),lty=c(3,2,4,1),  cex=1.1, bty="n")
  
 
# ***************************** A summary of these results. 
  
      
                records n.max n.start events median 0.95LCL 0.95UCL
ag=0, wbcMed=m1       7     7       7      7     17       3      56  23.7
ag=0, wbcMed=m2       9     9       9      9      4       2      30   6.6
ag=1, wbcMed=m1      10    10      10     10    104      16     134  56.3
ag=1, wbcMed=m2       7     7       7      7      5       1      26  15.8 
 
   
#        
#
#

m00 <-  exp(fz.lwM$coef%*%c(1,0,0)) #estimated mean for ag=-ve, wbc = m1.
  log(2)*m00 #estimated mean for ag=-ve, wbc = m1.Median 



ins <- read.table("~/public_html/Insul.csv",sep=",",header=T)
 attach(ins)  ; names(ins)
  kVolt <- as.factor(kV) # this useful necessary for checking models. 
  icens <- rep(1,length(ins$t))
 

insa.wei <- survreg(Surv(ins$t,icens)~kV,dist="weib")
insfa.wei <- survreg(Surv(ins$t,icens)~kVolt,dist="weib")
 plot(seq(28,38,2),insfa.wei$coef[2:7]) ;  kv<- seq(28,38,2)
 -0.5544469