0.fit_qinlabrls.Rmd
h qin
June 9, July 6-11, 2016
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.5source("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.001092234GompFlex$res.t## est L95% U95% se
## shape 0.1464128 0.09875196 0.1940735 0.02431718
## rate -6.4361044 -7.77186373 -5.1003451 0.68152238GompFlex$opt$hessian## shape rate
## shape 23564.3361 810.0597
## rate 810.0597 30.0000str(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]
}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# 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.007585summary(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.001817summary(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.01228summary(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-05plot( 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 924report$BestModel = ifelse(report$GompAICFlex < report$WeibAICFlex, "Gomp", "Weib")
report$BestModel = ifelse(abs(report$GompAICFlex - report$WeibAICFlex)<2, "<2", report$BestModel)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.774summary(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.937summary(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.9442summary(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.9442plot( 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 ?? 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