### ref: An Introduction to Statistical Modeling of Extreme Values
### Stuart Coles, pp. 38-43

#### This is a R code for log transformed failure time model. 
#### After log transformation, lamda=aw^b model is still not equal to
#### accelerated failure time model of corrison, but it will equal to 
#### accelerated failure time model of log transformed corrison. 

#### I also want to show log transformation does not change M.L.E and 
#### Deviance.

library(ismev)
library(stats4)

# data from the ismev library
# 

data(engine)

ft = engine$Time    # failure time
n = length(ft)      # sample size
w = engine$Corrosion   # corrosion measurement

plot(w,ft)

log.ft=log(ft)
log.w=log(w)

# MLE

# solve the nonlinear equation 

# minus the log likelihood

ll = function(a,b) { 
  -n*log(a) - b*sum(log(w))- sum(log.ft) + a*sum((w^b)*exp(log.ft))
}

model.mle = mle(minuslog=ll,start=list(a=1,b=1)) 
model.mle

a.mle = coef(model.mle)[1]
a.mle
b.mle = coef(model.mle)[2]
b.mle

ll.value = -ll(a.mle,b.mle)
ll.value

# CIs for a and b, the CI for b matches with accelerated ft model 
# of transformed corrosion level.  

I.O = matrix(c( n*a.mle^-2, sum(w^b.mle*ft*log(w)),
sum(w^b.mle*ft*log(w)), a.mle*sum(w^b.mle*ft*log(w)^2) ), c(2,2))
I.O 

I.O.inv = solve(I.O)  # produces the matrix inverse
I.O.inv

a.mle.se = sqrt(I.O.inv[1,1])
a.mle.se

conf.level = 0.95
cv = qnorm(1-(1-conf.level)/2)

a.ci = c(a.mle - cv*a.mle.se,a.mle + cv*a.mle.se)
a.ci


b.mle

b.mle.se = sqrt(I.O.inv[2,2])
b.mle.se

b.ci = c(b.mle - cv*b.mle.se,b.mle + cv*b.mle.se)
b.ci 


# Test if H0: b=0, H1: b=/= 0 
a.mle.small = n/sum(ft)
a.mle.small

# log-likelihood of null hypothysis
ll.value.small = n*log(a.mle.small)+sum(log.ft)-a.mle.small*sum(exp(log.ft))
ll.value.small

# Deviance is the same with the model before the log transformation
D = 2*(ll.value - ll.value.small)
D


# plot fitted model, I am not sure this is right.

a.mle = coef(model.mle)[1]
a.mle
b.mle = coef(model.mle)[2]
b.mle

w.index = seq(min(w),max(w),0.01)

ft.fit = (a.mle^-1)*w.index^(-b.mle)
log.ft.fit=log(ft.fit)

plot(w,log.ft,xlab="Corrosion Level",ylab="Transformed Lifetime",main="Fitted Model")
lines(w.index,log.ft.fit,type="l",col=3)

# From Model fit aspect, the transformed failure time model does not perform
# as good as the original model. The reason of using this model is to illustrate
# the log-transformation only shifts the likelihood. It does not change deviance
# therefore, the inference made from model will not change.