nested models

Cite Ziheng Yang's codeML, nested model test approaches

https://stats.stackexchange.com/questions/4717/what-is-the-difference-between-a-nested-and-a-non-nested-model

• Cox, D. R.: Tests of separate families of hypotheses, Proceedings 4th Berkeley Symposium in Mathematical Statistics and Probability, 1, 105–123 (1961). University of California Press.
• Cox, D. R.: Further results on test of separate families of hypotheses. Journal of the Royal Statistical Society B, 406–424 (1962).

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

sir2, SIR2 overexpression, LRT

Conclusions: 

(1) sir2D and WT have different R and G. 

(2) sir2D and SIR2 overexpression share the same G. 

Black: by4742; Red: sir2D, Blue: SIR2-overexpression.

nested binomial-aging model test,

Summary: This is proof-of-concept work, using simulated data to see whether nested-model test can discern change of parameters. 

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

source("lifespan.20140711.r")

#################################
##### LRT to exam two datasets on binomial aging model, 
### model H0,  R, t0, n are all the same
### model H1i  R1<>R2
### model H1g, 
### model H2, 
#  rawIt0n = c( R1, t01, n1, R2, t02, n2)
H0  <- function( rawIt0n )  { Rt0n <- c(rawIt0n[1],rawIt0n[2],rawIt0n[3],rawIt0n[1],rawIt0n[2],rawIt0n[3]) }  #all the same
H1r  <- function( rawIt0n ) { Rt0n <- c(rawIt0n[1],rawIt0n[2],rawIt0n[3],rawIt0n[4],rawIt0n[2],rawIt0n[3]) }  # R1<>R2
H1t  <- function( rawIt0n ) { Rt0n <- c(rawIt0n[1],rawIt0n[2],rawIt0n[3],rawIt0n[1],rawIt0n[5],rawIt0n[3]) }  # t01 <> t02
H1n  <- function( rawIt0n ) { Rt0n <- c(rawIt0n[1],rawIt0n[2],rawIt0n[3],rawIt0n[1],rawIt0n[2],rawIt0n[6]) }  # n1 <> n2

Hx.llh.binomialMortality.single.run <- function( INrawRt0n, model, lifespan1, lifespan2 ) {
  Rt0n = model( INrawRt0n )
  I1 = Rt0n[1]; t01 = Rt0n[2]; n1=Rt0n[3];
  I2 = Rt0n[4]; t02 = Rt0n[5]; n2=Rt0n[6];  
  x1 = lifespan1[!is.na(lifespan1)];
  x2 = lifespan2[!is.na(lifespan2)];

  log_s1 = (I1 * t01 /n1 )*(1 - (1 + x1/t01)^n1);
  log_m1 = log(I1) +  (n1-1) * log(1 + x1/t01 ); 
  log_s2 = (I2 * t02 /n2 )*(1 - (1 + x2/t02)^n2);
  log_m2 = log(I2) +  (n2-1) * log(1 + x2/t02 ); 
  my.lh = sum(log_s1)  + sum(log_m1) + sum(log_s2) + sum(log_m2)
  print (Rt0n ); #trace the convergence
  ret = - my.lh # because optim seems to maximize 
}

set.seed(20140711)
x1 = round( rbinomial.aging( c(0.01, 15, 4), 100) + 0.5, digits=0);
x2 = round( rbinomial.aging( c(0.01, 10, 4), 100) + 0.5, digits=0);

llh.H0  = optim( c(0.01, 15, 4, 0.01, 15, 4), Hx.llh.binomialMortality.single.run, model=H0,lifespan1=x1, lifespan2=x2, 
                 method="L-BFGS-B", lower=c(1E-5,0.05,0.1, 1E-5,0.05,0.1), upper=c(10,100,50,10,100,50) );  

llh.H1t  = optim( c(0.01, 15, 4, 0.01, 10, 4), Hx.llh.binomialMortality.single.run, model=H1t,lifespan1=x1, lifespan2=x2, 
                 method="L-BFGS-B", lower=c(1E-5,0.05,0.1, 1E-5,0.05,0.1), upper=c(10,100,50,10,100,50) );  

cbind(llh.H0$par, llh.H1t$par)
LH = c(llh.H0$value, llh.H1t$value)
deltaLH =  - LH + llh.H0$value; 
1 - pchisq( 2*deltaLH, df =1 );

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