Combine analytical solution and numerical solutions.

S_net  = S1 * S2

S1 has a analytic form.  S2 has to be found by simulation.  How do we combine these two in numerical simulation? We can simulate S2 and draw random number from S1.  The network lifespan is the minimal of the two:
   network_lifespan = min(lifespan1, lifespan2).
S_net can then be calculated from the distribution of network_lifespans.

Network heterogeneity and frailty

Binomial distribution should lead to a specific distribution of frailty. Vaupel79's frailty is a normalized ration of mortality rates. The denominator is basically a scaling factor, so distribution of this frailty should be decided by the numerator.

*** Useful R materials for teaching and learning

Useful R materials.

EdX computational genomics using R
http://genomicsclass.github.io/book/pages/classes.html
https://courses.edx.org/courses/HarvardX/PH525.1x/1T2015/info

Data Camp

https://www.datacamp.com/

 

R swirl 

https://www.datacamp.com/swirl-r-tutorial

TryR @ code school (no loops in chapter 1-7!)

• Stephen J. Eglen's PLOS article. A quick guide to teaching R programming to computational biology students. 

• 
  

http://cran.r-project.org/doc/contrib/Seefeld_StatsRBio.pdf

http://www.rseek.org. 

http://addictedtor.free.fr/graphiques/. 


Newly launched interactive cloud-based R for teaching: http://www.datamind.org/


Hong Qin's R and computational biology GitHub site




Why R? 

An Introduction to Statistical Learning

with Applications in R

http://www-bcf.usc.edu/~gareth/ISL/index.html 



Elements of statistical learning, 

http://www.r-bloggers.com/free-stanford-online-course-on-statistical-learning-with-r-starting-on-19-jan-2015/

 http://www.r-bloggers.com/in-depth-introduction-to-machine-learning-in-15-hours-of-expert-videos/
 

Why constant mortality rate is non-aging?

Subpopulation of an exponential distribution is the same distribution. So, as long as individuals can live to the next day, they will start the next day as if they are new-borns.  Everyday is a Groundhog Day in these non-aging systems.

Individuals from different ages have the same risk of dying. Or, some external factor are killing all individuals regardless of their ages.

What does Strehler-Mildvan correlation means?

What does not kill you makes you stronger.

how many viruses, bacteria, plants and animals on earth?

How much viruses on planet earth? 10^31 based on Carl Zimmer's blog.

http://www.virology.ws/2009/10/19/ten-cool-facts-about-viruses/

http://phenomena.nationalgeographic.com/2013/02/20/an-infinity-of-viruses/

Quote from Carl Zimmer's blog:
RNA viruses made up between 38 and 63% of the viruses in the sea water. In other words, about half of the viruses in the ocean are RNA viruses.


viruses invented DNA?

Total number of bacteria on planet earth
http://www.sdearthtimes.com/et0998/et0998s8.html
5x10^30 (Hmm, this is only half of the 10^31 viruses of Carl Zimmer's blog).

How many single-cell organisms on ocean floor? 2.9x10^29
http://www.nature.com/news/there-are-fewer-microbes-out-there-than-you-think-1.11275
(Who this would include prokaryote, archea, and eukayote).

I did not find quotes on how many plants and animals on earth.

June 17, 2015 Biomass of plants, animals, fungi, prokaryotes, viruses are summarized in the following plos biology paper.

http://journals.plos.org/plosbiology/article?id=10.1371%2Fjournal.pbio.1002168

Cox proportional harzard model

Assumptions of the Cox model, from http://www.stat.ubc.ca/~rollin/teach/643w04/lec/node69.html
    Though the Cox model is non-parametric to the extent that no assumptions are made about form of the baseline hazard, there are still a number of important issues which need be assessed before the model results can be safely applied.
    First and foremost is the issue of non-informative censoring. To satisfy this assumption, the design of the underlying study must ensure that the mechanisms giving rise to censoring of individual subjects are not related to the probability of an event occurring. For example, in clinical studies, care must be taken that continuation of follow-up not depend on a participants medical condition. Violation of this assumption can invalidates just about any sort of survival analysis, from Kaplan-Meier estimation to the Cox model.

     The second key assumption in the Cox model is that of proportional hazards. In a regression type setting this means that the survival curves for two strata (determined by the particular choices of values for the -variables) must have hazard functions that are proportional over time (i.e. constant relative hazard). This can be evaluated graphically using "log-log" plots in the two-sample comparison case. In that situation, and also for the Cox model, there are tests that can be applied to test proportionality.  

Other references:
http://courses.washington.edu/b515/l17.pdf

flexsurvreg() in package flexsurv

Notes on 20140818: I did not see the specification of gradient function in the call of optim()?!


# package(flexsurv)

> flexsurvreg
function (formula, data, dist, inits, fixedpars = NULL, cl = 0.95, 
    ...) 
{
    call <- match.call() #a useful base function
    indx <- match(c("formula", "data"), names(call), nomatch = 0)
    if (indx[1] == 0) 
        stop("A \"formula\" argument is required")
    temp <- call[c(1, indx)]
    temp[[1]] <- as.name("model.frame")
    m <- eval(temp, parent.frame())
    Y <- model.extract(m, "response")
    if (!inherits(Y, "Surv")) 
        stop("Response must be a survival object")
    Terms <- attr(m, "terms")
    X <- model.matrix(Terms, m)
    dat <- list(Y = Y, X = X[, -1, drop = FALSE], Xraw = m[, 
        -1, drop = FALSE])
    X <- dat$X
    if (missing(dist)) 
        stop("Distribution \"dist\" not specified")
    if (is.character(dist)) {
        match.arg(dist, names(flexsurv.dists))
        dlist <- flexsurv.dists[[dist]]
        dlist$name <- dist
    }
    else if (is.list(dist)) {
        dlist <- check.dlist(dist)
    }
    else stop("\"dist\" should be a string for a built-in distribution, or a list for a custom distribution")
    parnames <- dlist$pars
    ncovs <- ncol(X)
    nbpars <- length(parnames)
    npars <- nbpars + ncovs
    if (!missing(inits) && (!is.numeric(inits) || (length(inits) != 
        npars))) 
        stop("inits must be a numeric vector of length ", npars)
    if (missing(inits) || any(is.na(inits))) 
        default.inits <- c(dlist$inits(Y[, "time"]), rep(0, ncovs))
    if (missing(inits)) 
        inits <- default.inits
    else if (any(is.na(inits))) 
        inits[is.na(inits)] <- default.inits[is.na(inits)]
    for (i in 1:nbpars) inits[i] <- dlist$transforms[[i]](inits[i])
    outofrange <- which(is.nan(inits) | is.infinite(inits))
    if (any(outofrange)) {
        plural <- if (length(outofrange) > 1) 
            "s"
        else ""
        stop("Initial value", plural, " for parameter", plural, 
            " ", paste(outofrange, collapse = ","), " out of range")
    }
    cnames <- if (ncovs == 0) 
        NULL
    else colnames(X)
    names(inits) <- c(parnames, cnames)
    if (!is.null(fixedpars) && !is.logical(fixedpars) && (!is.numeric(fixedpars) || 
        any(!(fixedpars %in% 1:npars)))) {
        dots <- if (npars > 2) 
            "...,"
        else ""
        stop("fixedpars must be TRUE/FALSE or a vector of indices in 1,", 
            dots, npars)
    }
    if ((is.logical(fixedpars) && fixedpars == TRUE) || (is.numeric(fixedpars) && 
        all(fixedpars == 1:npars))) {
        minusloglik <- minusloglik.flexsurv(inits, t = Y[, "time"], 
            dead = Y[, "status"], X = X, dlist = dlist, inits = inits)
        for (i in 1:nbpars) inits[i] <- dlist$inv.transforms[[i]](inits[i])
        res <- matrix(inits, ncol = 1)
        dimnames(res) <- list(names(inits), "est")
        ret <- list(call = call, dlist = dlist, res = res, npars = 0, 
            loglik = -minusloglik, AIC = 2 * minusloglik, data = dat, 
            datameans = colMeans(dat$X), N = nrow(dat$Y), events = sum(dat$Y[, 
                "status"]), trisk = sum(dat$Y[, "time"]))
    }
    else {
        optpars <- inits[setdiff(1:npars, fixedpars)]  #nice trick to optimize a subset of pars
        opt <- optim(optpars, minusloglik.flexsurv, t = Y[, "time"], 
            dead = Y[, "status"], X = X, dlist = dlist, inits = inits, 
            fixedpars = fixedpars, hessian = TRUE, ...) #fixedpars is a useful optim feature
        est <- opt$par
        if (all(!is.na(opt$hessian)) && all(!is.nan(opt$hessian)) && 
            all(is.finite(opt$hessian)) && all(eigen(opt$hessian)$values > 
            0)) {
            cov <- solve(opt$hessian)
            se <- sqrt(diag(cov)) #qin, standard errors, see wiki entry 
            if (!is.numeric(cl) || length(cl) > 1 || !(cl > 0) || 
                !(cl < 1)) 
                stop("cl must be a number in (0,1)")
            lcl <- est - qnorm(1 - (1 - cl)/2) * se
            ucl <- est + qnorm(1 - (1 - cl)/2) * se
        }
        else {
            warning("Could not calculate asymptotic standard errors - Hessian is not positive definite. Optimisation has probably not converged to the maximum likelihood")
            lcl <- ucl <- NA
        }
        res <- cbind(est = inits, lcl = NA, ucl = NA)
        res[setdiff(1:npars, fixedpars), ] <- cbind(est, lcl, 
            ucl)
        colnames(res) <- c("est", paste(c("L", "U"), round(cl * 
            100), "%", sep = ""))
        for (i in 1:nbpars) res[i, ] <- dlist$inv.transforms[[i]](res[i, 
            ])
        ret <- list(call = match.call(), dlist = dlist, res = res, 
            npars = length(est), loglik = -opt$value, AIC = 2 * 
                opt$value + 2 * length(est), cl = cl, opt = opt, 
            data = dat, datameans = colMeans(dat$X), N = nrow(dat$Y), 
            events = sum(dat$Y[, "status"]), trisk = sum(dat$Y[, 
                "time"]))
    }
    class(ret) <- "flexsurvreg"
    ret
}

flexsurv, R package,

R package 'flexsurv' provides optim() based parametric fitting. Two parameter Gomertz and Weibull model are included.  The 'flexsurv' use Hessian to build confidence intervals around the maximum likelihood estimate.

I did not see an obvious way to do nested model test using 'flexsurv'.

Block, Li, Savits, 2003, initial and final behaviour of failure rate functions for mixtures and systems

The initial approximation of mixture failure rate is  straight forward. For each subpopulation pdf = failure rate *Viability.  Because viability approaches 1 initially, mixture of pdf leads to mixture of failure rates.

Reference:

Initial and Final Behaviour of Failure Rate Functions for Mixtures and Systems
Author(s): Henry W. Block, Yulin Li and Thomas H. Savits Source: Journal of Applied Probability, Vol. 40, No. 3 (Sep., 2003), pp. 721-740 Published by: Applied Probability Trust

Stable URL: http://www.jstor.org/stable/3215946 .

Gamma function, a smooth curve that connect the factorial points

From wikipedia:
The gamma function can be seen as a solution to the following interpolation problem:

"Find a smooth curve that connects the points (x, y) given by y = (x − 1)! at the positive integer values for x."

 

Although there are multiple solutions, gamma function is a popular choice.

Resources on scientific writing

Writing steps:
  Big ideas
  Outlines
  Parageraph arrangement
  Subject, verb, object
  Plurals, singular, tense, pronouns

http://www.americanscientist.org/issues/pub/the-science-of-scientific-writing

http://www.biomedcentral.com/authors/authoracademy

http://www.biomedcentral.com/authors/report

Presentation tips (draft)

• Pay attention to seating of audience and do not block their views.  
• Some default line plots generated by R can be very dim on screens, especially for the dashed lines. 

notes on dbSNP, in progress

Multiple genome locations can be mapped to a single SNP.  One reasons is the ambiguity of alignment.
http://www.biostars.org/p/2323/


https://cgsmd.isi.edu/dbsnpq/downloads.php

Notes, Lai & Xie, 2006, stochastic ageing and dependence for reliability

page 48, mixture failure rate r(t):
 r(t) = weighted sum of subpopulation pdf / weighted sum of subpopulation survival function
 This can be found by the definition of r(t).  Failure rate is mortality rate in biology.

page 57, initial failure mixture rate r(0+):
 r(0+) = weighted sum of subpopulation r_i(0+),
shown by Block, Li, Savit, 2003, initial and final behavior of failure rate function for mixtures and systems, Journal of Applied Probability, 40, 721-740
 I am surprised that this basic formula has a 2003 citation.  This formula was in GG91 book, and they cited the Barlow, Proschan, (Hunter) book.