library(survival)
attach(aml)


#  Excercise 3 series is actually not the same as AML data:
#   5   5   8   8+  12  16+ 23  27  30  33  43  


t.i <- c(5,5,8,8,12,16,23,27,30,33,43)
st  <- c(1,1,1,0,1,0,1,1,1,1,1)
exc3.surv <-  Surv(t.i,st)


# 3.1 ALEPH


# Log-lik 2nd derivatives according to Weibull parameterization

#  d^2 log-L / d alpha^2:
 - n.u / alpha^2 - sum( (lambda*t.i)^alpha * (log(lambda*t.i))^2)

#  d^2 log-L / d alpha d lambda:
   n.u / lambda - sum(alpha*lambda^(alpha-1)*t.i^alpha*log(lambda*t.i)) - sum((lambda*t.i)^alpha)/lambda

#  d^2 log-L / d lambda^2:
   -n.u*alpha/lambda^2 - sum((alpha-1)*alpha*lambda^(alpha-2)*t.i^alpha)

# To derive tests that correspond to the tests implemented in R we express  	
# Log-lik 2nd derivatives according to R parameterization:  theta = -log(alpha) and mu = -log(lambda):

#  d^2 log-L / d theta^2:
 (- n.u / alpha^2 - sum( (lambda*t.i)^alpha * (log(lambda*t.i))^2)) * alpha^2

#  d^2 log-L / d theta d mu:
  ( n.u / lambda - sum(alpha*lambda^(alpha-1)*t.i^alpha*log(lambda*t.i)) - sum((lambda*t.i)^alpha)/lambda) * alpha * lambda

#  d^2 log-L / d mu^2:
  ( -n.u*alpha/lambda^2 - sum((alpha-1)*alpha*lambda^(alpha-2)*t.i^alpha)) * lambda^2




#  Wald statistic computation


expo.fit <- survreg(exc3.surv~1,dist="expo")
weib.fit <- survreg(exc3.surv~1,dist="weibull")

summary(weib.fit)
weib.fit$var
sqrt(diag(weib.fit$var))
solve(weib.fit$var)

(n.u    <-  sum(st))
(alpha  <-  1 / weib.fit$scale)
(lambda <-  exp(-weib.fit$coef))


#  d^2 log-L / d theta^2:
 (- n.u / alpha^2 - sum( (lambda*t.i)^alpha * (log(lambda*t.i))^2)) * alpha^2

#  d^2 log-L / d theta d mu:
  ( n.u / lambda - sum(alpha*lambda^(alpha-1)*t.i^alpha*log(lambda*t.i)) - sum((lambda*t.i)^alpha)/lambda) * alpha * lambda

#  d^2 log-L / d mu^2:
  ( -n.u*alpha/lambda^2 - sum((alpha-1)*alpha*lambda^(alpha-2)*t.i^alpha)) * lambda^2


# Score statistic computation

(n.u    <-  sum(st))
(alpha  <-  1 )
(lambda <-  exp(-expo.fit$coef))

# score statistic - d log-L / d theta = d log-L / d alpha * ( d alpha / d theta) :
(score.stat  <-  (n.u / alpha + n.u * log(lambda) + sum(log(t.i)[st==1]) - sum( (lambda*t.i)^alpha * log(lambda*t.i))) * ( - alpha))

# Log-lik 2nd derivatives:

#  d^2 log-L / d theta^2:
( i.22 <-  - (- n.u / alpha^2 - sum( (lambda*t.i)^alpha * (log(lambda*t.i))^2)) * alpha^2)

#  d^2 log-L / d theta d mu:
( i.12 <-  - ( n.u / lambda - sum(alpha*lambda^(alpha-1)*t.i^alpha*log(lambda*t.i)) - sum((lambda*t.i)^alpha)/lambda) * alpha * lambda)

#  d^2 log-L / d mu^2:
( i.11 <-   -  ( -n.u*alpha/lambda^2 - sum((alpha-1)*alpha*lambda^(alpha-2)*t.i^alpha)) * lambda^2)

# information matrix:

(i.m <- cbind(c(i.11,i.12),c(i.12,i.22)))

score.stat / sqrt(i.m[2,2])



# LR test


summary(weib.fit)
summary(expo.fit)

anova(expo.fit,weib.fit)




# 3.2 BET

# Weibull qqplot


surv.fit <- survfit(exc3.surv)
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))
summary(weib.fit)




# 3.2 GIMEL

# exponential model CI for S(30) and h(30)

summary(expo.fit)

# est and CI for lambda (i.e. for h(30) )

(est <- exp(- expo.fit$coeff ))
(ul  <- exp(- (expo.fit$coeff - qnorm(.975)*sqrt(expo.fit$var)))[1,1])
(ll  <- exp(- (expo.fit$coeff + qnorm(.975)*sqrt(expo.fit$var)))[1,1])


# est and CI for S(30)

exp(-est*30)
exp(-ul*30)
exp(-ll*30)