Likelihood Ratios for Medical Decision Making

diagnostic.decision <- function(Sensitivity, Specificity,
                                Prevalence = 0.01,
                                AE.Risk = 0.005,
                                Predicate.Sensitivity = 0.9, 
                                Predicate.Specificity = 0.9, 
                                Pretest.Prob = c(0.1, 0.25, 0.5, 0.75, 0.9),
                                Cost.Benefit.Ratio = 1
                                )
{
    ## Purpose: Provide metrics to evaluate a medical test. 
    ## Arguments:
    ##   Sensitivity, Specificity: between 0 and 1.
    ##   Prevalence: Disease Prevalence in the Target Population
    ##   Pretest.Prob: a vector of pre-test probability
    ##   AE.Risk: Risk of an adverse event due to treatment for False Positive patients.
    ##   Predicate.Sensitivity, Predicate.Specificity: Performance of the comparator.
    ##   Cost.Benefit.Ratio: Cost/Benefit = Cost(False Positive) / Cost(False Negative); default to 1. 
    ## Return: Likelihood Ratios, Post-test Probability, Plots, etc.

    ## Author: Feiming Chen
    ## ________________________________________________

    dev.close()
    p2o <- function(p) p/(1 - p)        # convert probability to odds
    o2p <- function(o) o/(1 + o)        # convert odds to probability
    r <- function(x) round(x, 2)           # rounding function

    cat("\n----------- Diagnostic Performance for Subject Device ------------\n")
    cat("Sensitivity =", r(Sensitivity), "\n")   
    cat("Specificity =", r(Specificity), "\n\n") 

    cat("Risk Difference Evaluation:\n")
    J <- Sensitivity + Specificity - 1  # Risk Difference = Youden's Index = J
    cat("J =", r(J), "\n\n") 

    LR.pos <- Sensitivity / (1 - Specificity) # Positive Likelihood Ratio (Risk Ratio)
    LR.neg <- (1 - Sensitivity) / Specificity # Negative Likelihood Ratio
    cat("Risk Ratio Evaluation:\n")
    cat("LR+ =", r(LR.pos), "\n")   
    cat("LR- =", r(LR.neg), "\n\n") 

    cat("Predictive Values:\n")
    PPV <- o2p(Sensitivity / (1 - Specificity) * p2o(Prevalence))
    NPV <- o2p(Specificity / (1 - Sensitivity) / p2o(Prevalence))
    cat("PPV = ", r(PPV), "\n")
    cat("NPV = ", r(NPV), "\n")

    cat("\nDecision Analysis:\n")
    cat("Cost/Benefit Ratio =", Cost.Benefit.Ratio, "\n")
    Risk.Threshold <- o2p(Cost.Benefit.Ratio)
    cat("Implied Risk Threshold =", r(Risk.Threshold), "\n")

    cat("\nRegret Graph Analysis:\n")
    abcissas.1 <- o2p(Cost.Benefit.Ratio / LR.pos)
    abcissas.2 <- o2p(Cost.Benefit.Ratio / LR.neg)
    cat("Disease Prevalence =", Prevalence, "\n")
    cat("Range of Disease Prevalence (Pre-Test Probability) to Choose a Test Strategy:\n")
    cat("0 --", r(abcissas.1), ": Test for None (Presume Negative; Treatment None)\n")
    cat(r(abcissas.1),"--", r(abcissas.2), ": Subject Device Test\n")
    cat(r(abcissas.2), "-- 1 : Test for All (Presume Positive; Treatment All)\n")    

    cat("\n----------- Diagnostic Performance for Predicate Device ------------\n")
    cat("Sensitivity =", r(Predicate.Sensitivity), "\n")   
    cat("Specificity =", r(Predicate.Specificity), "\n\n") 

    cat("Risk Difference Evaluation:\n")
    J <- Predicate.Sensitivity + Predicate.Specificity - 1  # Risk Difference = Youden's Index = J
    cat("J =", r(J), "\n\n") 

    LR2.pos <- Predicate.Sensitivity / (1 - Predicate.Specificity) # Positive Likelihood Ratio (Risk Ratio)
    LR2.neg <- (1 - Predicate.Sensitivity) / Predicate.Specificity # Negative Likelihood Ratio
    cat("Risk Ratio Evaluation:\n")
    cat("LR+ =", r(LR2.pos), "\n")   
    cat("LR- =", r(LR2.neg), "\n\n") 

    cat("\n----------- Compare the Subject and Predicate Device ------------\n")
    cat("Likelihood Ratio Graph:\n")
    plot(c(1 - Predicate.Specificity, 1 - Specificity), c(Predicate.Sensitivity, Sensitivity), 
         xlim = c(0, 1), ylim = c(0, 1), asp = 1, xlab = "False Positive Rate", ylab = "True Positive Rate", pch = 19, col = c("red", "blue"), bty = "n", 
         main = "Likelihood Ratio Graph")
    legend("bottomright", legend = c("Predicate", "Subject"), pch = 19, bg="lightyellow", col = c("red", "blue"), title = "Device")
    abline(0, LR.pos)                 
    min.spec <- 1 - 1/LR.pos           
    cat("Minimum Specificity for a new test to be superior to the subject device in terms of PPV and NPV =", r(min.spec), "\n")

    min.sens <- 1 - LR.neg
    abline(min.sens, LR.neg)
    cat("Minimum Sensitivity for a new test to be superior to the subject device in terms of PPV and NPV =", r(min.sens), "\n")

    cat("\n----------- Expected Diagnostic Yield ------------\n")
    cat("Assume 1,000,000 patients from the target populaiton\n")
    cat("Assume the Disease Prevalence =", Prevalence, "\n")
    cat("Assume the Risk of an Adverse Event due to Treatment for False Positive Patients =", AE.Risk, "\n")
    cat("Predicate Device Performance: Sensitivity =", Predicate.Sensitivity, ", Specificity =", Predicate.Specificity, "\n")
    cat("Subject Device Performance: Sensitivity =", Sensitivity, ", Specificity =", Specificity, "\n\n")
    n <- 1000000
    n1 <- round(n * Prevalence)         # number with disease
    n2 <- n - n1                        # number without disease

    EDY <- matrix(c(n1, n2,
                    n1 * Sensitivity, n2 * (1 - Specificity),
                    n1 * Predicate.Sensitivity, n2 * (1 - Predicate.Specificity),
                    n1 * Sensitivity - n1 * Predicate.Sensitivity, 
                    n2 * (1 - Specificity) - n2 * (1 - Predicate.Specificity)
                    ), 2)
    FP.per.TP <- apply(EDY, 2, function(x) round(x[2]/(x[1] + 0.01)))
    FP.per.TP.Pct <- FP.per.TP / FP.per.TP[1] * 100
    EDY1 <- round(rbind(EDY, FP.per.TP, FP.per.TP.Pct))
    dimnames(EDY1) <-  list(c("Disease", "No.Disease", "Number.of.FP.Per.TP", "Percent.of.Worst.Case"), c("Expected.Number", "Subject+", "Predicate+", "Difference"))
    print(EDY1)
    cat("\nSafety Evaluation:\n")
    SAF <- matrix(round(c(n2,
                          n2 * (1 - Specificity) * AE.Risk,
                          n2 * (1 - Predicate.Specificity) * AE.Risk,
                          n2 * (1 - Specificity) * AE.Risk - n2 * (1 - Predicate.Specificity) * AE.Risk, 
                          (n2 * (1 - Specificity) * AE.Risk - n2 * (1 - Predicate.Specificity) * AE.Risk) * 100 /
                          (n1 * Sensitivity - n1 * Predicate.Sensitivity + 0.001))), 1,
                  dimnames = list(c("No.Disease"), c("Expected.Number", "Subject+", "Predicate+", "Difference", "Number.of.AE.Per.100.Extra.TP")))
    print(SAF)


    cat("\n----------- Medical Decision Making ------------\n")
    cat("Pre-test Probability =", Pretest.Prob, "\n")
    Pretest.Odds <- p2o(Pretest.Prob)
    cat("Pre-test Odds =", r(Pretest.Odds), "\n\n")
    Posttest.Odds.pos <- Pretest.Odds * LR.pos # with positive diagnostic result
    Posttest.Odds.neg <- Pretest.Odds * LR.neg # with negative diagnostic result
    Posttest.Prob.pos <- o2p(Posttest.Odds.pos)
    Posttest.Prob.neg <- o2p(Posttest.Odds.neg)

    cat("IF the diagnostic test is Positive:\n")
    cat("Post-test Odds =", r(Posttest.Odds.pos), "\n")
    cat("Post-test Probability =", r(Posttest.Prob.pos), "\n\n")

    cat("IF the diagnostic test is Negative:\n")
    cat("Post-test Odds =", r(Posttest.Odds.neg), "\n")
    cat("Post-test Probability =", r(Posttest.Prob.neg), "\n\n")

    N <- length(Pretest.Prob)
    x <- rep(0, N)
    
    dev.new()
    plot(c(-1.1, 1.1), c(0, 1), type="n", xlab="", ylab="Disease Probability", axes=F,
         main="Pre-Test Probability is updated to Post-Test Probability,\nafter a Diagnostic Test.",
         sub=paste("Diagnostic Test: Sensitivity =", r(Sensitivity), ", Specificity =", r(Specificity), 
                   ", LR+ =", r(LR.pos), ", LR- =", r(LR.neg)))
    arrows(x, Pretest.Prob, x-1, Posttest.Prob.neg, length=0.15)
    arrows(x, Pretest.Prob, x+1, Posttest.Prob.pos, length=0.15)
    points(x, Pretest.Prob, pch=19)
    points(x-1, Posttest.Prob.neg, pch=19)
    points(x+1, Posttest.Prob.pos, pch=19)
    axis(side=1, at=c(-1, 0, 1), labels=c("Post-Test (Test = NEG)", "Pre-Test", "Post-Test (Test = POS)"))
    axis(side=2, at=seq(0, 1, 0.1))
    o <- 0.2
    text(x, Pretest.Prob, labels=Pretest.Prob, pos=3, offset=o)
    text(x-1, Posttest.Prob.neg, labels=round(Posttest.Prob.neg, 2), pos=2, offset=o)
    text(x+1, Posttest.Prob.pos, labels=round(Posttest.Prob.pos, 2), pos=4, offset=o)

    list(EDY = EDY1, SAF = SAF)
}
if (F) {                                # Unit Test
    a <- diagnostic.decision(Sensitivity = 0.64, 
                             Specificity = 0.84, 
                             Prevalence = 0.0022,
                             AE.Risk = 0.0068,
                             Predicate.Sensitivity = 0.224,
                             Predicate.Specificity = 0.881)

}

----------- Diagnostic Performance for Subject Device ------------
Sensitivity = 0.64 
Specificity = 0.84 

Risk Difference Evaluation:
J = 0.48 

Risk Ratio Evaluation:
LR+ = 4 
LR- = 0.43 

Predictive Values:
PPV =  0.01 
NPV =  1 

Decision Analysis:
Cost/Benefit Ratio = 1 
Implied Risk Threshold = 0.5 

Regret Graph Analysis:
Disease Prevalence = 0.0022 
Range of Disease Prevalence (Pre-Test Probability) to Choose a Test Strategy:
0 -- 0.2 : Test for None (Presume Negative; Treatment None)
0.2 -- 0.7 : Subject Device Test
0.7 -- 1 : Test for All (Presume Positive; Treatment All)

----------- Diagnostic Performance for Predicate Device ------------
Sensitivity = 0.22 
Specificity = 0.88 

Risk Difference Evaluation:
J = 0.1 

Risk Ratio Evaluation:
LR+ = 1.88 
LR- = 0.88 


----------- Compare the Subject and Predicate Device ------------
Likelihood Ratio Graph:
Minimum Specificity for a new test to be superior to the subject device in terms of PPV and NPV = 0.75 
Minimum Sensitivity for a new test to be superior to the subject device in terms of PPV and NPV = 0.57 

----------- Expected Diagnostic Yield ------------
Assume 1,000,000 patients from the target populaiton
Assume the Disease Prevalence = 0.0022 
Assume the Risk of an Adverse Event due to Treatment for False Positive Patients = 0.0068 
Predicate Device Performance: Sensitivity = 0.224 , Specificity = 0.881 
Subject Device Performance: Sensitivity = 0.64 , Specificity = 0.84 

                      Expected.Number Subject+ Predicate+ Difference
Disease                          2200     1408        493        915
No.Disease                     997800   159648     118738      40910
Number.of.FP.Per.TP               454      113        241         45
Percent.of.Worst.Case             100       25         53         10

Safety Evaluation:
           Expected.Number Subject+ Predicate+ Difference Number.of.AE.Per.100.Extra.TP
No.Disease          997800     1086        807        278                            30

----------- Medical Decision Making ------------
Pre-test Probability = 0.1 0.25 0.5 0.75 0.9 
Pre-test Odds = 0.11 0.33 1 3 9 

IF the diagnostic test is Positive:
Post-test Odds = 0.44 1.33 4 12 36 
Post-test Probability = 0.31 0.57 0.8 0.92 0.97 

IF the diagnostic test is Negative:
Post-test Odds = 0.05 0.14 0.43 1.29 3.86 
Post-test Probability = 0.05 0.12 0.3 0.56 0.79 

Bland Altman Limit of Agreement Calculation with the presence of repeated measurements where the true value vary

Refer to this blog post for the function var.comp.one.way.model

bland.altman.repeated.vary <- function(fmla, dat)
{
    ## Purpose: Calculate LoA (Limit of Agreement) when data has repeated measurements (true value varies)
    ## Arguments:
    ##    fmla: something like "Difference ~ Subject".
    ##    dat: data frame with "Difference" and "Subject".
    ## Return: LoA
    ## Author: Feiming Chen
    ## Reference:
    ##     Bland, J. Martin, and Douglas G. Altman. "Agreement between methods of measurement with multiple observations per
    ##     individual." Journal of biopharmaceutical statistics 17.4 (2007): 571-582.
    ## ________________________________________________

    ans <- var.comp.one.way.model(fmla, dat)
    s <- ans$VCA$SD[3]
    m <- ans$Mean
    cat("\nMean =", m, ", SD =", s, "\n")
    loa <- c(m - 1.96*s, m + 1.96*s) 
    cat("\n95% Limits of Agreement: (", loa[1], ",", loa[2], ")\n")
    invisible(list(m=m, s=s, loa=loa))
}
if (F) {                                # Unit Test
    dat <- data.frame(Subject = factor(c(1, 1, 1, 2, 2, 2, 3, 3, 3, 3)), A = rnorm(10, 10), B = rnorm(10, 10))
    dat$Difference <- dat$A - dat$B
    fmla <- Difference ~ Subject
    bland.altman.repeated.vary(fmla, dat)
    ## Mean = 0.0997215 , SD = 1.851 
    ## 95% Limits of Agreement: ( -3.528239 , 3.727681 )
}

Calculate Nonparametric 95% Reference Interval

calc.reference.interval <- function(val, wgt)
{
    ## Purpose: Calculate Nonparametric 95% Reference Interval. 
    ## Notes: 
    ##    1) Each Reference Interval needs at least 120 reference observations (per CLSI C28-A3c).
    ##    2) A general criterion for determining sample size is that the width of the 90% confidence interval (CI)
    ##       for a reference limit should be acceptably small relative to the width of the 95% reference interval
    ##       itself (Harris and Boyd); these authors recommend that the width of the 90% CI be less than 0.2 times
    ##       the width of the reference interval. 
    ## Arguments:
    ##   val: numerical measurements
    ##   wgt: weight (frequency) of these measurements
    ## Return: a Nonparametric 95% Reference Interval. 
    ## Author: Feiming Chen
    ## ________________________________________________

    x <- rep(val, wgt)
    cat("N =", length(x), "\n")
    cat("Reference Interval:\n")
    print(q <- quantile(x, probs = c(0.025, 0.975)))

    ## Confidence Intervals for the Nonparametric Method (Using Bootstrap with 2000 Reps)
    require(boot)
    b <- boot(x, 
              function(x, i) {
                  quantile(x[i], probs = c(0.025, 0.975))
              }, R = 2000)

    cat("\n90% Bootstrap CI for Lower 95% Reference Limit:")
    L <- boot.ci(b, conf=0.90, type="basic", index=1)
    print(L2 <- L$basic[1,c(4,5)])

    cat("\n90% Bootstrap CI for Upper 95% Reference Limit:")
    U <- boot.ci(b, conf=0.90, type="basic", index=2)
    print(U2 <- U$basic[1,c(4,5)])

    cat("\nRatio of Limit CI Width to Reference Interval Width (Recommend: < 0.2):\n")
    q2 <- diff(q)                       # Width of Reference Interval
    L3 <- diff(L2)                      # Width of Lower Limit 90% CI
    U3 <- diff(U2)                      # Width of Upper Limit 90% CI
    cat("Lower Limit CI Width Ratio:", round(L3 / q2, 3), "\n")
    cat("Upper Limit CI Width Ratio:", round(U3 / q2, 3), "\n")

    list(Reference.Limit = q, 
         Lower.Limit.CI = L2,
         Upper.Limit.CI = U2)
}
if (F) {                                # Unit Test
    ## Data from CLSI C28-A3c (Reference Interval) Table 4 (Frequency Distributions of Calcium in 240 Medical Students by Sex)
    ## Only Women's Data. CLSI Reported "Nonparametric" Method: N = 120, Calcium.Women = (8.9, 10.2),
    ##                    Lower Limit 90% CI = (8.8, 9.1), Upper Limit 90% CI = (10.1, 10.3)
    val <- c(8.8,8.9,9.0,9.1,9.2,9.3,9.4,9.5,9.6,9.7,9.8,9.9,10.0,10.1,10.2,10.3,10.4,10.5,10.6)
    wgt <- c(1,2,1,3,11,11,8,16,16,26,8,7,3,2,3,2,0,0,0)
    ans <- calc.reference.interval(val, wgt)
    ##   2.5%   97.5% 
    ## 8.9975 10.2000 
}

Variance Component Analysis for One-Way Random Model

var.comp.one.way.model <- function(fmla, dat)
{
    ## Purpose: Variance Component Analysis for One-Way Random Model. 
    ##          Application: Repeatability and Reproducibility Analysis (R&R for Precision Study).
    ##          Use Method of Moment estimation.  
    ## Arguments:
    ##    fmla: a formula with one response and one factor (One-way Random Model).
    ##    dat: a data frame
    ## Return: Variance Component Analysis such as
    ##         DF: Degrees of Freedom
    ##         VC: Variance Component (with proportion of total in percentage)
    ##         SD: Standard Deviation (from the squared root of Variance Component)
    ##         CV%: Coefficient of Variation as Percentage
    ## Author: Feiming Chen
    ## ________________________________________________

    a <- aov(fmla, dat)
    b <- summary(a)
    cat("Analysis of Variance:\n")
    print(b)                   # Analysis of Variance
    MEAN <- mean(a$model[[1]])  # Mean 
    N <- nrow(a$model)          # Number of observations
    cat("\nMean: ", round(MEAN, 3), " (N =", N, ")\n", sep="")
    
    d <- b[[1]]
    VC.error <- d$"Mean Sq"[2]          # Mean Sum of Squares for Error = Variance Component for Error (Repeatability)
    Q2 <- d$"Sum Sq"[1]                 # Sum of Squares for Treatment
    DF.treatment <- d$Df[1]             # DF for Treatment
    N.i <- table(a$model[[2]])          # Number of Replicates in Each Treatment Choice
    VC.treatment <-  (Q2 - DF.treatment * VC.error) / (N - sum(N.i^2) / N)

    VC <- c(VC.treatment, VC.error, VC.treatment + VC.error)
    VC.pct.total <- VC / VC[3] * 100
    SD <- sqrt(VC)
    CV.pct <- SD / MEAN * 100
    
    ans <- data.frame(Name = c(attr(a$terms, "term.labels"), "Error", "Total"), 
                      VC = round(VC, 3), 
                      VC.Pct.Total = round(VC.pct.total, 3), 
                      SD = round(SD, 3), 
                      CV.Pct = round(CV.pct, 3))

    cat("\nVariance Component Analysis:\n")
    ans
}
if (F) {                                # Unit Test
    dat <- data.frame(Yield=c(3.9, 4.05, 4.25, 3.60, 4.2, 4.05, 3.85, 4.15, 4.6, 4.15, 4.4, 3.35, 3.8), 
                      Variety=factor(c(1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4)))
    fmla <- Yield ~ Variety
    var.comp.one.way.model(fmla, dat)

    ## Validation
    library(VCA)
    anovaVCA(fmla, dat)
}

Compare the speed of product vs exp/sum/log calculations

## Compare the speed of product vs exp/sum/log calculations

library(microbenchmark)
x <- runif(1000) + 0.5
all.equal(prod(x), exp(sum(log(x))))
## [1] TRUE                                

microbenchmark::microbenchmark(
                    prod(x), 
                    exp(sum(log(x)))
)

## Unit: microseconds
##             expr    min     lq     mean  median      uq     max neval
##          prod(x)  2.075  2.107  2.52904  2.1945  2.5285  13.964   100
## exp(sum(log(x))) 51.525 52.075 62.89344 52.6135 62.5120 400.676   100

We can see it is much faster to use the straightforward calculation of "prod".

Use R Environment Object to Share and Modify (Big) Data Within Functions

We can use the R environment object just like a list.  One advantage is that it will be passed by reference to functions.  So it can be modified inside functions without the need to make additional copies.  Refer to Environment for more detailed explanation.  Following is an example usage.

## Use Reference Semantics to pass data object to functions by reference.
dat <- new.env(parent = emptyenv())

## Put in some data
dat$iris <- iris
dat$cars <- cars

## Checking the original object. 
ls(dat)                                
## [1] "cars" "iris"

## An example function
f <- function(d) d$N <- nrow(d$iris) + nrow(d$cars)

## Passing the data by reference: It can be modified inside the function.
f(dat)

## Checking the modified object. 
ls(dat)                                
## [1] "cars" "iris" "N"   
dat$N
## [1] 200
identical(dat$N, nrow(dat$cars) + nrow(dat$iris))
## [1] TRUE

## In practice, we can use f1(dat), f2(dat), f3(dat), ...
## to continuously update the information/states maintained in "dat".

Examining R Profiling Data

rprof <- function(srcfile)
{ 
    ## Purpose: R Profiling with R "proftools" package. 
    ## Arguments:
    ##   srcfile: path to source file. 
    ## Return: plots and summary of profiling. 
    ## Author: Feiming Chen, Date: 23 May 2018, 10:31
    ## ________________________________________________

    library(proftools)
    pd0 <- profileExpr(source(srcfile))
    pd <- filterProfileData(pd0, select = "withVisible", skip = 4)

    cat("\nSummarize by function -------------------------------------\n")
    print(head(funSummary(pd, srclines = FALSE), 10))

    cat("\nSummarize by call -----------------------------------------\n")
    print(head(callSummary(pd), 10))

    cat("\nHot Execution Paths ---------------------------------------\n")
    print(hotPaths(pd, total.pct = 10.0))

    plotProfileCallGraph(pd, maxnodes = 15)

    flameGraph(pd)

    calleeTreeMap(pd)

    pd
}
if (F) {                                # Unit Test
    srcfile <- system.file("samples", "bootlmEx.R", package = "proftools")
    rprof(srcfile)
}

The output of the Unit Test is as follows:

Summarize by function -------------------------------------
                    total.pct gc.pct self.pct gcself.pct
statistic               80.57   7.21     0.66       0.00
lapply                  80.57   7.21     1.75       0.00
FUN                     80.57   7.21     1.97       0.44
boot                    80.57   7.21     0.00       0.00
glm                     49.13   5.24     0.87       0.22
eval                    32.75   1.97     1.97       0.22
model.frame.default     27.51   1.75     1.31       0.00
sapply                  26.42   3.06     0.87       0.00
predict                 24.45   1.53     0.44       0.00
predict.glm             24.02   1.53     0.44       0.00

Summarize by call -----------------------------------------
                                       total.pct gc.pct self.pct gcself.pct
lapply -> FUN                              80.57   7.21     1.09       0.22
FUN -> statistic                           80.57   7.21     0.66       0.00
statistic -> glm (bootlmEx.R:33)           49.13   5.24     0.87       0.22
boot (bootlmEx.R:39) -> lapply             38.21   3.93     0.00       0.00
boot (bootlmEx.R:44) -> lapply             37.55   3.06     0.00       0.00
eval -> eval                               30.13   1.97     0.87       0.22
glm (bootlmEx.R:33) -> eval                27.73   1.75     0.87       0.00
statistic -> predict (bootlmEx.R:35)       24.45   1.53     0.44       0.00
predict (bootlmEx.R:35) -> predict.glm     24.02   1.53     0.44       0.00
sapply -> lapply                           23.14   3.06     1.75       0.00

Hot Execution Paths ---------------------------------------
 path                            total.pct self.pct
 boot (bootlmEx.R:39)            38.21     0.00    
 . lapply                        38.21     0.00    
 . . FUN                         38.21     0.00    
 . . . statistic                 38.21     0.00    
 . . . . glm (bootlmEx.R:33)     24.67     0.44    
 . . . . . eval                  13.76     0.44    
 . . . . . . eval                13.32     0.22    
 . . . . predict (bootlmEx.R:35) 12.23     0.00    
 . . . . . predict.glm           12.01     0.44    
 . . . . . . predict.lm          11.35     0.22    
 boot (bootlmEx.R:44)            37.55     0.00    
 . lapply                        37.55     0.00    
 . . FUN                         37.55     0.00    
 . . . statistic                 37.55     0.00    
 . . . . glm (bootlmEx.R:33)     24.45     0.44    
 . . . . . eval                  13.97     0.44    
 . . . . . . eval                13.54     0.00    
 . . . . predict (bootlmEx.R:35) 12.23     0.22    

 . . . . . predict.glm           12.01     0.00    
 . . . . . . predict.lm          11.35     0.44    
 lm (bootlmEx.R:52)              10.92     0.00