Saturday, October 19, 2013

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.

Wednesday, October 16, 2013

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.


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 $x$-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'.


Tuesday, October 15, 2013

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

Saturday, October 12, 2013

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. 

Friday, October 11, 2013

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. 





Thursday, October 10, 2013

Q-and-A on network reliability model of aging

1. How does studying someone's individual characteristics of aging relate to studying the mortality rate?
Mortality rate describes the probability of dying. Declining vigor of individual is modeled by their 'increasing chance of dying".

2. Is yeast being studied instead of using an actual human cell?
Yeast cells are used as models for human cells.

3. Is NRMCA the program we will be using to study cells?
NRMCA is a model. We will have to write codes based on NRMCA for simulations.

4. What does it mean to be an emergent property of gene networks? ex. "Our central hypothesis is that cellular aging is an emergent property of gene networks (pg. 2)".
The question is basically what is an emergent property. It is a concept that has not been well-defined. It is generally agreed that emergent property is a property found at higher system level, but not at its component levels.

5. "The age of the network is determined by the minimal age of the modules" - does this mean that we are studying the age of network but that is determined by the components making up this network? Thus we are not studying specific cells but instead networks of cells.
Networks are made of essential modules in NRMCA, each module contain one essential genes and n number of non-essential genes. A network is equivalent to a cell.

6. Does essential genes with fewer interactions are the limiting steps mean that fewer interactions implies a higher mortality rate?
Yes.

7. Which equations are most relevant to my research?
This question needs to be more specific.
 
Vocabulary
1. Stochastically: things occur by chances

2. Stochastic heterogeneity: noises due to random chances
3. Ziemman Function: this is a mathematical function
4. ubiquitous: everywhere
5. SIR2: a gene's name
6. rDNA circles: ribosomal DNA cicles
7. TOR pathway: TOR is a gene's name. TOR pathway refers to its signal cascading pathway.
8. canalization: Literally means all things come together. A metaphorical word for the biological phenomena that different genotypes lead to the same phenotype.

Sunday, October 6, 2013

Petrovski, PlosGenetics, mutation tolerance in human genes, data resources

PlosGenetics, Genic Intolerance to Functional Variation and the Interpretation of Personal Genomes, Slavé Petrovski et al, 2013

The paper measured mutations using exome sequencing, conduct network analysis.

http://channelopathist.net/2013/10/06/mutation-intolerance-why-some-genes-withstand-mutations-and-others-dont/?utm_source=buffer&utm_campaign=Buffer&utm_content=buffer94829&utm_medium=twitter

http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1003709