Monday, July 11, 2016

binomial aging model fitting, natural isolates

20160711-> Fitting RLS dataset merged by strains.
20160706-0711. Finished batch fitting of all individual RLS data sets. The fitting results showed that ‘n’ are often in the range of [5,7], though it is very noisy.
rm(list=ls())
setwd("~/github/0.network.aging.prj.bmc/0a.rls.fitting")
library('flexsurv')
## Warning: package 'flexsurv' was built under R version 3.2.4
## Loading required package: survival
## Warning: package 'survival' was built under R version 3.2.5
source("lifespan.r")

parse the strains from files

files = list.files(path="../qinlab_rls/", pattern="rls.tab")
tmp1 = gsub("\\d{6}.", "", files)
redundant_strains = gsub(".rls.tab", "", tmp1)
strains = sort( unique( redundant_strains ))
strains
##  [1] "101S"           "BY4716"         "BY4741"         "BY4742"        
##  [5] "BY4743"         "JSBY4741"       "M1-2"           "M13"           
##  [9] "M14"            "M2-8"           "M22"            "M32"           
## [13] "M34"            "M5"             "M8"             "RM112N"        
## [17] "S288c"          "SGU57"          "sir2D.4741a"    "sir2D.4742"    
## [21] "sir2DSIR2.4742" "SK1"            "W303"           "YPS128"        
## [25] "YPS163"
report = data.frame(cbind(strains))
report$samplesize = NA; report$R=NA; report$t0=NA; report$n=NA; report$G=NA; 

Explore the fitting outcomes of ‘flexsurv’.

i=2
  tb = read.table( paste("../qinlab_rls/",files[i],sep=''), sep="\t")
  GompFlex = flexsurvreg(formula = Surv(tb[,1]) ~ 1, dist = 'gompertz')
  WeibFlex = flexsurvreg(formula = Surv(tb[,1]) ~ 1, dist = 'weibull')
str(GompFlex)
## List of 28
##  $ call          : language flexsurvreg(formula = Surv(tb[, 1]) ~ 1, dist = "gompertz")
##  $ dlist         :List of 6
##   ..$ name          : chr "gompertz"
##   ..$ pars          : chr [1:2] "shape" "rate"
##   ..$ location      : chr "rate"
##   ..$ transforms    :List of 2
##   .. ..$ :function (x)  
##   .. ..$ :function (x, base = exp(1))  
##   ..$ inv.transforms:List of 2
##   .. ..$ :function (x)  
##   .. ..$ :function (x)  
##   ..$ inits         :function (t, mf, mml, aux)  
##  $ aux           : NULL
##  $ ncovs         : int 0
##  $ ncoveffs      : int 0
##  $ mx            :List of 2
##   ..$ rate : int(0) 
##   ..$ shape: NULL
##  $ basepars      : int [1:2] 1 2
##  $ covpars       : NULL
##  $ AIC           : num 213
##  $ data          :List of 3
##   ..$ Y  : num [1:30, 1:6] 10 10 11 14 15 15 16 22 24 25 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:30] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:6] "time" "status" "start" "stop" ...
##   ..$ m  :'data.frame':  30 obs. of  2 variables:
##   .. ..$ Surv(tb[, 1]): Surv [1:30, 1:2] 10  10  11  14  15  15  16  22  24  25  ...
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : NULL
##   .. .. .. ..$ : chr [1:2] "time" "status"
##   .. .. ..- attr(*, "type")= chr "right"
##   .. ..$ (weights)    : num [1:30] 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 Surv(tb[, 1]) ~ 1
##   .. .. .. ..- attr(*, "variables")= language list(Surv(tb[, 1]))
##   .. .. .. ..- attr(*, "factors")= int(0) 
##   .. .. .. ..- attr(*, "term.labels")= chr(0) 
##   .. .. .. ..- attr(*, "order")= int(0) 
##   .. .. .. ..- attr(*, "intercept")= int 1
##   .. .. .. ..- attr(*, "response")= int 1
##   .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. .. ..- attr(*, "predvars")= language list(Surv(tb[, 1]))
##   .. .. .. ..- attr(*, "dataClasses")= Named chr "nmatrix.2"
##   .. .. .. .. ..- attr(*, "names")= chr "Surv(tb[, 1])"
##   .. ..- attr(*, "covnames")= chr(0) 
##   .. ..- attr(*, "covnames.orig")= chr(0) 
##   ..$ mml:List of 2
##   .. ..$ rate : num [1:30, 1] 1 1 1 1 1 1 1 1 1 1 ...
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:30] "1" "2" "3" "4" ...
##   .. .. .. ..$ : chr "(Intercept)"
##   .. .. ..- attr(*, "assign")= int 0
##   .. ..$ shape: NULL
##  $ datameans     : num(0) 
##  $ N             : int 30
##  $ events        : int 30
##  $ trisk         : num 810
##  $ concat.formula:Class 'formula' length 3 Surv(tb[, 1]) ~ 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "covnames")= chr(0) 
##   .. ..- attr(*, "covnames.orig")= chr(0) 
##  $ all.formulae  :List of 1
##   ..$ rate:Class 'formula' length 3 Surv(tb[, 1]) ~ 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##  $ dfns          :List of 8
##   ..$ p    :function (q, shape, rate = 1, lower.tail = TRUE, log.p = FALSE)  
##   ..$ d    :function (x, shape, rate = 1, log = FALSE)  
##   ..$ h    :function (x, shape, rate = 1, log = FALSE)  
##   ..$ H    :function (x, shape, rate = 1, log = FALSE)  
##   ..$ r    :function (n, shape = 1, rate = 1)  
##   ..$ DLd  :function (t, shape, rate)  
##   ..$ DLS  :function (t, shape, rate)  
##   ..$ deriv: logi TRUE
##  $ res           : num [1:2, 1:4] 0.146413 0.001603 0.098752 0.000421 0.194074 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "shape" "rate"
##   .. ..$ : chr [1:4] "est" "L95%" "U95%" "se"
##  $ res.t         : num [1:2, 1:4] 0.1464 -6.4361 0.0988 -7.7719 0.1941 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "shape" "rate"
##   .. ..$ : chr [1:4] "est" "L95%" "U95%" "se"
##  $ cov           : num [1:2, 1:2] 0.000591 -0.015967 -0.015967 0.464473
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "shape" "rate"
##   .. ..$ : chr [1:2] "shape" "rate"
##  $ coefficients  : Named num [1:2] 0.146 -6.436
##   ..- attr(*, "names")= chr [1:2] "shape" "rate"
##  $ npars         : int 2
##  $ fixedpars     : NULL
##  $ optpars       : int [1:2] 1 2
##  $ loglik        : num -104
##  $ logliki       : num [1:30] -5.01 -5.01 -4.87 -4.46 -4.33 ...
##  $ cl            : num 0.95
##  $ opt           :List of 6
##   ..$ par        : Named num [1:2] 0.146 -6.436
##   .. ..- attr(*, "names")= chr [1:2] "shape" "rate"
##   ..$ value      : num 104
##   ..$ counts     : Named int [1:2] 35 10
##   .. ..- attr(*, "names")= chr [1:2] "function" "gradient"
##   ..$ convergence: int 0
##   ..$ message    : NULL
##   ..$ hessian    : num [1:2, 1:2] 23564 810 810 30
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "shape" "rate"
##   .. .. ..$ : chr [1:2] "shape" "rate"
##  - attr(*, "class")= chr "flexsurvreg"
GompFlex$res
##               est         L95%        U95%          se
## shape 0.146412756 0.0987519623 0.194073550 0.024317178
## rate  0.001602638 0.0004214271 0.006094643 0.001092234
GompFlex$res.t
##              est        L95%       U95%         se
## shape  0.1464128  0.09875196  0.1940735 0.02431718
## rate  -6.4361044 -7.77186373 -5.1003451 0.68152238
GompFlex$opt$hessian
##            shape     rate
## shape 23564.3361 810.0597
## rate    810.0597  30.0000
str(WeibFlex)
## List of 28
##  $ call          : language flexsurvreg(formula = Surv(tb[, 1]) ~ 1, dist = "weibull")
##  $ dlist         :List of 6
##   ..$ name          : chr "weibull.quiet"
##   ..$ pars          : chr [1:2] "shape" "scale"
##   ..$ location      : chr "scale"
##   ..$ transforms    :List of 2
##   .. ..$ :function (x, base = exp(1))  
##   .. ..$ :function (x, base = exp(1))  
##   ..$ inv.transforms:List of 2
##   .. ..$ :function (x)  
##   .. ..$ :function (x)  
##   ..$ inits         :function (t, mf, mml, aux)  
##  $ aux           : NULL
##  $ ncovs         : int 0
##  $ ncoveffs      : int 0
##  $ mx            :List of 2
##   ..$ scale: int(0) 
##   ..$ shape: NULL
##  $ basepars      : int [1:2] 1 2
##  $ covpars       : NULL
##  $ AIC           : num 218
##  $ data          :List of 3
##   ..$ Y  : num [1:30, 1:6] 10 10 11 14 15 15 16 22 24 25 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:30] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:6] "time" "status" "start" "stop" ...
##   ..$ m  :'data.frame':  30 obs. of  2 variables:
##   .. ..$ Surv(tb[, 1]): Surv [1:30, 1:2] 10  10  11  14  15  15  16  22  24  25  ...
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : NULL
##   .. .. .. ..$ : chr [1:2] "time" "status"
##   .. .. ..- attr(*, "type")= chr "right"
##   .. ..$ (weights)    : num [1:30] 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 Surv(tb[, 1]) ~ 1
##   .. .. .. ..- attr(*, "variables")= language list(Surv(tb[, 1]))
##   .. .. .. ..- attr(*, "factors")= int(0) 
##   .. .. .. ..- attr(*, "term.labels")= chr(0) 
##   .. .. .. ..- attr(*, "order")= int(0) 
##   .. .. .. ..- attr(*, "intercept")= int 1
##   .. .. .. ..- attr(*, "response")= int 1
##   .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. .. ..- attr(*, "predvars")= language list(Surv(tb[, 1]))
##   .. .. .. ..- attr(*, "dataClasses")= Named chr "nmatrix.2"
##   .. .. .. .. ..- attr(*, "names")= chr "Surv(tb[, 1])"
##   .. ..- attr(*, "covnames")= chr(0) 
##   .. ..- attr(*, "covnames.orig")= chr(0) 
##   ..$ mml:List of 2
##   .. ..$ scale: num [1:30, 1] 1 1 1 1 1 1 1 1 1 1 ...
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:30] "1" "2" "3" "4" ...
##   .. .. .. ..$ : chr "(Intercept)"
##   .. .. ..- attr(*, "assign")= int 0
##   .. ..$ shape: NULL
##  $ datameans     : num(0) 
##  $ N             : int 30
##  $ events        : int 30
##  $ trisk         : num 810
##  $ concat.formula:Class 'formula' length 3 Surv(tb[, 1]) ~ 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "covnames")= chr(0) 
##   .. ..- attr(*, "covnames.orig")= chr(0) 
##  $ all.formulae  :List of 1
##   ..$ scale:Class 'formula' length 3 Surv(tb[, 1]) ~ 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##  $ dfns          :List of 8
##   ..$ p    :function (q, shape, scale = 1, lower.tail = TRUE, log.p = FALSE)  
##   ..$ d    :function (x, shape, scale = 1, log = FALSE)  
##   ..$ h    :function (x, ...)  
##   ..$ H    :function (x, ...)  
##   ..$ r    :function (n, shape, scale = 1)  
##   ..$ DLd  :function (t, shape, scale)  
##   ..$ DLS  :function (t, shape, scale)  
##   ..$ deriv: logi TRUE
##  $ res           : num [1:2, 1:4] 3.77 30 2.76 27.17 5.13 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "shape" "scale"
##   .. ..$ : chr [1:4] "est" "L95%" "U95%" "se"
##  $ res.t         : num [1:2, 1:4] 1.33 3.4 1.02 3.3 1.64 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "shape" "scale"
##   .. ..$ : chr [1:4] "est" "L95%" "U95%" "se"
##  $ cov           : num [1:2, 1:2] 0.02502 0.00225 0.00225 0.00255
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "shape" "scale"
##   .. ..$ : chr [1:2] "shape" "scale"
##  $ coefficients  : Named num [1:2] 1.33 3.4
##   ..- attr(*, "names")= chr [1:2] "shape" "scale"
##  $ npars         : int 2
##  $ fixedpars     : NULL
##  $ optpars       : int [1:2] 1 2
##  $ loglik        : num -107
##  $ logliki       : Named num [1:30] -5.13 -5.13 -4.87 -4.24 -4.06 ...
##   ..- attr(*, "names")= chr [1:30] "1" "2" "3" "4" ...
##  $ cl            : num 0.95
##  $ opt           :List of 6
##   ..$ par        : Named num [1:2] 1.33 3.4
##   .. ..- attr(*, "names")= chr [1:2] "shape" "scale"
##   ..$ value      : num 107
##   ..$ counts     : Named int [1:2] 2 1
##   .. ..- attr(*, "names")= chr [1:2] "function" "gradient"
##   ..$ convergence: int 0
##   ..$ message    : NULL
##   ..$ hessian    : num [1:2, 1:2] 43.4 -38.3 -38.3 425.3
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "shape" "scale"
##   .. .. ..$ : chr [1:2] "shape" "scale"
##  - attr(*, "class")= chr "flexsurvreg"

Now, fit all RLS data sets by strains

for( i in 1:length(report[,1])){
#for( i in 3:4){
  my.files = files[grep(strains[i], files)]
  tb = read.table( paste("../qinlab_rls/",my.files[1],sep=''), sep="\t")
  if( length(my.files)> 1){
    for( fi in 2:length(my.files)) {
      tmp.tb = read.table( paste("../qinlab_rls/",my.files[fi],sep=''), sep="\t")
      tb = rbind( tb, tmp.tb)
    }
  }
  report$samplesize[i] = length(tb[,1])

  GompFlex = flexsurvreg(formula = Surv(tb[,1]) ~ 1, dist = 'gompertz')
  WeibFlex = flexsurvreg(formula = Surv(tb[,1]) ~ 1, dist = 'weibull')

  report$avgLS[i] =  mean(tb[,1])
  report$stdLS[i] =  sd(tb[,1])
  report$CV[i] = report$stdLS[i] / report$avgLS[i]

  report$GompGFlex[i] = GompFlex$res[1,1]
  report$GompRFlex[i] = GompFlex$res[2,1]
  report$GompLogLikFlex[i] = round(GompFlex$loglik, 1)
  report$GompAICFlex[i] = round(GompFlex$AIC)

  report$WeibShapeFlex[i] = WeibFlex$res[1,1]
  report$WeibRateFlex[i] = WeibFlex$res[2,1]
  report$WeibLogLikFlex[i] = round(WeibFlex$loglik, 1)  
  report$WeibAICFlex[i] = round(WeibFlex$AIC)

  #set initial values
  Rhat = report$GompRFlex[i]; # 'i' was missing. a bug costed HQ a whole afternoon.
  Ghat = report$GompGFlex[i];
  nhat = 6;  
  t0= (nhat-1)/Ghat;
  fitBinom = optim ( c(Rhat, t0, nhat),  llh.binomialMortality.single.run, 
                     lifespan=tb[,1], 
                     #method='SANN') #SANN needs control  
                     method="L-BFGS-B", 
                     lower=c(1E-10, 1, 1), upper=c(1,200,20) );  
  report[i, c("R", "t0", "n")] = fitBinom$par[1:3]
  report$G[i] = (report$n[i] - 1)/report$t0[i]
}
Show the results
report[ grep("BY", report$strains), ]
##    strains samplesize           R        t0        n          G    avgLS
## 2   BY4716         30 0.004214506 68.883068 7.897675 0.10013600 31.76667
## 3   BY4741         91 0.004561640 69.177799 7.921589 0.10005506 30.40659
## 4   BY4742         91 0.004131546  3.446732 2.533993 0.44505728 20.67033
## 5   BY4743         90 0.003172625 63.416617 8.013376 0.11059209 33.40000
## 6 JSBY4741         30 0.004232442 51.481720 6.021519 0.09753984 35.40000
##       stdLS        CV  GompGFlex   GompRFlex GompLogLikFlex GompAICFlex
## 2 12.082142 0.3803403 0.07246963 0.005281817         -118.2         240
## 3 12.213816 0.4016831 0.07215985 0.005922721         -358.1         720
## 4  9.446992 0.4570315 0.06536841 0.016813684         -339.8         684
## 5 10.790966 0.3230828 0.07864221 0.004029916         -349.9         704
## 6 11.352047 0.3206793 0.09711885 0.001921016         -114.5         233
##   WeibShapeFlex WeibRateFlex WeibLogLikFlex WeibAICFlex
## 2      2.909731     35.65364         -116.4         237
## 3      2.672627     34.08892         -357.0         718
## 4      2.286304     23.28727         -330.0         664
## 5      3.342059     37.18081         -341.3         687
## 6      3.661056     39.35057         -114.5         233
Strehler-Mildvan correlation in natural isolates
# my.strains=c("101S", "BY4743","M1-2","M13","M14","M2-8","M22","M32","M34","M5","M8","RM112N","S288c","SGU57","SK1","W303","YPS128","YPS163")

my.strains=c("101S", "M1-2","M13","M14","M2-8","M22","M32","M34","M5","M8","RM112N","S288c","SGU57", "YPS128","YPS163")
report2 = report[ match(my.strains, report$strains), ]
summary(lm(log10(report2$GompRFlex) ~ report2$GompGFlex))
## 
## Call:
## lm(formula = log10(report2$GompRFlex) ~ report2$GompGFlex)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50796 -0.22010 -0.03155  0.20504  0.36320 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -1.6976     0.3278  -5.179 0.000177 ***
## report2$GompGFlex  -8.7148     2.7614  -3.156 0.007585 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2789 on 13 degrees of freedom
## Multiple R-squared:  0.4338, Adjusted R-squared:  0.3902 
## F-statistic:  9.96 on 1 and 13 DF,  p-value: 0.007585
summary(lm(log10(report2$R) ~ report2$GompGFlex))
## 
## Call:
## lm(formula = log10(report2$R) ~ report2$GompGFlex)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.16896 -0.09699  0.01555  0.08156  0.21625 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -1.9873     0.1451 -13.699 4.21e-09 ***
## report2$GompGFlex  -4.7701     1.2223  -3.903  0.00182 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1234 on 13 degrees of freedom
## Multiple R-squared:  0.5395, Adjusted R-squared:  0.5041 
## F-statistic: 15.23 on 1 and 13 DF,  p-value: 0.001817
summary(lm(log10(report2$R) ~ report2$G))
## 
## Call:
## lm(formula = log10(report2$R) ~ report2$G)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.191584 -0.112059  0.002454  0.085419  0.274808 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.0949     0.1574 -13.313 5.96e-09 ***
## report2$G    -3.1369     1.0797  -2.905   0.0123 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1417 on 13 degrees of freedom
## Multiple R-squared:  0.3937, Adjusted R-squared:  0.3471 
## F-statistic: 8.441 on 1 and 13 DF,  p-value: 0.01228
summary(lm(report2$GompGFlex ~ report2$G))
## 
## Call:
## lm(formula = report2$GompGFlex ~ report2$G)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0158806 -0.0099107 -0.0007003  0.0080448  0.0311292 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.02065    0.01524   1.355    0.199    
## report2$G    0.67115    0.10459   6.417 2.28e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01372 on 13 degrees of freedom
## Multiple R-squared:  0.7601, Adjusted R-squared:  0.7416 
## F-statistic: 41.18 on 1 and 13 DF,  p-value: 2.281e-05
plot( report2$t0 ~ report2$n, col='red', xlim=c(5.5, 9), ylim=c(25, 80), pch=19, xlab='n', ylab='t0')
text(report2$n +0.25, report2$t0+0.5, report2$strains, cex=0.75)

report2
##    strains samplesize           R       t0        n          G    avgLS
## 1     101S         85 0.001963826 36.96827 7.286027 0.17003846 31.34118
## 7     M1-2         54 0.002597778 40.44275 7.581458 0.16273519 27.83333
## 8      M13         70 0.002863446 40.56801 7.650792 0.16394178 26.54286
## 9      M14         60 0.004156533 54.65063 6.031192 0.09206101 36.55000
## 10    M2-8        105 0.003438767 42.75291 8.026908 0.16436094 24.80952
## 11     M22         60 0.004669075 46.38267 6.035812 0.10857096 31.83333
## 12     M32         60 0.001856854 35.31281 7.690560 0.18946549 27.96667
## 13     M34         58 0.002238125 31.64343 7.085632 0.19231896 27.01724
## 14      M5        166 0.003348780 74.58364 7.898541 0.09249403 36.62651
## 15      M8         60 0.001863662 31.46248 6.021022 0.15958760 34.93333
## 16  RM112N         59 0.002702475 55.93484 6.023393 0.08980795 44.06780
## 17   S288c         41 0.004995123 57.42117 7.948488 0.12100918 26.26829
## 18   SGU57         58 0.006334776 55.69961 7.714361 0.12054593 23.86207
## 24  YPS128         69 0.002045751 41.96634 6.960180 0.14202285 35.00000
## 25  YPS163        130 0.001820701 36.99506 6.820331 0.15732725 34.43077
##        stdLS        CV  GompGFlex    GompRFlex GompLogLikFlex GompAICFlex
## 1   7.512772 0.2397093 0.13407292 0.0012654562         -295.0         594
## 7   9.201620 0.3305971 0.12310344 0.0024494207         -193.6         391
## 8   9.142488 0.3444425 0.12284897 0.0029261133         -250.1         504
## 9  12.804164 0.3503191 0.09148573 0.0019510502         -233.7         471
## 10  8.133614 0.3278424 0.11608390 0.0043612883         -373.3         751
## 11 10.271182 0.3226549 0.10778993 0.0021454528         -222.9         450
## 12  6.888868 0.2463242 0.14043277 0.0017384578         -205.6         415
## 13  8.206740 0.3037593 0.15692797 0.0012896468         -198.0         400
## 14 12.938747 0.3532618 0.06684902 0.0041472501         -670.1        1344
## 15  6.905823 0.1976858 0.15888831 0.0003653141         -201.9         408
## 16 13.006450 0.2951464 0.08938574 0.0010363304         -232.0         468
## 17 10.254327 0.3903690 0.08686882 0.0064512047         -154.2         312
## 18 10.538898 0.4416590 0.08956763 0.0076748634         -216.5         437
## 24  9.719598 0.2777028 0.11866014 0.0011041546         -252.2         508
## 25  8.591449 0.2495282 0.13387276 0.0007889455         -459.6         923
##    WeibShapeFlex WeibRateFlex WeibLogLikFlex WeibAICFlex
## 1       4.778139     34.15561         -291.0         586
## 7       3.538533     30.92356         -195.6         395
## 8       3.114056     29.26548         -259.1         522
## 9       3.361407     40.80155         -236.8         478
## 10      3.371375     27.61417         -368.5         741
## 11      3.596392     35.37118         -223.9         452
## 12      4.419104     30.62828         -201.4         407
## 13      3.985745     29.77062         -203.7         411
## 14      3.098392     40.98709         -658.3        1321
## 15      5.872213     37.70008         -200.1         404
## 16      4.038465     48.63895         -233.7         471
## 17      2.792964     29.42942         -153.3         311
## 18      2.455703     26.85910         -218.2         440
## 24      4.294122     38.52990         -253.0         510
## 25      4.825655     37.62671         -460.2         924
Gompert versus Weibull? AIC: smaller is better (for information loss)
report$BestModel = ifelse(report$GompAICFlex < report$WeibAICFlex, "Gomp", "Weib")
report$BestModel = ifelse(abs(report$GompAICFlex - report$WeibAICFlex)<2, "<2", report$BestModel)
CV ~ Gomp and Weibull? How does noises influence likelihood of Gompertz and Weibull fitting?
summary(lm(report$CV ~ report$BestModel ))
## 
## Call:
## lm(formula = report$CV ~ report$BestModel)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.133680 -0.055476 -0.003523  0.048974  0.138701 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           0.32019    0.04341   7.376 2.21e-07 ***
## report$BestModelGomp  0.03043    0.05012   0.607    0.550    
## report$BestModelWeib  0.01117    0.04816   0.232    0.819    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07518 on 22 degrees of freedom
## Multiple R-squared:  0.02302,    Adjusted R-squared:  -0.0658 
## F-statistic: 0.2592 on 2 and 22 DF,  p-value: 0.774
summary(lm(report$CV ~ report$WeibLogLikFlex ))
## 
## Call:
## lm(formula = report$CV ~ report$WeibLogLikFlex)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.138661 -0.041551 -0.005703  0.044865  0.152245 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            3.343e-01  3.688e-02   9.063 4.73e-09 ***
## report$WeibLogLikFlex -1.042e-05  1.304e-04  -0.080    0.937    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07438 on 23 degrees of freedom
## Multiple R-squared:  0.0002772,  Adjusted R-squared:  -0.04319 
## F-statistic: 0.006377 on 1 and 23 DF,  p-value: 0.937
summary(lm(report$CV ~ report$GompLogLikFlex  ))
## 
## Call:
## lm(formula = report$CV ~ report$GompLogLikFlex)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.138744 -0.041556 -0.005757  0.044670  0.152321 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            3.346e-01  3.650e-02   9.166 3.85e-09 ***
## report$GompLogLikFlex -9.067e-06  1.281e-04  -0.071    0.944    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07438 on 23 degrees of freedom
## Multiple R-squared:  0.0002176,  Adjusted R-squared:  -0.04325 
## F-statistic: 0.005006 on 1 and 23 DF,  p-value: 0.9442
summary(lm(report$CV ~ (report$GompLogLikFlex - report$WeibLogLikFlex)))
## 
## Call:
## lm(formula = report$CV ~ (report$GompLogLikFlex - report$WeibLogLikFlex))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.138744 -0.041556 -0.005757  0.044670  0.152321 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            3.346e-01  3.650e-02   9.166 3.85e-09 ***
## report$GompLogLikFlex -9.067e-06  1.281e-04  -0.071    0.944    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07438 on 23 degrees of freedom
## Multiple R-squared:  0.0002176,  Adjusted R-squared:  -0.04325 
## F-statistic: 0.005006 on 1 and 23 DF,  p-value: 0.9442
plot( report$GompLogLikFlex ~ report$CV, col="red", pch=3, xlim=c(0, 0.8), ylim=c(-220, -80))
points( report$CV, report$WeibLogLikFlex, col="blue", pch=4)
m1 = lm( report$GompLogLikFlex ~ report$CV)
m2 = lm( report$WeibLogLikFlex ~ report$CV)
abline( m1, col="red", lty=2)
abline( m2, col='blue', lty=1)
text(0.6, -210, "nearly the same!?")

The QIN-RLS data suggested that noisy system signal perfer Gompertz model, based on GG01 theory. Notice that CV measures distrubition of system signals and are different from white noises (residues)
report$DeltaGompWeiLLH = report$GompLogLikFlex - report$WeibLogLikFlex
plot( report$DeltaGompWeiLLH ~ report$CV )
m3 = lm( report$DeltaGompWeiLLH ~ report$CV)
abline(m3, col='red')

TODO: Calculate the white noises (fitting errors) using the fitting residues for Gompertz and Weibull models.
# report$residues ?? 
Output
write.csv(report, file = 'sandbox/_report_qinlab_strains_rls.csv', row.names = FALSE)

write.csv(report2, file = 'sandbox/_report_qinlab_natural_isolates_rls.csv', row.names = FALSE)

No comments:

Post a Comment