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    ############
#################################

gradient decent for optimization

Notes from Andrew Ng's Machine learning, Coursera. 

Compare binomial, gompertz, and weibull model, fitting with simulated binomial-aging-model lifespan

# ~/0.network.aging.prj/0a.network.fitting/_binomial_aging20140711.r

#test binomial aging model
rm(list=ls());

source("lifespan.20131221.r")
require(flexsurv)


##### log likelihood function, 3-parameter binomial mortality rate model
# http://hongqinlab.blogspot.com/2013/12/binomial-mortailty-model.html
# m = R ( 1 + t/t0)^(n-1)
# s = exp( (R t0/n)*(1 - (1+t/t0)^n ) )  

llh.binomialMortality.single.run <- function( Rt0n, lifespan ) {
  I = Rt0n[1]; t0 = Rt0n[2]; n=Rt0n[3];
  my.data = lifespan[!is.na(lifespan)];
  log_s = (I * t0 /n )*(1 - (1 + my.data/t0)^n);
  log_m = log(I) +  (n-1) * log(1 + my.data/t0 ); 
  my.lh = sum(log_s)  + sum(log_m);
  print (Rt0n ); #trace the convergence
  ret = - my.lh # because optim seems to maximize 
}

#random number from binomial-aging model
rbinomial.aging <- function ( Rt0n, n){
    x.uniform = runif(n)
    inverse.binomialaging.CDF = function(R,t0,n,y) {t0*((1- n*log(y)/(R*t0))^(1/n) -1) }
    x.binomialaging = inverse.binomialaging.CDF(Rt0n[1],Rt0n[2], Rt0n[3], x.uniform)
    return(x.binomialaging)
}

set.seed(20140711)
x = rbinomial.aging( c(0.001,10, 5), 1000)
fitBinom = optim ( c(0.001,10, 3),  llh.binomialMortality.single.run , lifespan=x, 
              method="L-BFGS-B", lower=c(1E-10,0.05,1E-1), upper=c(10,100,50) );  

fitGomp = optim ( c(0.005, 0.08),  llh.G.single.run , lifespan=x,
                  method="L-BFGS-B", lower=c(1E-10,1E-2), upper=c(10,10) );  

fitWeib2 = flexsurvreg(formula = Surv(x) ~ 1, dist="weibull")
fitGomp2 = flexsurvreg(formula = Surv(x) ~ 1, dist="gompertz")
fitGomp2$coeff

c(-fitGomp$value, fitGomp2$loglik, fitWeib2$loglik, -fitBinom$value)

Flexsurv can over-fit the Gompertz with negative values, and Weibull gave the best value!!! 
Even though the random number are generated by binomial aging model, its fitting gave the lowest likelihood?!

TODO: I should to provide a gradient function for optim(method='L-BFGS-B')
Or try optimx()

For binomial-aging, $n$ is under-estimated by optim(), but $m0$ is over-estimated. When the entire data is fitted with Gompertz model, $G$ is too low. How about fitting within initial stage? 

set.seed(20140711)

x = round( rbinomial.aging( c(0.001,10, 4), 5000), digits=1 )#t0=10, n=4

tb = calculate.mortality.rate(x)

plot( tb$s ~ tb$t)

plot( tb$mortality.rate ~ tb$t )

plot( tb$mortality.rate ~ tb$t, log='y', main='linear for Gompertz' )

plot( tb$mortality.rate ~ tb$t, log='yx', main='linear for Weibull' )

set.seed(20140711)
x = rbinomial.aging( c(0.01,25, 4), 100)  #large m0, t0=25, n=4
x = round( x + 0.5, digits=0)
fitBinom = optim ( c(0.01,25, 3.5),  llh.binomialMortality.single.run , lifespan=x,
              method="L-BFGS-B", lower=c(1E-10,0.05,1E-1), upper=c(10,100,50) );

fitGomp = optim ( c(0.005, 0.08),  llh.G.single.run , lifespan=x,
                  method="L-BFGS-B", lower=c(1E-10,1E-2), upper=c(10,10) );

fitWeib2 = flexsurvreg(formula = Surv(x) ~ 1, dist="weibull")
fitGomp2 = flexsurvreg(formula = Surv(x) ~ 1, dist="gompertz")
fitGomp2$coeff

c(-fitGomp$value, fitGomp2$loglik, fitWeib2$loglik, -fitBinom$value)

Good, binomial aging model gives the best likelihood.

Plot confirmed that Gompertz looks better than Weibull. 

Align two vector, first try, not perfect results

See "align_curves.20140316.R" at github/project.H2O2.tolerance/Feb14,2011.H2O2.titration.bioscreen.star.strain

#20140316 align growth curves by minimizing sum of errors 

rm(list=ls())
debug = 9; 

#
pairwise_sum_of_errors = function( X1, X2, Start1,End1,Start2, End2 ) {
  if( (End1 - Start1) == (End2 - Start2) ) {
    return( sum( ( X1[Start1:End1]-X2[Start2:End2] )^2))
  } else {
    print("Error: X1 and X2 should have the same lengths")
    return( NA); 
  }
}

#testing 
x1 = 1:100;  x2 = 5:104
pairwise_sum_of_errors( x1, x2, 5, 100, 1,96)
sum((x1[5:100]-x2[1:96])^2)
pairwise_sum_of_errors( x1, x2, 5, 100, 1,100)

#align two vectors with equal lengths
#Let's assume start2 lags behind start1 
#start2 = start1 + delta # 2nd start lags behind 1st start by delta
pairwise_cost_of_shifted_vectors = function( delta, X1, X2,  low.threshold=0.01, debug=10 ){
  delta = floor( (delta+0.5) ) #delta must be an integer
  X1[X1<low.threshold]=0;
  X2[X2<low.threshold]=0;
  Start1 = 1;
  Start2 = Start1 + delta;  #left truncation at delta
  End2 = length(X2)
  End1 = length(X1) - delta; #right truncation at delta
  if(debug){ print(paste('Start1', Start1, "Start2", Start2, "End1", End1, "End2", End2))
  }
  if( (End1 - Start1) == (End2 - Start2) ) {
    return( sum( ( X1[Start1:End1]-X2[Start2:End2] )^2))
  } else {
    print("Error: X1 and X2 should have the same lengths")
    return( NA); 
  }
}

#testing 1
x1 = 5:104; x2 = 1:100;  
pairwise_sum_of_errors( x1, x2, 1,96, 5,100 )
pairwise_cost_of_shifted_vectors( 4, x1, x2)

delta = 4
res1 = optim(c(delta), fn=pairwise_cost_of_shifted_vectors, X1=x1, X2=x2, method="Brent", lower=c(1), upper=c(100))
#return 3 not 4? 

#testing 2
x1 = 50:149; x2 = 1:100;  
pairwise_sum_of_errors( x1, x2, 1,50, 50,99 )
pairwise_cost_of_shifted_vectors( 49, x1, x2)

delta = 50
res2 = optim(c(delta), fn=pairwise_cost_of_shifted_vectors, X1=x1, X2=x2, method="Brent", lower=c(1), upper=c(100))
#return 50 not 49?

Verify optim() estimation using 2 and 3 parameter Gompertz model.

Using 2 and 3 parameter Gompertz model. When M is large, ignoring M will inflate I, the initial mortality rate. 




rm(list=ls())
source("/Users/hongqin/lib/R/lifespan.r")

##two parameter Gompertz model 
N=1000
I =0.005
G = 0.10
t= seq(1, 100, by=1)
s = G.s(c(I,G,0), t)
plot( s ~ t)
mu = I * exp(G*t)
plot ( s*mu ~ t )
pmf = s * mu  # prob mass function
pmf = pmf / sum(pmf)
summary ( pmf * N)
hist( round( pmf * N) )
lifespanT = rep( t, round(pmf*N))
hist(lifespanT, br=20)
summary(lifespanT)


##### log likelihood function, simple gompertz mortality model
#s = exp( (I/G) *(1 - exp(G* my.data)) )  ;
#m = I * exp( G * my.data );   
llh.gompertz.single.run <- function( IG, lifespan, trace=0 ) {
  #print(IG)
  I = IG[1]; G = IG[2]; 
  my.data = lifespan[!is.na(lifespan)];
  log_s = (I/G) *(1 - exp(G* my.data))
  #if( I< 0 ) { I = 1E-10 }
  log_m = log(I) +  G * my.data ; 
  my.lh = sum(log_s)  + sum(log_m);
  if (trace) {
    print (IG ); #trace the convergence
  }
  ret = - my.lh # because optim seems to maximize 
}

lifespan = sample(lifespanT, 30)
shat = calculate.s( lifespan )
plot(shat$s ~ shat$t )

ret1a = optim ( c(I,G)*2, fn=llh.gompertz.single.run, lifespan=lifespan, lower=c(1E-10, 1E-5), method="L-BFGS-B" );
ret1a$par / c(I,G)

ret1a = optim ( c(I, G), fn=llh.gompertz.single.run, lifespan=lifespan, lower=c(1E-10, 1E-5), method="L-BFGS-B" );
ret1a$par / c(I,G)

ret1b = optim ( c(I, G)*c(0.01, 1.5), fn=llh.gompertz.single.run, lifespan=lifespan, lower=c(1E-10, 1E-5), method="L-BFGS-B" );
ret1b$par / c(I,G)

ret1c = optim ( c(I,G)*c(0.01,1.5), fn=llh.gompertz.single.run, lifespan=lifespan, trace=1  );
ret1c$par / c(I, G)

ret1d = optim ( c(0.05, 0.15), fn=llh.gompertz.single.run, lifespan=lifespan, trace=1  );
ret1d$par / c(I,G)

ret1e = optim ( c(0.05, 0.1), fn=llh.gompertz.single.run, lifespan=lifespan  );
ret1e$par / ret1d$par

#################
## 3  parameter Gompertz model 
N=1000
I =0.005
G = 0.10
M = 0.05

t= seq(1, 100, by=1)
s = GM.s(c(I,G,0), t)
plot( s ~ t)

mu = I * exp(G*t) + M

plot ( s*mu ~ t )
pmf = s * mu  # prob mass function
pmf = pmf / sum(pmf)
summary ( pmf * N)
hist( round( pmf * N) )
lifespanT = rep( t, round(pmf*N))
hist(lifespanT, br=20)
summary(lifespanT)



lifespan = sample(lifespanT, 30)
shat = calculate.s( lifespan )
plot(shat$s ~ shat$t )

ret1a = optim ( c(I,G, M)*2, fn=llh.GM.single.run, lifespan=lifespan, lower=c(1E-10, 1E-5), method="L-BFGS-B" );
ret1a$par / c(I,G, M)


ret1b = optim ( c(I, G), fn=llh.G.single.run, lifespan=lifespan, lower=c(1E-10, 1E-5), method="L-BFGS-B" );
ret1b$par / c(I,G)
#ignore larger M can inflate R



ret1b = optim ( c(I, G)*c(0.01, 1.5), fn=llh.gompertz.single.run, lifespan=lifespan, lower=c(1E-10, 1E-5), method="L-BFGS-B" );
ret1b$par / c(I,G)

ret1c = optim ( c(I,G)*c(0.01,1.5), fn=llh.gompertz.single.run, lifespan=lifespan, trace=1  );
ret1c$par / c(I, G)

ret1d = optim ( c(0.05, 0.15), fn=llh.gompertz.single.run, lifespan=lifespan, trace=1  );
ret1d$par / c(I,G)

ret1e = optim ( c(0.05, 0.1), fn=llh.gompertz.single.run, lifespan=lifespan  );
ret1e$par / ret1d$par

Investigate optim() by simulation of two parameter Gompertz model

To see how optim() works in R, I tested it using simulated data. Method L-BFGS-B appears to give the most accurate fitting result. For methods tested, the initial values seem to really influence the accuracy of the fitting result.



rm(list=ls())
source("/Users/hongqin/lib/R/lifespan.r")
N=1000
I =0.005
G = 0.10
t= seq(1, 100, by=1)
s = G.s(c(I,G,0), t)
plot( s ~ t)
mu = I * exp(G*t)
plot ( s*mu ~ t )
pmf = s * mu  # prob mass function
pmf = pmf / sum(pmf)
summary ( pmf * N)
hist( round( pmf * N) )
lifespanT = rep( t, round(pmf*N))
hist(lifespanT, br=20)
summary(lifespanT)


##### log likelihood function, simple gompertz mortality model
#s = exp( (I/G) *(1 - exp(G* my.data)) )  ;
#m = I * exp( G * my.data );  
llh.gompertz.single.run <- function( IG, lifespan, trace=0 ) {
  #print(IG)
  I = IG[1]; G = IG[2];
  my.data = lifespan[!is.na(lifespan)];
  log_s = (I/G) *(1 - exp(G* my.data))
  #if( I< 0 ) { I = 1E-10 }
  log_m = log(I) +  G * my.data ;
  my.lh = sum(log_s)  + sum(log_m);
  if (trace) {
    print (IG ); #trace the convergence
  }
  ret = - my.lh # because optim seems to maximize
}

lifespan = sample(lifespanT, 30)
shat = calculate.s( lifespan )
plot(shat$s ~ shat$t )

ret1a = optim ( c(I,G)*2, fn=llh.gompertz.single.run, lifespan=lifespan, lower=c(1E-10, 1E-5), method="L-BFGS-B" );
ret1a$par / c(I,G)

ret1a = optim ( c(I, G), fn=llh.gompertz.single.run, lifespan=lifespan, lower=c(1E-10, 1E-5), method="L-BFGS-B" );
ret1a$par / c(I,G)

ret1b = optim ( c(I, G)*c(0.01, 1.5), fn=llh.gompertz.single.run, lifespan=lifespan, lower=c(1E-10, 1E-5), method="L-BFGS-B" );
ret1b$par / c(I,G)

ret1c = optim ( c(I,G)*c(0.01,1.5), fn=llh.gompertz.single.run, lifespan=lifespan, trace=1  );
ret1c$par / c(I, G)

ret1d = optim ( c(0.05, 0.15), fn=llh.gompertz.single.run, lifespan=lifespan, trace=1  );
ret1d$par / c(I,G)

ret1e = optim ( c(0.05, 0.1), fn=llh.gompertz.single.run, lifespan=lifespan  );
ret1e$par / ret1d$par

R code for nested model test, likelihood ratio test, using AgingData 0517.xls

2013 Nov 6 Note: Permutation based test are probably more 'accurate' at estimating p-values than Chi-square. 


rm( list=ls())
list.files()

file = "AgingData_0517.csv"
tb = read.csv( file );


##### log likelihood function, simple gompertz mortality model
   #s = exp( (I/G) *(1 - exp(G* my.data)) )  ;
   #m = I * exp( G * my.data );  
 llh.gompertz.single.run <- function( IG, lifespan ) {
   #print(IG)
   I = IG[1]; G = IG[2];
   my.data = lifespan[!is.na(lifespan)];
   log_s = (I/G) *(1 - exp(G* my.data))
   if( I< 0 ) { I = 1E-10 }
   log_m = log(I) +  G * my.data ;
   my.lh = sum(log_s)  + sum(log_m);
   print (IG ); #trace the convergence
   ret = - my.lh # because optim seems to maximize
 }

#####

ret1a = optim ( c(0.0034, 0.1), fn=llh.gompertz.single.run, lifespan=tb[,1], lower=c(1E-10, 1E-5), method="L-BFGS-B" );
ret1b = optim ( c(0.05, 0.2), fn=llh.gompertz.single.run, lifespan=tb[,1], lower=c(1E-10, 1E-5), method="L-BFGS-B" );
ret1c = optim ( c(0.039, 0.1), fn=llh.gompertz.single.run, lifespan=tb[,1]  );
ret1d = optim ( c(0.05, 0.15), fn=llh.gompertz.single.run, lifespan=tb[,1]  );
ret1c$par / ret1d$par
ret1e = optim ( c(0.1, 0.1), fn=llh.gompertz.single.run, lifespan=tb[,1]  );
ret1e$par / ret1d$par

#set up a template for fitting parameters
label = c("R","G", "mean", "median", "stddev", "CV", 'p_wilcox', 'p_ks')
allFitPar = data.frame( matrix(NA, nrow=length(label), ncol=length(tb[1,])) )
rownames(allFitPar) = label
names(allFitPar) = names(tb)

 
#loop over
for (col in names(tb)) {
  cur_lifespan = tb[,col]
  retTmp = optim ( c(0.039, 0.1), fn=llh.gompertz.single.run, lifespan=cur_lifespan  );
  allFitPar[1:2, col] = retTmp$par
  allFitPar["mean", col] = mean(cur_lifespan, na.rm=T)
  allFitPar["median", col] = median(cur_lifespan, na.rm=T)
  allFitPar["stddev", col] = sqrt( var(cur_lifespan, na.rm=T))
  allFitPar["CV", col] = allFitPar["stddev", col] / allFitPar["mean", col]
  tmp = wilcox.test( tb[,'wt'], tb[,col] ) 
  allFitPar['p_wilcox', col] = tmp$p.value
  tmp2 = ks.test( tb[,'wt'], tb[,col] ) 
  allFitPar['p_ks', col] = tmp2$p.value
 
}
allFitPar


#################################
##### LRT
### 1) model H0, same G and same I
### 2) model H1i, same G, different I
### 3) model H1g, different G, same I
### 4) model H2, different G, different I
                                  #    I1       G1       I2         G2
 H0  <- function( rawIG ) { IG <- c(rawIG[1], rawIG[2], rawIG[1], rawIG[2] ) }  #all the same
 H1i <- function( rawIG ) { IG <- c(rawIG[1], rawIG[2], rawIG[3], rawIG[2] ) }  #different I
 H1g <- function( rawIG ) { IG <- c(rawIG[1], rawIG[2], rawIG[1], rawIG[4] ) }  # different G
 H2  <- function( rawIG ) { IG <- c(rawIG[1], rawIG[2], rawIG[3], rawIG[4] ) }  # all different

 Hx.llh.gompertz.single.run <- function( rawIG, model, lifespan1, lifespan2 ) {
   IG = model(rawIG);
   I1 = IG[1]; G1 = IG[2]; I2 = IG[3]; G2 = IG[4];
   my.data1 = lifespan1[!is.na(lifespan1)];
   my.data2 = lifespan2[!is.na(lifespan2)];
   log_s1 = (I1/G1) *(1 - exp(G1* my.data1))
   log_s2 = (I2/G2) *(1 - exp(G2* my.data2))
   log_m1 = log(I1) +  G1 * my.data1 ;
   log_m2 = log(I2) +  G2 * my.data2 ;
   my.lh = sum(log_s1) + sum(log_m1) + sum(log_s2) + sum(log_m2)
   print (IG ); #trace the convergence
   ret = - my.lh # because optim seems to maximize
 }

## LRT to exam whether wt, wtCR share the same G, I
llh.H0   = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H0,   lifespan1=tb$wt, lifespan2=tb$wtCR );
llh.H1i  = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H1i,  lifespan1=tb$wt, lifespan2=tb$wtCR );
llh.H1g  = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H1g,  lifespan1=tb$wt, lifespan2=tb$wtCR );
llh.H2  = optim( c(0.01,0.2,0.01,0.1),   Hx.llh.gompertz.single.run, model=H2,  lifespan1=tb$wt, lifespan2=tb$wtCR );

cbind(llh.H0$par, llh.H1i$par, llh.H1g$par, llh.H2$par);

LH = c(llh.H0$value, llh.H1i$value, llh.H1g$value, llh.H2$value);
LH
deltaLH =  - LH + llh.H0$value;
deltaLH
1 - pchisq( 2*deltaLH, df =c(1,1,1,2) );
#p=0.06 not significant for all models.

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

## LRT to exam whether ybr053cD, and CR share the same G, I
llh.H0   = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H0,   lifespan1=tb[,3], lifespan2=tb[,4] );
llh.H1i  = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H1i,  lifespan1=tb[,3], lifespan2=tb[,4] );
llh.H1g  = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H1g,  lifespan1=tb[,3], lifespan2=tb[,4] );
llh.H2   = optim( c(0.01,0.2,0.01,0.1),   Hx.llh.gompertz.single.run, model=H2,  lifespan1=tb[,3], lifespan2=tb[,4] );

cbind(llh.H0$par, llh.H1i$par, llh.H1g$par, llh.H2$par);

LH = c(llh.H0$value, llh.H1i$value, llh.H1g$value, llh.H2$value);
LH
deltaLH =  - LH + llh.H0$value;
deltaLH
1 - pchisq( 2*deltaLH, df =c(1,1,1,2) );
#[1] 1.00000000 0.01498437 0.01368829 0.04515191

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

## LRT to exam whether ybr053cD, and CR share the same G, I
llh.H0   = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H0,   lifespan1=tb[,5], lifespan2=tb[,6] );
llh.H1i  = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H1i,  lifespan1=tb[,5], lifespan2=tb[,6] );
llh.H1g  = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H1g,  lifespan1=tb[,5], lifespan2=tb[,6] );
llh.H2   = optim( c(0.01,0.2,0.01,0.1),   Hx.llh.gompertz.single.run, model=H2,  lifespan1=tb[,5], lifespan2=tb[,6] );
cbind(llh.H0$par, llh.H1i$par, llh.H1g$par, llh.H2$par);

LH = c(llh.H0$value, llh.H1i$value, llh.H1g$value, llh.H2$value);
LH
deltaLH =  - LH + llh.H0$value;
deltaLH
1 - pchisq( 2*deltaLH, df =c(1,1,1,2) );
#[1] 1.00000000 0.04595566 0.31817030 0.04188332

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


## LRT to exam whether ybr053cD, and CR share the same G, I
llh.H0   = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H0,   lifespan1=tb[,7], lifespan2=tb[,8] );
llh.H1i  = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H1i,  lifespan1=tb[,7], lifespan2=tb[,8] );
llh.H1g  = optim( c(0.01,0.2,0.01,0.1),  Hx.llh.gompertz.single.run, model=H1g,  lifespan1=tb[,7], lifespan2=tb[,8] );
llh.H2   = optim( c(0.01,0.2,0.01,0.1),   Hx.llh.gompertz.single.run, model=H2,  lifespan1=tb[,7], lifespan2=tb[,8] );

cbind(llh.H0$par, llh.H1i$par, llh.H1g$par, llh.H2$par);

LH = c(llh.H0$value, llh.H1i$value, llh.H1g$value, llh.H2$value);
LH
deltaLH =  - LH + llh.H0$value;
deltaLH
1 - pchisq( 2*deltaLH, df =c(1,1,1,2) );
#[1] 1.00000000 0.04595566 0.31817030 0.04188332

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

optim() usage in R.

optim(par, fn, gr = NULL, ...,
      method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"),
      lower = -Inf, upper = Inf,
      control = list(), hessian = FALSE) 

'par' is the default parameter for 'fn'.
Additional parameters for 'fn' can be passed on by specifying them after 'fn', such as
optim( par=c(0.05, 0.1), fn=llh, lifespan=current_lifespan)

Reference:
http://www.biostat.umn.edu/~ph7440/pubh7440/OptimLinearModel.R