# Survival model with frailty term.

# Ardo, UCL 2013

# Prelim:
library(statmod)
library(msm)
digits<-3

#####################################
# Gaussian Quadrature for integration:
nnodes<- 25
method<-"GQ"
quad<-gauss.quad(nnodes,"legendre")
nodes<-quad$nodes
weights<-quad$weights
# Integral function:
integraal<-function(integrand,mu,range){
     a<- mu-range;  b<- mu+range
     x<-(b-a)/2*nodes+(a+b)/2
     f<-rep(NA,nnodes)
     for(i in 1:nnodes){f[i]<-integrand(x[i])}
     return((b-a)/2*weights%*%f)
}
# Width parameter for integration call:
ww<-4

# Prelim MsM:
eps<-1/12
dead<-2
subjects<-as.numeric(names(table(data$id)))
N<-length(subjects)
J<-5
cat("\nState table :\n")
print(statetable.msm(state,id,data=data))
centre.id<-data[data$wave1==1,]$id.centre

# To avoid precision problems, a big number is used in the likelihood:
BigNumber<-2000

# Prepare data:
da <- new.env()
for(i in 1:N){
  data.i<-data[data$id==subjects[i],]
  eval(parse(text=paste(paste("da$d",as.character(i),sep=""),"<-data.i",sep="")))
}
da<-as.list(da)

# Likelihood:
loglikelihood<-function(p){

 # Intensities parameters:
 beta<-as.vector(p[1:5])
 sigma<-exp(p[6])
  
 # Contribution per unit:
 loglik<-0
 for(j in 1:J){
  index<- which(centre.id==j)
  integrand<-function(u.j){
    loglik.j<-0
    for(i in index){
     # Data for subject ij:
     eval(parse(text=paste("data.i<-",paste("da$d",as.character(i),sep=""),sep="")))
     O<-data.i$state; L.i<-length(O)
     t<-data.i$age-60
     cova<-cbind(1,t,data.i$educ,data.i$ybrth-1900,data.i$educ*t)
  
    # Likelihood contribution:
    if(L.i>2){
     for(k in 2:(L.i-1)){
      # Q and P matrix:
      Q<- exp(beta%*%cova[k-1,]+u.j)
      P<- exp(-Q*(t[k]-t[k-1]))
      # No censoring, no death:
      loglik.j<-loglik.j+log( P )
      }
   }
   # Last interval:
   k<-L.i
   Q<- exp(beta%*%cova[k-1,]+u.j)
   if(O[L.i]==dead){
     P<- exp(-Q*(t[k]-eps-t[k-1]))
     P.eps<- exp(-Q*eps)
     # Death:
     loglik.j<-loglik.j+log( (1-P.eps)*P )
   }else{
     P<- exp(-Q*(t[k]-t[k-1]))
     # Censoring:
     loglik.j<-loglik.j+log( P )
     }
   }
   return(exp(loglik.j+BigNumber+log(dnorm(u.j,mean=0,sd=sigma))))
  }
  int<-integraal(integrand,mu=0,range=ww*sigma)
  loglik<-loglik+log(int)-BigNumber
  }
  # Monitoring:
  cat("-2*Loglik = ", -2*loglik," with sigma = ",sigma,"\n")
  -loglik
}

# Maximise and report:
p<-c( -5,0.2,-0.4,0,0,-2)
opt.method<-"Nelder-Mead"
#max<-optim(par=p, fn=loglikelihood, method = opt.method,control=list(maxit=20000),hessian=TRUE)
cat("\n-2loglik =", 2*max$val,"\n")
cat("Optimisation method =", opt.method,"\n")
cat("Integration by ", method,"with",nnodes," nodes, and width parameter",ww,"\n")
cat("Convergence code =", max$con,"\n")
cat("Epsilon in msm model =", eps,"\n")
cat("Parameters and SEs:\n")
p<-max$par
p.se<-sqrt(diag(solve(max$hessian)))
print(round(cbind(p,p.se,"  Wald ChiSq"=(p/p.se)^2,"Pr>ChiSq"=1-pchisq((p/p.se)^2,df=1)),digits))
L<-length(p)
# Sigma:
sigma<-exp(p[L])
# Delta method for SE:
sigma.se<- p.se[L]*exp(p[L])
cat("sigma = ",round(sigma,digits),"with SE",round(sigma.se,digits),"\n")


# MAP estimation of frailties:
beta<-as.vector(p[1:5])
sigma<-exp(p[6])
cat("MAP estimation of frailties...\n")
frailty<-rep(NA,J)
info<-rep(NA,J)
for(j in 1:J){
index<- which(centre.id==j)
loglikelihood.j<-function(p){
   u.j<-p
   loglik.j<-0
   for(i in index){
     # Data for subject ij:
     eval(parse(text=paste("data.i<-",paste("da$d",as.character(i),sep=""),sep="")))
     O<-data.i$state; L.i<-length(O)
     t<-data.i$age-60
     cova<-cbind(1,t,data.i$educ,data.i$ybrth-1900,data.i$educ*t)
  
    # Likelihood contribution:
    if(L.i>2){
     for(k in 2:(L.i-1)){
      # Q and P matrix:
      Q<- exp(beta%*%cova[k-1,]+u.j)
      P<- exp(-Q*(t[k]-t[k-1]))
      # No censoring, no death:
      loglik.j<-loglik.j+log( P )
      }
   }
   # Last interval:
   k <- L.i
   Q<- exp(beta%*%cova[k-1,]+u.j)
   if(O[L.i]==dead){
     P<- exp(-Q*(t[k]-eps-t[k-1]))
     P.eps<- exp(-Q*eps)
     # Death:
     loglik.j<-loglik.j+log( (1-P.eps)*P )
   }else{
     P<- exp(-Q*(t[k]-t[k-1]))
     # Censoring:
     loglik.j<-loglik.j+log( P )
     }
  }
  loglik.j<-loglik.j+log(dnorm(u.j,mean=0,sd=sigma))
  return(-loglik.j)
}
start.par<-0
max.j<-optim(par=start.par, fn=loglikelihood.j, method = "BFGS",control=list(trace=0,maxit=20000),hessian=FALSE)
frailty[j]<-max.j$par
info[j]<-max.j$convergence
}
tab<-frailty
names(tab)<-centres
print(round(tab,digits))