LRT test on Gompertz-Makeham model,

# I used single-fitting as initial values for the double-fitting. It worked. 
# Hessian does seems to be available when lower bounds are specifized in optim()



#20140818 Nested model test on LS1 RLS using the Gompertz-Makehma model
# Conclusions:
# H0: all I, G, M are the same between control and LS1
# H1: M are different between control and LS1
# p = 1.285958e-05, indicating that H1 is significantly different from the hull hypothesis

# H0 I = 1.000000e-05, G= 0.35, M = 0.021
# H1: I = 5.597089e-07, G = 0.45, M1 = 0.0063, M2 = 0.052


rm(list=ls());
source("~/lib/R/lifespan.r")
require(nlme)
require(survival)
require(xlsx)
require(flexsurv)

tb = read.xlsx( "LS1 MD data.xlsx", 1 );
ctl = tb[,2]
ls1 = tb[,3]

fitCtlGom = flexsurvreg(formula = Surv(ctl) ~ 1, dist="gompertz")
fitCtlGom$res

S = calculate.s(ctl)
s= S$s; t=S$t;
g1 = gnls( s ~  exp( (I/G) *(1 - exp(G* t)) - M*t ), start=list( I=0.003, G=0.16, M=0.01) )
retGM = optim ( g1$coef,  llh.GM.single.run, lifespan=ctl,
               lower = c(1E-10, 1E-5, 0), upper = c(0.2, 2, 0.1), method="L-BFGS-B");
retG = optim ( c(3.6E-5, 0.3),  llh.G.single.run, lifespan=ctl,
               lower = c(1E-10, 1E-5), upper = c(0.2, 2 ), method="L-BFGS-B");

fitLS1Gom = flexsurvreg(formula = Surv(ls1) ~ 1, dist="gompertz")
fitLS1Gom$res

S2 = calculate.s(ls1)
s = S2$s; t = S2$t
g2 = gnls( s ~  exp( (I/G) *(1 - exp(G* t)) - M*t ), start=list( I=0.003, G=0.16, M=0.01) )
retGM2 = optim ( g2$coef,  llh.GM.single.run, lifespan=ls1, lower = c(1E-10, 1E-5, 0), upper = c(0.2, 2, 0.1), method="L-BFGS-B");
retG2 = optim ( c(0.0166, 0.1),  llh.G.single.run, lifespan=ls1,
               lower = c(1E-10, 1E-5), upper = c(0.2, 2 ), method="L-BFGS-B");

retGM$par
retGM2$par
initIGM = c(retGM$par, retGM2$par)

#################################
##### LRT to exam whether two data on I,G,M.
#                            rawIGM =c( I1,     G1,          M1,      I2,         G2,     M2 )
H0   <- function( rawIGM ) { IG <- c(rawIGM[1], rawIGM[2], rawIGM[3], rawIGM[1],  rawIGM[2], rawIGM[3]) }  #all the same
H1m  <- function( rawIGM ) { IG <- c(rawIGM[1], rawIGM[2], rawIGM[3], rawIGM[1],  rawIGM[2], rawIGM[6]) } # M different

Hx.llh.GM2sample.single.run <- function( rawIGM, model, lifespan1, lifespan2 ) {
  IGM = model(rawIGM);
  I1 = IGM[1]; G1 = IGM[2]; M1=IGM[3]; I2 = IGM[4]; G2 = IGM[5]; M2=IGM[6];
  my.lh1 = llh.GM.single.run(IGM[1:3], lifespan1)
  my.lh2 = llh.GM.single.run(IGM[4:6], lifespan2)
  my.lh = my.lh1 + my.lh2
  print (IGM ); #trace the convergence
  ret = my.lh
}

rawIGM = c(retGM$par, retGM2$par)
llh.H0   = optim( initIGM, Hx.llh.GM2sample.single.run, model=H0,   lifespan1=ctl, lifespan2=ls1,
                  method="L-BFGS-B", lower=c(1E-5,1E-5,1E-5, 1E-5,1E-5,1E-5) );
                 
llh.H1m  = optim( initIGM, Hx.llh.GM2sample.single.run, model=H1m,  lifespan1=ctl, lifespan2=ls1,
                  method="L-BFGS-B", lower=c(1E-7,1E-7,1E-7, 1E-7,1E-7,1E-7) );

rbind( llh.H0$par, llh.H1m$par)
deltaLH = llh.H0$value - rbind( llh.H0$value, llh.H1m$value)

1 - pchisq( 2*deltaLH, df =1 );


#################################
###### End of LRT    ############
#################################