library(survival)
rm(x)
rm(time)
rm(status)
attach(aml)


#-------------------------




par(mfcol=c(1,2))
x <- seq(.0001,3,length=40)
plot(x,dexp(x,3),type="l",ylab="f(t)",xlab="t")
lines(x,dexp(x,1))

lines(x,dexp(x,.5))
text(.5,2.0,"Lambda=3")
text(.2,1.1,"Lambda=1")
text(.2,0.35,"Lambda=0.5")


x <- seq(.0001,6,length=40)
plot(x,1-pexp(x,.5),ylim=c(0,1),type="l",ylab="S(t)",xlab="t")

lines(x,1-pexp(x,1))


lines(x,1-pexp(x,3))
text(.1,0.1,"Lambda=3")
text(2,0.7,"Lambda=0.5")


#----    plot Weibull pdf and hazard ----

par(mfcol=c(1,2))

x <- seq(.0001,4,length=40)
plot(x,dweibull(x,shape=3,scale=1),type="l",ylab="f(t) - scale=1",xlab="t")
lines(x,dweibull(x,shape=1.5,scale=1))
lines(x,dweibull(x,shape=1,scale=1))
lines(x,dweibull(x,shape=0.5,scale=1))
text(1.5,1.1,"shape=3")
text(0.5,0.25,"shape=0.5")

plot(x,3*x^2,type="l",ylab="h(t) - scale=1",xlab="t",ylim=c(0,2.5))
lines(x,1.5*x^.5)
lines(x,x/x)
lines(x,.5*x^(-.5))
text(1.2,2,"shape=3")
text(2.5,2,"shape=1.5")
text(3,0.25,"shape=0.5")


#----  plot standard extreme value pdf and survival ----
 y <- seq(-4,2,length=40)

par(mfcol=c(1,2))

plot(y,exp(y-exp(y)),type="l",ylab="standard extreme value f(t)",xlab="t")
plot(y,exp(-exp(y)),type="l",ylab="standard extreme value S(t)",xlab="t")


#----  plot Log-logistic pdf and hazard ----


x <- seq(.0001,3,length=40)

d.log.log <- function(x,shape=1,scale=1)
shape*scale* (scale*x)^(shape-1) * (1+(scale*x)^shape)^(-2)
 h.log.log <- function(x,shape=1,scale=1)
shape*scale* (scale*x)^(shape-1)  / ( 1 + (scale*x)^shape)

par(mfcol=c(1,2))
plot(x,d.log.log(x,shape=3,scale=1),type="l",ylab="f(t) - scale=1",xlab="t",ylim=c(0,2),main="Density")
lines(x,d.log.log(x,shape=1.5,scale=1))
lines(x,d.log.log(x,shape=.5,scale=1))
lines(x,d.log.log(x,shape=.25,scale=1))
text(1.5,0.0,"shape=3")
text(1.5,0.15,"shape=1.5")
text(1.5,0.4,"shape=0.5")
text(1.5,0.8,"shape=0.25")

plot(x,h.log.log(x,shape=3,scale=1),type="l",ylab="h(t) - scale=1",xlab="t",ylim=c(0,6),main="Hazard")
lines(x,h.log.log(x,shape=1.5,scale=1))
lines(x,h.log.log(x,shape=.5,scale=1))
lines(x,h.log.log(x,shape=.25,scale=1))
text(1.5,0,"shape=3")
text(1.5,1.8,"shape=0.25")

#----  plot Gamma and log-Gamma ----

h.gamma <- function(x = 1:10, shape = 1, scale = 1)
{
	f.g <- dgamma(x, shape, scale)
	S.g <- 1 - pgamma(x, shape, scale)
	return(f.g/S.g)
}

d.log.gamma <- function(z = 1:10, shape = 1)
	exp(shape * z - exp(z))  /gamma(shape)

x <- seq(.0001,5,length=40)
 par(mfcol=c(1,2))
plot(x,dgamma(x,shape=3,scale=1),type="l",ylab="f(t)- scale=1",xlab="t",ylim=c(0,1.2))
lines(x,dgamma(x,shape=2,scale=1))
lines(x,dgamma(x,shape=1,scale=1))
lines(x,dgamma(x,shape=0.5,scale=1))
text(4,0.2,"shape=3")
text(1,1.1,"shape=0.5")

plot(x,h.gamma(x,shape=3,scale=1),type="l",ylab="h(t)- scale=1",xlab="t",ylim=c(0,2))
lines(x,h.gamma(x,shape=2,scale=1))
lines(x,h.gamma(x,shape=1,scale=1))
lines(x,h.gamma(x,shape=0.5,scale=1))
text(3.5,0.4,"shape=3")
text(0.8,0.6,"shape=2")
text(0.5,1.1,"shape=1")
text(1,1.8,"shape=0.5")

par(mfcol=c(1,1))

x <- seq(-4,3,length=40)

plot(x,d.log.gamma(x,shape=2),type="l",ylab="f(z)- scale=1",xlab="z",ylim=c(0,0.6))

lines(x,d.log.gamma(x,shape=1))

lines(x,d.log.gamma(x,shape=0.5))
text(1.8,0.4,"shape=2")
text(-1.2,0.35,"shape=1")
text(-3.5,0.15,"shape=0.5")



# ----  QQ plot example 1:  qqnorm ----

par(mfcol=c(1,1))
x <- rnorm(12,mean=5,sd=3)
print(round(x,2))
(a <- qqnorm(x))
abline(5,3)

sort(a$x)
qnorm(ppoints(12))

# ----  weibull QQ plot for aml1 data ----

attach(aml)

surv.fit <- survfit(Surv(time[x=="Maintained"],status[x=="Maintained"]),type="kaplan-meier")
summary(surv.fit)
surv.vec <-  summary(surv.fit)$surv
time.vec <- summary(surv.fit)$time
p.vec <- 1 - (surv.vec + c(1,surv.vec[-length(surv.vec)]))/2

y.vec <- log(time.vec)
x.weibull <- log(qweibull(p.vec,shape=1,scale=1))
plot(x.weibull,y.vec)
abline(lm(y.vec~x.weibull))
summary(lm(y.vec~x.weibull))




#---- difference between QQ plots ------

y.weib <- sort(log(rweibull(200,shape=1,scale=1)))
y.log.norm <- sort(rnorm(200,mean=0,sd=1))
y.log.logist <- sort(rlogis(200,location=0,scale=1))

p.vec <- (2*(1:200)-1)/400
x.weib <-     log(qweibull(p.vec,shape=1,scale=1))
x.log.norm <- qnorm(p.vec,mean=0,sd=1)
x.log.logist <- qlogis(p.vec,location=0,scale=1)
 
ylm <- range(y.weib,y.log.norm,y.log.logist)

par(mfrow=c(1,1))
 plot(x.weib,y.weib,ylim=ylm,xlab="Extreme value quantiles",ylab="empirical quantiles",
			main="Weibull fit",type="l",col=1)
lines(x.weib,y.log.norm,col=2)
lines(x.weib,y.log.logist,col=3)
abline(0,1)
legend(-6,4,legend=c("Weibull sample","Log-norm sample","Log-logist sample"),lty=1,col=1:3)


plot(x.log.norm,y.weib,ylim=ylm,xlab="Normal quantiles",ylab="empirical quantiles",
			main="Log-normal fit",type="l",lty=1)
lines(x.log.norm,y.log.norm,lty=2)
lines(x.log.norm,y.log.logist,lty=3)

abline(0,1)
legend(-3,4,legend=c("Weibull sample","Log-norm sample","Log-logist sample"),lty=1:3)

plot(x.log.logist,y.weib,ylim=ylm,xlab="Logistic quantiles",ylab="empirical quantiles",
			main="Log-logistic fit",type="l",lty=1)
lines(x.log.logist,y.log.norm,lty=2)
lines(x.log.logist,y.log.logist,lty=3)

abline(0,1)

legend(-4,4,legend=c("Weibull sample","Log-norm sample","Log-logist sample"),lty=1:3)
			



# -----  Computations for exponential MLE ------


sum.t <- sum(time[x=="Maintained"])
n <- length(time[x=="Maintained"])
ul <- qchisq(0.975,22)
ll <- qchisq(0.025,22)

weeks[group==1]

n/sum.t
n/sum.t * ll / (2*n)
n/sum.t * ul / (2*n) 

sum.t / n 
(n/sum.t * ll / (2*n))^(-1)
(n/sum.t * ul / (2*n))^(-1) 




# ---- Weibull MLE 3D plots -----


Surv(time[x=="Maintained"],status[x=="Maintained"])
u.vec <- time[x=="Maintained" & status == 1]
c.vec <- time[x=="Maintained" & status == 0]

aml1.weib.log.lik <- function(alpha,lambda)
{
	sum(log(dweibull(u.vec,shape=alpha,scale=1/lambda)))+
	sum(log(1-pweibull(c.vec,shape=alpha,scale=1/lambda)))
}

aml1.weib.prof.lik <- function(alpha)
( length(u.vec) / sum(c(u.vec,c.vec)^alpha))^(1/alpha)
 alpha.ser <- seq(0.2,3,length=40)


lambda.ser <- seq(0.002,0.06,length=40)


log.lik.mat <- array(NA,dim=c(40,40))
for(i in 1:40) for(j in 40:1) 
	log.lik.mat[i,j] <- aml1.weib.log.lik(alpha.ser[i],lambda.ser[j])

contour(alpha.ser,lambda.ser,log.lik.mat,levels=c(-50,-45,-40,-37,-36))

prof.lik.ser <- rep(NA,40)
for(i in 1:40)
	prof.lik.ser[i] <- aml1.weib.prof.lik(alpha.ser[i])
lines(alpha.ser,prof.lik.ser)	


# ---- Weibull CI -----

contour(alpha.ser,lambda.ser,-2*(log.lik.mat+35.7),	levels=round(c(qchisq(0.80,2),qchisq(0.95,2)),2))

	
#---- 3 MLE for AML1 data ----


summary(survreg(Surv(time[x=="Maintained"],status[x=="Maintained"])~1,dist="weib"))
a <- summary(survreg(Surv(time[x=="Maintained"],status[x=="Maintained"])~1,dist="weib"))

a$coef  		# mu parameter
exp(-a$coef)	# lambda	

a$scale		# sigma parameter


summary(survreg(Surv(time[x=="Maintained"],status[x=="Maintained"])~1,dist="loglogistic"))
summary(survreg(Surv(time[x=="Maintained"],status[x=="Maintained"])~1,dist="lognormal"))



plot(survfit(Surv(time[x=="Maintained"],status[x=="Maintained"])~1),conf.int=F)
x.ser <- seq(0,161,length=30)
lines(x.ser,1-pweibull(x.ser,shape=exp(0.5),scale=exp(3.64)),lty=2)
lines(x.ser,1-pweibull(x.ser,shape=exp(0.0314),scale=exp(4.0997)),lty=2)
lines(x.ser,1-plogis(log(x.ser),location=3.515,scale=0.542),lty=3)
lines(x.ser,1-pnorm(log(x.ser),mean=3.61,sd=0.961),lty=4)

legend(100,.9,c("K-M","Weibull - QQ","Weibull - MLE","Loglogistic","Lognorm"),lty=c(1,2,2,3,4))




#---- return to the QQ plots ------
 

surv.fit <- survfit(Surv(weeks[group==1],status[group==1]),type="kaplan-meier")
summary(surv.fit)

surv.vec <- 

summary(surv.fit)$surv

time.vec <- 

summary(surv.fit)$time
p.vec <- 1 - (surv.vec + c(1,surv.vec[-length(surv.vec)]))/2
 y.vec <- log(time.vec)
x.weibull <- log(qweibull(p.vec,shape=1,scale=1))

x.loglog  <- qlogis(p.vec,location=0,scale=1)

x.lognorm <- qnorm(p.vec,mean=0,sd=1)

par(mfcol=c(1,3))
plot(x.weibull,y.vec) 

abline(4.0997,0.969)

plot(x.loglog,y.vec) 

abline(3.515,0.542)

plot(x.lognorm,y.vec) 

abline(3.61,0.961)
par(mfcol=c(1,1))