# MsM model with frailties.

# Ardo, UCL 2013


# Prelim:
digits<-3

# Prelim MsM:
eps<-1/12
dead<-3; censored<- -2
subjects<-as.numeric(names(table(data$id)))
Q<-rbind(c(0,1,1), c(1,0,1),c(0,0,0))
N<-length(subjects)
J<-5
cat("\nState table :\n")
print(statetable.msm(state,id,data=data))
BigNumber<-2000
# Parameter for integration:
ww<-4

# Center id per subject:
index<-which(data$wave1==1)
idN.centre<-data$id.centre[index]

# P matrix given Q for three-state reversible:
matrixexp<-function(Q,t){
    # See Nicky J. Welton (2004)
    P<-diag(3)
    la1<-Q[1,2]+Q[1,3]
    la2<-Q[2,1]+Q[2,3]
    h<-sqrt((la1-la2)^2+4*Q[1,2]*Q[2,1])
    P[1,1]<-( (-la1+la2+h)*exp(-(la1+la2-h)*t/2)+(la1-la2+h)*exp(-(la1+la2+h)*t/2) )/(2*h)
    P[1,2]<-( (-la1+la2+h)*(la1-la2+h)*(exp(-(la1+la2-h)*t/2)-exp(-(la1+la2+h)*t/2)) )/(4*h*Q[2,1])
    P[1,3]<-1-P[1,1]-P[1,2]
    P[2,1]<-( Q[2,1]*(exp(-(la1+la2-h)*t/2)-exp(-(la1+la2+h)*t/2)) )/h
    P[2,2]<-((la1-la2+h)*exp(-(la1+la2-h)*t/2)+(-la1+la2+h)*exp(-(la1+la2+h)*t/2))/(2*h)
    P[2,3]<-1-P[2,1]-P[2,2]
    P[3,]<-c(0,0,1)
    P
}
# Build Q:
defineQ<-function(beta,beta12,age,educ,ybrth,educ.age,u.j,v.i){
      Q[1,2]<- exp(beta[1,1]+beta[1,2]*age+beta12[1]*educ+beta12[2]*ybrth+beta12[3]*educ.age+u.j)
      Q[1,3]<- exp(beta[2,1]+beta[2,2]*age+v.i);  Q[1,1]<- -sum(Q[1,])
      Q[2,1]<- exp(beta[3,1]+beta[3,2]*age);  
      Q[2,3]<- exp(beta[4,1]+beta[4,2]*age+v.i);  Q[2,2]<- -sum(Q[2,])
      return(Q)
}

# Prepare data:
if(0){
# Build data:
da <- new.env()
for(i in 1:N){
  data.i<-data[data$id==subjects[i],]
  j<-data.i$id.centre[1]
  eval(parse(text=paste(paste("da$d",as.character(i),".",as.character(j),sep=""),"<-data.i",sep="")))
}
da<-as.list(da)
cat("Data prepared.\n")
}

# Likelihood:
loglikelihood<-function(p){

 # Intensities parameters:
 beta<-cbind(p[1:4],p[5:8])
 sigma.j<-exp(p[9])
 sigma.i<-exp(p[10])
 beta12<-p[11:13]
 
 # Contribution per unit:
 loglik<-0
 for(j in 1:J){
  # Integrand for cluster effect u.j:
  integrand1<-function(u.j){
    loglik.j<-0
    for(i in which(idN.centre==j)){
    # Data for subject ij:
    eval(parse(text=paste("data.ij<-",paste("da$d",as.character(i),".",as.character(j),sep=""),sep="")))
    O<-data.ij$state; L.ij<-length(O)
    t<-data.ij$time
    t.age<-data.ij$age-60
    educ<-data.ij$educ[1]
    t.ybrth<-data.ij$ybrth[1]-1900
    educ.t.age<-educ*t.age
    tprobs<-rep(NA,L.ij-1)
    
    # Integrand for effect v.i:
    integrand2<-function(v.i){
     # Likelihood contribution:
     if(L.ij>2){
      for(k in 2:(L.ij-1)){
       # Q and P matrix:
       Q.k<-defineQ(beta,beta12,t.age[k-1],educ,t.ybrth,educ.t.age[k-1],u.j,v.i)
       P<-matrixexp(Q=Q.k,t=t[k]-t[k-1])
       # No censoring, no death:
       tprobs[k-1]<- P[O[k-1],O[k]] 
       }
     }
     # Q matrix:
     k <- L.ij
     Q.k<-defineQ(beta,beta12,t.age[k-1],educ,t.ybrth,educ.t.age[k-1],u.j,v.i)
     if(O[L.ij]==dead){
       # P matrices:
       P<-matrixexp(Q=Q.k,t=t[k]-eps-t[k-1])
       P.eps<-matrixexp(Q=Q.k,t=eps)
       # Death:
       tprobs[k-1]<-P[O[k-1],1]*P.eps[1,dead]+P[O[k-1],2]*P.eps[2,dead]
     }else{
       # P matrix:
       P<-matrixexp(Q=Q.k,t=t[k]-t[k-1])
       # Censoring:
       tprobs[k-1]<- P[O[k-1],1]+P[O[k-1],2]
     }
     return(prod(tprobs)*dnorm(v.i,mean=0,sd=sigma.i))
   }
   int2<-integraal(integrand2,mu=0,range=ww*sigma.i)
   loglik.j<-loglik.j+log(int2)
   }
   return(exp(loglik.j+BigNumber+log(dnorm(u.j,mean=0,sd=sigma.j))))
 }
 int1<-integraal(integrand1,mu=0,range=ww*sigma.j)
 loglik<-loglik+log(int1)-BigNumber
 }
 # Monitoring:
 cat("-2*Loglik = ", -2*loglik,"\n")
 cat("Parameters: ", round(c(beta,beta12,sigma.j,sigma.i),2),"\n\n")
 -loglik
}

# Starting values from fixed-effects model estimated in msm():
#p<-c(p.msm[1:8],-.5,-.5,p.msm[9:11])

# Starting values from BUGS:
#p<-c(-2.195,-3.992,	-2.104,	-1.727, 0.0439,0.06273,-0.02665,0.02428,log(0.2567),log(0.06767),-0.3663,-0.0698,0.008621)
# Starting values from fixed-effect model:
#p<-c( -2.61,-3.99,-2.14,-1.69,0.05,0.06,-0.02,0.02,log(1),log(1), -0.33,-0.05,0.008)
p<-max$par

# Maximise:
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 =", method,"with",nnodes,"nodes.\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<-10
sigma.j <-exp(p[L-1])
sigma.i <-exp(p[L])
# Delta method for sigmas:
sigma.j.se<- p.se[L-1]*exp(p[L-1])
sigma.i.se<- p.se[L]*exp(p[L])
cat("sigma.j = ",round(sigma.j,digits),"with SE",round(sigma.j.se,digits),"\n")
cat("sigma.i = ",round(sigma.i,digits),"with SE",round(sigma.i.se,digits),"\n\n")


# MAP estimation of frailties:
beta<-cbind(p[1:4],p[5:8])
sigma.j<-exp(p[9])
sigma.i<-exp(p[10])
beta12<-p[11:13]
cat("MAP estimation of centre frailties...\n")
frailty<-rep(NA,J)
frailty.se<-rep(NA,J)
info<-rep(NA,J)
for(j in 1:J){
cat("Centre",j,"...\n")
index<- which(idN.centre==j)
loglikelihood.j<-function(p){
   u.j<-p
   loglik.j<-0
   for(i in index){
    # Data for subject ij:
    eval(parse(text=paste("data.ij<-",paste("da$d",as.character(i),".",as.character(j),sep=""),sep="")))
    O<-data.ij$state; L.ij<-length(O)
    t<-data.ij$time
    t.age<-data.ij$age-60
    educ<-data.ij$educ[1]
    t.ybrth<-data.ij$ybrth[1]-1900
    educ.t.age<-educ*t.age
    tprobs<-rep(NA,L.ij-1)
    
    # Integrand for effect v.i:
    integrand2<-function(v.i){
     # Likelihood contribution:
     if(L.ij>2){
      for(k in 2:(L.ij-1)){
       # Q and P matrix:
       Q.k<-defineQ(beta,beta12,t.age[k-1],educ,t.ybrth,educ.t.age[k-1],u.j,v.i)
       P<-matrixexp(Q=Q.k,t=t[k]-t[k-1])
       # No censoring, no death:
       tprobs[k-1]<- P[O[k-1],O[k]] 
       }
     }
     # Q matrix:
     k <- L.ij
     Q.k<-defineQ(beta,beta12,t.age[k-1],educ,t.ybrth,educ.t.age[k-1],u.j,v.i)
     if(O[L.ij]==dead){
       # P matrices:
       P<-matrixexp(Q=Q.k,t=t[k]-eps-t[k-1])
       P.eps<-matrixexp(Q=Q.k,t=eps)
       # Death:
       tprobs[k-1]<-P[O[k-1],1]*P.eps[1,dead]+P[O[k-1],2]*P.eps[2,dead]
     }else{
       # P matrix:
       P<-matrixexp(Q=Q.k,t=t[k]-t[k-1])
       # Censoring:
       tprobs[k-1]<- P[O[k-1],1]+P[O[k-1],2]
     }
     return(prod(tprobs)*dnorm(v.i,mean=0,sd=sigma.i))
   }
   loglik.j<-loglik.j+log(integraal(integrand2,mu=0,range=3*sigma.i))
  }
 loglik.j<-loglik.j+log(dnorm(u.j,mean=0,sd=sigma.j))
 # Monitoring:
 #cat("Loglik.j = ", loglik.j,"\n")
 return(-loglik.j)
}
start.par<-0
max.j<-optim(par=start.par, fn=loglikelihood.j, method = "BFGS",control=list(trace=0,maxit=20000),hessian=TRUE)
frailty[j]<-max.j$par
frailty.se[j]<-sqrt(1/max.j$hessian)
info[j]<-max.j$convergence
}
tab<-frailty
names(tab)<-centres
print(round(tab,digits))
cat("With SEs:",round(frailty.se,digits),"\n")