# Data

## Coefficient of variation (CV)

Motivating the coefficient of variation (CV) for beginners:

• Boss: Measure it 5 times.
• You: 8, 8, 9, 6, and 8
• B: SD=1. Make it three times more precise!
• Y: 0.20 0.20 0.23 0.15 0.20 meters. SD=0.3!
• B: All you did was change to meters! Report the CV instead!
• Y: Damn it.
R> sd(c(8, 8, 9, 6, 8))
[1] 1.095445
R> sd(c(8, 8, 9, 6, 8)*2.54/100)
[1] 0.02782431


## Transform sample values to their percentiles

set.seed(1234)
x <- rnorm(10)
x
# [1] -1.2070657  0.2774292  1.0844412 -2.3456977  0.4291247  0.5060559
# [7] -0.5747400 -0.5466319 -0.5644520 -0.8900378
ecdf(x)(x)
# [1] 0.2 0.7 1.0 0.1 0.8 0.9 0.4 0.6 0.5 0.3

rank(x)
# [1]  2  7 10  1  8  9  4  6  5  3


# Box(Box and whisker) plot in R

See

An example for a graphical explanation.

> x=c(0,4,15, 1, 6, 3, 20, 5, 8, 1, 3)
> summary(x)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0       2       4       6       7      20
> sort(x)
[1]  0  1  1  3  3  4  5  6  8 15 20
> boxplot(x, col = 'grey')

# https://en.wikipedia.org/wiki/Quartile#Example_1
> summary(c(6, 7, 15, 36, 39, 40, 41, 42, 43, 47, 49))
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
6.00   25.50   40.00   33.18   42.50   49.00

• The lower and upper edges of box is determined by the first and 3rd quartiles (2 and 7 in the above example).
• 2 = median(c(0, 1, 1, 3, 3, 4)) = (1+3)/2
• 7 = median(c(4, 5, 6, 8, 15, 20)) = (6+8)/2
• IQR = 7 - 2 = 5
• The thick dark horizon line is the median (4 in the example).
• Outliers are defined by (the empty circles in the plot)
• Observations larger than 3rd quartile + 1.5 * IQR (7+1.5*5=14.5) and
• smaller than 1st quartile - 1.5 * IQR (2-1.5*5=-5.5).
• Note that the cutoffs are not shown in the Box plot.
• Whisker (defined using the cutoffs used to define outliers)
• Upper whisker is defined by the largest "data" below 3rd quartile + 1.5 * IQR (8 in this example), and
• Lower whisker is defined by the smallest "data" greater than 1st quartile - 1.5 * IQR (0 in this example).
• See another example below where we can see the whiskers fall on observations.

Note the wikipedia lists several possible definitions of a whisker. R uses the 2nd method (Tukey boxplot) to define whiskers.

## Create boxplots from a list object

Normally we use a vector to create a single boxplot or a formula on a data to create boxplots.

But we can also use split() to create a list and then make boxplots.

## geom_boxplot

Without jitter

ggplot(dfbox, aes(x=sample, y=expr)) +
geom_boxplot() +
theme(axis.text.x=element_text(color = "black", angle=30, vjust=.8,
hjust=0.8, size=6),
plot.title = element_text(hjust = 0.5)) +
labs(title="", y = "", x = "")


With jitter

ggplot(dfbox, aes(x=sample, y=expr)) +
geom_boxplot(outlier.shape=NA) + #avoid plotting outliers twice
geom_jitter(position=position_jitter(width=.2, height=0)) +
theme(axis.text.x=element_text(color = "black", angle=30, vjust=.8,
hjust=0.8, size=6),
plot.title = element_text(hjust = 0.5)) +
labs(title="", y = "", x = "")


What do hjust and vjust do when making a plot using ggplot? The value of hjust and vjust are only defined between 0 and 1: 0 means left-justified, 1 means right-justified.

# stem and leaf plot

stem(). See R Tutorial.

Note that stem plot is useful when there are outliers.

> stem(x)

The decimal point is 10 digit(s) to the right of the |

0 | 00000000000000000000000000000000000000000000000000000000000000000000+419
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
11 |
12 | 9

> max(x)
[1] 129243100275
> max(x)/1e10
[1] 12.92431

> stem(y)

The decimal point is at the |

0 | 014478
1 | 0
2 | 1
3 | 9
4 | 8

> y
[1] 3.8667356428 0.0001762708 0.7993462430 0.4181079732 0.9541728562
[6] 4.7791262101 0.6899313108 2.1381289177 0.0541736818 0.3868776083

> set.seed(1234)
> z <- rnorm(10)*10
> z
[1] -12.070657   2.774292  10.844412 -23.456977   4.291247   5.060559
[7]  -5.747400  -5.466319  -5.644520  -8.900378
> stem(z)

The decimal point is 1 digit(s) to the right of the |

-2 | 3
-1 | 2
-0 | 9665
0 | 345
1 | 1


# Don't invert that matrix

## Different matrix decompositions/factorizations

set.seed(1234)
x <- matrix(rnorm(10*2), nr= 10)
cmat <- cov(x); cmat
# [,1]       [,2]
# [1,]  0.9915928 -0.1862983
# [2,] -0.1862983  1.1392095

# cholesky decom
d1 <- chol(cmat)
t(d1) %*% d1  # equal to cmat
d1  # upper triangle
# [,1]       [,2]
# [1,] 0.9957875 -0.1870864
# [2,] 0.0000000  1.0508131

# svd
d2 <- svd(cmat)
## $dose ## NLvMH NvL MvH ## None -0.5 0.5 0.0 ## Low -0.5 -0.5 0.0 ## Med 0.5 0.0 0.5 ## High 0.5 0.0 -0.5 library(dplyr) dose.means <- summarize(group_by(data, dose), y.mean=mean(y)) dose.means ## Source: local data frame [4 x 2] ## ## dose y.mean ## 1 None 112.6267 ## 2 Low 121.3500 ## 3 Med 126.7839 ## 4 High 113.5517 # The coefficient estimate for the first contrast (3.18) equals the average of # the last two groups (126.78 + 113.55 /2 = 120.17) minus the average of # the first two groups (112.63 + 121.35 /2 = 116.99).  ## Multicollinearity > op <- options(contrasts = c("contr.helmert", "contr.poly")) > npk.aov <- aov(yield ~ block + N*P*K, npk) > alias(npk.aov) Model : yield ~ block + N * P * K Complete : (Intercept) block1 block2 block3 block4 block5 N1 P1 K1 N1:P1 N1:K1 P1:K1 N1:P1:K1 0 1 1/3 1/6 -3/10 -1/5 0 0 0 0 0 0 > options(op)  ## Exposure Independent variable = predictor = explanatory = exposure variable ## Confounders, confounding ## Causal inference ## Confidence interval vs prediction interval Confidence intervals tell you about how well you have determined the mean E(Y). Prediction intervals tell you where you can expect to see the next data point sampled. That is, CI is computed using Var(E(Y|X)) and PI is computed using Var(E(Y|X) + e). ## Heteroskedasticity ## Linear regression with Map Reduce ## Relationship between multiple variables ## Model fitting evaluation ## Generalized least squares # Quantile regression # Non- and semi-parametric regression ## Mean squared error ## Splines ## k-Nearest neighbor regression • k-NN regression in practice: boundary problem, discontinuities problem. • Weighted k-NN regression: want weight to be small when distance is large. Common choices - weight = kernel(xi, x) ## Kernel regression • Instead of weighting NN, weight ALL points. Nadaraya-Watson kernel weighted average: ${\displaystyle {\hat {y}}_{q}=\sum c_{qi}y_{i}/\sum c_{qi}={\frac {\sum {\text{Kernel}}_{\lambda }({\text{distance}}(x_{i},x_{q}))*y_{i}}{\sum {\text{Kernel}}_{\lambda }({\text{distance}}(x_{i},x_{q}))}}}$. • Choice of bandwidth ${\displaystyle \lambda }$ for bias, variance trade-off. Small ${\displaystyle \lambda }$ is over-fitting. Large ${\displaystyle \lambda }$ can get an over-smoothed fit. Cross-validation. • Kernel regression leads to locally constant fit. • Issues with high dimensions, data scarcity and computational complexity. # Principal component analysis See PCA. # Partial Least Squares (PLS) ${\displaystyle X=TP^{\mathrm {T} }+E}$ ${\displaystyle Y=UQ^{\mathrm {T} }+F}$ where X is an ${\displaystyle n\times m}$ matrix of predictors, Y is an ${\displaystyle n\times p}$ matrix of responses; T and U are ${\displaystyle n\times l}$ matrices that are, respectively, projections of X (the X score, component or factor matrix) and projections of Y (the Y scores); P and Q are, respectively, ${\displaystyle m\times l}$ and ${\displaystyle p\times l}$ orthogonal loading matrices; and matrices E and F are the error terms, assumed to be independent and identically distributed random normal variables. The decompositions of X and Y are made so as to maximise the covariance between T and U (projection matrices). PLS, PCR (principal components regression) and ridge regression tend to behave similarly. Ridge regression may be preferred because it shrinks smoothly, rather than in discrete steps. # High dimension Partial least squares prediction in high-dimensional regression Cook and Forzani, 2019 ## Feature selection # Independent component analysis ICA is another dimensionality reduction method. ## ICA vs PCA ## ICS vs FA # Correspondence analysis # Non-negative matrix factorization # t-SNE t-Distributed Stochastic Neighbor Embedding (t-SNE) is a technique for dimensionality reduction that is particularly well suited for the visualization of high-dimensional datasets. library(ggplot2) library(Rtsne) URL <- "https://assets.datacamp.com/production/repositories/1680/datasets/68b37d6c5f7f6768d5e11796687993b6f3da1f72/mnist-sample-200.RData" load(url(URL)) dim(mnist_sample) # 200 obs, 28x28 pixels # [1] 200 785 mnist_sample[1:2, 1:5] # label pixel0 pixel1 pixel2 pixel3 # 1 5 0 0 0 0 # 2 0 0 0 0 0 set.seed(1234) # tsne_output <- Rtsne(mnist[, -1]) tsne <- Rtsne(mnist_sample[, -1], PCA = FALSE, dims = 3, perplexity = 50, max_iter = 1300) head(tsne$Y)
tsne$itercosts tsne_out <- Rtsne(mnist_sample[, -1], perplexity = 5) tsne_plot <- data.frame(tsne_x= tsne_out$Y[,1],
tsne_y = tsne_out$Y[,2], digit = as.factor(mnist_sample$label))
ggplot(tsne_plot, aes(x= tsne_x, y = tsne_y, color = digit)) +
ggtitle("MNIST embedding of the first 200 digits") +
geom_text(aes(label = digit)) +
theme(legend.position="none")

library(data.table)
centroids <- as.data.table(tsne_out$Y) setnames(centroids, c("X", "Y")) centroids[, label := as.factor(mnist_sample$label)]

centroids[, mean_X := mean(X), by = label]
centroids[, mean_Y := mean(Y), by = label]
centroids <- unique(centroids, by = "label")

ggplot(centroids, aes(x= mean_X, y = mean_Y, color = label)) +
ggtitle("Centroids coordinates") +
geom_text(aes(label = label)) +
theme(legend.position = "none")


# Calibration

• Search by image: graphical explanation of calibration problem
• https://www.itl.nist.gov/div898/handbook/pmd/section1/pmd133.htm Calibration and calibration curve.
• Y=voltage (observed), X=temperature (true/ideal). The calibration curve for a thermocouple is often constructed by comparing thermocouple (observed)output to relatively (true)precise thermometer data.
• when a new temperature is measured with the thermocouple, the voltage is converted to temperature terms by plugging the observed voltage into the regression equation and solving for temperature.
• It is important to note that the thermocouple measurements, made on the secondary measurement scale, are treated as the response variable and the more precise thermometer results, on the primary scale, are treated as the predictor variable because this best satisfies the underlying assumptions (Y=observed, X=true) of the analysis.
• Calibration interval
• In almost all calibration applications the ultimate quantity of interest is the true value of the primary-scale measurement method associated with a measurement made on the secondary scale.
• It seems the x-axis and y-axis have similar ranges in many application.
• An Exercise in the Real World of Design and Analysis, Denby, Landwehr, and Mallows 2001. Inverse regression
• How to determine calibration accuracy/uncertainty of a linear regression?
• Linear Regression and Calibration Curves
• Regression and calibration Shaun Burke
• calibrate package
• investr: An R Package for Inverse Estimation. Paper
• The index of prediction accuracy: an intuitive measure useful for evaluating risk prediction models by Kattan and Gerds 2018. The following code demonstrates Figure 2.
# Odds ratio =1 and calibrated model
set.seed(666)
x = rnorm(1000)
z1 = 1 + 0*x
pr1 = 1/(1+exp(-z1))
y1 = rbinom(1000,1,pr1)
mean(y1) # .724, marginal prevalence of the outcome
dat1 <- data.frame(x=x, y=y1)
newdat1 <- data.frame(x=rnorm(1000), y=rbinom(1000, 1, pr1))

# Odds ratio =1 and severely miscalibrated model
set.seed(666)
x = rnorm(1000)
z2 =  -2 + 0*x
pr2 = 1/(1+exp(-z2))
y2 = rbinom(1000,1,pr2)
mean(y2) # .12
dat2 <- data.frame(x=x, y=y2)
newdat2 <- data.frame(x=rnorm(1000), y=rbinom(1000, 1, pr2))

library(riskRegression)
lrfit1 <- glm(y ~ x, data = dat1, family = 'binomial')
IPA(lrfit1, newdata = newdat1)
#     Variable     Brier           IPA     IPA.gain
# 1 Null model 0.1984710  0.000000e+00 -0.003160010
# 2 Full model 0.1990982 -3.160010e-03  0.000000000
# 3          x 0.1984800 -4.534668e-05 -0.003114664
1 - 0.1990982/0.1984710
# [1] -0.003160159

lrfit2 <- glm(y ~ x, family = 'binomial')
IPA(lrfit2, newdata = newdat1)
#     Variable     Brier       IPA     IPA.gain
# 1 Null model 0.1984710  0.000000 -1.859333763
# 2 Full model 0.5674948 -1.859334  0.000000000
# 3          x 0.5669200 -1.856437 -0.002896299
1 - 0.5674948/0.1984710
# [1] -1.859334

From the simulated data, we see IPA = -3.16e-3 for a calibrated model and IPA = -1.86 for a severely miscalibrated model.

See ROC.

# Generalized Linear Model

Lectures from a course in Simon Fraser University Statistics.

## Quasi Likelihood

Quasi-likelihood is like log-likelihood. The quasi-score function (first derivative of quasi-likelihood function) is the estimating equation.

## Deviance, stats::deviance() and glmnet::deviance.glmnet() from R

• It is a generalization of the idea of using the sum of squares of residuals (RSS) in ordinary least squares to cases where model-fitting is achieved by maximum likelihood. See What is Deviance? (specifically in CART/rpart) to manually compute deviance and compare it with the returned value of the deviance() function from a linear regression. Summary: deviance() = RSS in linear models.
• https://www.rdocumentation.org/packages/stats/versions/3.4.3/topics/deviance
• Likelihood ratio tests and the deviance http://data.princeton.edu/wws509/notes/a2.pdf#page=6
• Deviance(y,muhat) = 2*(loglik_saturated - loglik_proposed)
• Interpreting Residual and Null Deviance in GLM R
• Null Deviance = 2(LL(Saturated Model) - LL(Null Model)) on df = df_Sat - df_Null. The null deviance shows how well the response variable is predicted by a model that includes only the intercept (grand mean).
• Residual Deviance = 2(LL(Saturated Model) - LL(Proposed Model)) = ${\displaystyle 2(LL(y|y)-LL({\hat {\mu }}|y))}$, df = df_Sat - df_Proposed=n-p. ==> deviance() has returned.
• Null deviance > Residual deviance. Null deviance df = n-1. Residual deviance df = n-p.
## an example with offsets from Venables & Ripley (2002, p.189)
utils::data(anorexia, package = "MASS")

anorex.1 <- glm(Postwt ~ Prewt + Treat + offset(Prewt),
family = gaussian, data = anorexia)
summary(anorex.1)

# Call:
#   glm(formula = Postwt ~ Prewt + Treat + offset(Prewt), family = gaussian,
#       data = anorexia)
#
# Deviance Residuals:
#   Min        1Q    Median        3Q       Max
# -14.1083   -4.2773   -0.5484    5.4838   15.2922
#
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  49.7711    13.3910   3.717 0.000410 ***
#   Prewt        -0.5655     0.1612  -3.509 0.000803 ***
#   TreatCont    -4.0971     1.8935  -2.164 0.033999 *
#   TreatFT       4.5631     2.1333   2.139 0.036035 *
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# (Dispersion parameter for gaussian family taken to be 48.69504)
#
# Null deviance: 4525.4  on 71  degrees of freedom
# Residual deviance: 3311.3  on 68  degrees of freedom
# AIC: 489.97
#
# Number of Fisher Scoring iterations: 2

deviance(anorex.1)
# [1] 3311.263

• In glmnet package. The deviance is defined to be 2*(loglike_sat - loglike), where loglike_sat is the log-likelihood for the saturated model (a model with a free parameter per observation). Null deviance is defined to be 2*(loglike_sat -loglike(Null)); The NULL model refers to the intercept model, except for the Cox, where it is the 0 model. Hence dev.ratio=1-deviance/nulldev, and this deviance method returns (1-dev.ratio)*nulldev.
x=matrix(rnorm(100*2),100,2)
y=rnorm(100)
fit1=glmnet(x,y)
deviance(fit1)  # one for each lambda
#  [1] 98.83277 98.53893 98.29499 98.09246 97.92432 97.78472 97.66883
#  [8] 97.57261 97.49273 97.41327 97.29855 97.20332 97.12425 97.05861
# ...
# [57] 96.73772 96.73770
fit2 <- glmnet(x, y, lambda=.1) # fix lambda
deviance(fit2)
# [1] 98.10212
deviance(glm(y ~ x))
# [1] 96.73762
sum(residuals(glm(y ~ x))^2)
# [1] 96.73762


# Simulate data

## Density plot

# plot a Weibull distribution with shape and scale
func <- function(x) dweibull(x, shape = 1, scale = 3.38)
curve(func, .1, 10)

func <- function(x) dweibull(x, shape = 1.1, scale = 3.38)
curve(func, .1, 10)


The shape parameter plays a role on the shape of the density function and the failure rate.

• Shape <=1: density is convex, not a hat shape.
• Shape =1: failure rate (hazard function) is constant. Exponential distribution.
• Shape >1: failure rate increases with time

## Signal to noise ratio

${\displaystyle {\frac {\sigma _{signal}^{2}}{\sigma _{noise}^{2}}}={\frac {Var(f(X))}{Var(e)}}}$ if Y = f(X) + e

Some examples of signal to noise ratio

## Effect size, Cohen's d and volcano plot

${\displaystyle \theta ={\frac {\mu _{1}-\mu _{2}}{\sigma }},}$

# Multiple comparisons

Take an example, Suppose 550 out of 10,000 genes are significant at .05 level

1. P-value < .05 ==> Expect .05*10,000=500 false positives
2. False discovery rate < .05 ==> Expect .05*550 =27.5 false positives
3. Family wise error rate < .05 ==> The probablity of at least 1 false positive <.05

According to Lifetime Risk of Developing or Dying From Cancer, there is a 39.7% risk of developing a cancer for male during his lifetime (in other words, 1 out of every 2.52 men in US will develop some kind of cancer during his lifetime) and 37.6% for female. So the probability of getting at least one cancer patient in a 3-generation family is 1-.6**3 - .63**3 = 0.95.

## False Discovery Rate

Suppose ${\displaystyle p_{1}\leq p_{2}\leq ...\leq p_{n}}$. Then

${\displaystyle {\text{FDR}}_{i}={\text{min}}(1,n*p_{i}/i)}$.

So if the number of tests (${\displaystyle n}$) is large and/or the original p value (${\displaystyle p_{i}}$) is large, then FDR can hit the value 1.

However, the simple formula above does not guarantee the monotonicity property from the FDR. So the calculation in R is more complicated. See How Does R Calculate the False Discovery Rate.

Below is the histograms of p-values and FDR (BH adjusted) from a real data (Pomeroy in BRB-ArrayTools).

And the next is a scatterplot w/ histograms on the margins from a null data.

## q-value

q-value is defined as the minimum FDR that can be attained when calling that feature significant (i.e., expected proportion of false positives incurred when calling that feature significant).

If gene X has a q-value of 0.013 it means that 1.3% of genes that show p-values at least as small as gene X are false positives.

Another view: q-value = FDR adjusted p-value. A p-value of 5% means that 5% of all tests will result in false positives. A q-value of 5% means that 5% of significant results will result in false positives. here.

## SAM/Significance Analysis of Microarrays

The percentile option is used to define the number of falsely called genes based on 'B' permutations. If we use the 90-th percentile, the number of significant genes will be less than if we use the 50-th percentile/median.

In BRCA dataset, using the 90-th percentile will get 29 genes vs 183 genes if we use median.

## Multivariate permutation test

In BRCA dataset, using 80% confidence gives 116 genes vs 237 genes if we use 50% confidence (assuming maximum proportion of false discoveries is 10%). The method is published on EL Korn, JF Troendle, LM McShane and R Simon, Controlling the number of false discoveries: Application to high dimensional genomic data, Journal of Statistical Planning and Inference, vol 124, 379-398 (2004).

## Empirical Bayes Normal Means Problem with Correlated Noise

The package cashr and the source code of the paper

# offset() function

## Offset in Poisson regression

1. We need to model rates instead of counts
2. More generally, you use offsets because the units of observation are different in some dimension (different populations, different geographic sizes) and the outcome is proportional to that dimension.

An example from here

Y  <- c(15,  7, 36,  4, 16, 12, 41, 15)
N  <- c(4949, 3534, 12210, 344, 6178, 4883, 11256, 7125)
x1 <- c(-0.1, 0, 0.2, 0, 1, 1.1, 1.1, 1)
x2 <- c(2.2, 1.5, 4.5, 7.2, 4.5, 3.2, 9.1, 5.2)

glm(Y ~ offset(log(N)) + (x1 + x2), family=poisson) # two variables
# Coefficients:
# (Intercept)           x1           x2
#     -6.172       -0.380        0.109
#
# Degrees of Freedom: 7 Total (i.e. Null);  5 Residual
# Null Deviance:	    10.56
# Residual Deviance: 4.559 	AIC: 46.69
glm(Y ~ offset(log(N)) + I(x1+x2), family=poisson)  # one variable
# Coefficients:
# (Intercept)   I(x1 + x2)
#   -6.12652      0.04746
#
# Degrees of Freedom: 7 Total (i.e. Null);  6 Residual
# Null Deviance:	    10.56
# Residual Deviance: 8.001 	AIC: 48.13


## Offset in Cox regression

An example from biospear::PCAlasso()

coxph(Surv(time, status) ~ offset(off.All), data = data)
# Call:  coxph(formula = Surv(time, status) ~ offset(off.All), data = data)
#
# Null model
#   log likelihood= -2391.736
#   n= 500

# versus without using offset()
coxph(Surv(time, status) ~ off.All, data = data)
# Call:
# coxph(formula = Surv(time, status) ~ off.All, data = data)
#
#          coef exp(coef) se(coef)    z    p
# off.All 0.485     1.624    0.658 0.74 0.46
#
# Likelihood ratio test=0.54  on 1 df, p=0.5
# n= 500, number of events= 438
coxph(Surv(time, status) ~ off.All, data = data)loglik # [1] -2391.702 -2391.430 # initial coef estimate, final coef  ## Offset in linear regression # Overdispersion Var(Y) = phi * E(Y). If phi > 1, then it is overdispersion relative to Poisson. If phi <1, we have under-dispersion (rare). ## Heterogeneity The Poisson model fit is not good; residual deviance/df >> 1. The lack of fit maybe due to missing data, covariates or overdispersion. Subjects within each covariate combination still differ greatly. Consider Quasi-Poisson or negative binomial. ## Test of overdispersion or underdispersion in Poisson models ## Negative Binomial The mean of the Poisson distribution can itself be thought of as a random variable drawn from the gamma distribution thereby introducing an additional free parameter. ## Binomial # Count data ## Zero counts ## Bias Bias in Small-Sample Inference With Count-Data Models Blackburn 2019 # Survival data analysis # Logistic regression ## Simulate binary data from the logistic model set.seed(666) x1 = rnorm(1000) # some continuous variables x2 = rnorm(1000) z = 1 + 2*x1 + 3*x2 # linear combination with a bias pr = 1/(1+exp(-z)) # pass through an inv-logit function y = rbinom(1000,1,pr) # bernoulli response variable #now feed it to glm: df = data.frame(y=y,x1=x1,x2=x2) glm( y~x1+x2,data=df,family="binomial")  ## Building a Logistic Regression model from scratch ## Odds ratio Calculate the odds ratio from the coefficient estimates; see this post. require(MASS) N <- 100 # generate some data X1 <- rnorm(N, 175, 7) X2 <- rnorm(N, 30, 8) X3 <- abs(rnorm(N, 60, 30)) Y <- 0.5*X1 - 0.3*X2 - 0.4*X3 + 10 + rnorm(N, 0, 12) # dichotomize Y and do logistic regression Yfac <- cut(Y, breaks=c(-Inf, median(Y), Inf), labels=c("lo", "hi")) glmFit <- glm(Yfac ~ X1 + X2 + X3, family=binomial(link="logit")) exp(cbind(coef(glmFit), confint(glmFit)))  ## AUC  predict.glm() ROCR::prediction() ROCR::performance() glmobj ------------> predictTest -----------------> ROCPPred ---------> AUC newdata labels  # Medical applications ## Subgroup analysis Other related keywords: recursive partitioning, randomized clinical trials (RCT) ## Interaction analysis # Statistical Learning ## LDA (Fisher's linear discriminant), QDA ## Bagging Chapter 8 of the book. • Bootstrap mean is approximately a posterior average. • Bootstrap aggregation or bagging average: Average the prediction over a collection of bootstrap samples, thereby reducing its variance. The bagging estimate is defined by ${\displaystyle {\hat {f}}_{bag}(x)={\frac {1}{B}}\sum _{b=1}^{B}{\hat {f}}^{*b}(x).}$ ## Boosting ### AdaBoost AdaBoost.M1 by Freund and Schapire (1997): The error rate on the training sample is ${\displaystyle {\bar {err}}={\frac {1}{N}}\sum _{i=1}^{N}I(y_{i}\neq G(x_{i})),}$ Sequentially apply the weak classification algorithm to repeatedly modified versions of the data, thereby producing a sequence of weak classifiers ${\displaystyle G_{m}(x),m=1,2,\dots ,M.}$ The predictions from all of them are combined through a weighted majority vote to produce the final prediction: ${\displaystyle G(x)=sign[\sum _{m=1}^{M}\alpha _{m}G_{m}(x)].}$ Here ${\displaystyle \alpha _{1},\alpha _{2},\dots ,\alpha _{M}}$ are computed by the boosting algorithm and weight the contribution of each respective ${\displaystyle G_{m}(x)}$. Their effect is to give higher influence to the more accurate classifiers in the sequence. ### Dropout regularization ### Gradient boosting ## Gradient descent Gradient descent is a first-order iterative optimization algorithm for finding the minimum of a function (Wikipedia). The error function from a simple linear regression looks like {\displaystyle {\begin{aligned}Err(m,b)&={\frac {1}{N}}\sum _{i=1}^{n}(y_{i}-(mx_{i}+b))^{2},\\\end{aligned}}} We compute the gradient first for each parameters. {\displaystyle {\begin{aligned}{\frac {\partial Err}{\partial m}}&={\frac {2}{n}}\sum _{i=1}^{n}-x_{i}(y_{i}-(mx_{i}+b)),\\{\frac {\partial Err}{\partial b}}&={\frac {2}{n}}\sum _{i=1}^{n}-(y_{i}-(mx_{i}+b))\end{aligned}}} The gradient descent algorithm uses an iterative method to update the estimates using a tuning parameter called learning rate. new_m &= m_current - (learningRate * m_gradient) new_b &= b_current - (learningRate * b_gradient)  After each iteration, derivative is closer to zero. Coding in R for the simple linear regression. ### Gradient descent vs Newton's method ## Classification and Regression Trees (CART) ### Construction of the tree classifier • Node proportion ${\displaystyle p(1|t)+\dots +p(6|t)=1}$ where ${\displaystyle p(j|t)}$ define the node proportions (class proportion of class j on node t. Here we assume there are 6 classes. • Impurity of node t ${\displaystyle i(t)}$ is a nonnegative function ${\displaystyle \phi }$ of the ${\displaystyle p(1|t),\dots ,p(6|t)}$ such that ${\displaystyle \phi (1/6,1/6,\dots ,1/6)}$ = maximumm ${\displaystyle \phi (1,0,\dots ,0)=0,\phi (0,1,0,\dots ,0)=0,\dots ,\phi (0,0,0,0,0,1)=0}$. That is, the node impurity is largest when all classes are equally mixed together in it, and smallest when the node contains only one class. • Gini index of impurity ${\displaystyle i(t)=-\sum _{j=1}^{6}p(j|t)\log p(j|t).}$ • Goodness of the split s on node t ${\displaystyle \Delta i(s,t)=i(t)-p_{L}i(t_{L})-p_{R}i(t_{R}).}$ where ${\displaystyle p_{R}}$ are the proportion of the cases in t go into the left node ${\displaystyle t_{L}}$ and a proportion ${\displaystyle p_{R}}$ go into right node ${\displaystyle t_{R}}$. A tree was grown in the following way: At the root node ${\displaystyle t_{1}}$, a search was made through all candidate splits to find that split ${\displaystyle s^{*}}$ which gave the largest decrease in impurity; ${\displaystyle \Delta i(s^{*},t_{1})=\max _{s}\Delta i(s,t_{1}).}$ • Class character of a terminal node was determined by the plurality rule. Specifically, if ${\displaystyle p(j_{0}|t)=\max _{j}p(j|t)}$, then t was designated as a class ${\displaystyle j_{0}}$ terminal node. ### R packages ## Partially additive (generalized) linear model trees ## Supervised Classification, Logistic and Multinomial ## Variable selection ### Review Variable selection – A review and recommendations for the practicing statistician by Heinze et al 2018. ### Variable selection and variable importance plot ### Variable selection and cross-validation ### Mallow Cp Mallows's Cp addresses the issue of overfitting. The Cp statistic calculated on a sample of data estimates the mean squared prediction error (MSPE). ${\displaystyle E\sum _{j}({\hat {Y}}_{j}-E(Y_{j}\mid X_{j}))^{2}/\sigma ^{2},}$ The Cp statistic is defined as ${\displaystyle C_{p}={SSE_{p} \over S^{2}}-N+2P.}$ • https://en.wikipedia.org/wiki/Mallows%27s_Cp • Used in Yuan & Lin (2006) group lasso. The degrees of freedom is estimated by the bootstrap or perturbation methods. Their paper mentioned the performance is comparable with that of 5-fold CV but is computationally much faster. ### Variable selection for mode regression http://www.tandfonline.com/doi/full/10.1080/02664763.2017.1342781 Chen & Zhou, Journal of applied statistics ,June 2017 ### lmSubsets lmSubsets: Exact variable-subset selection in linear regression. 2020 ## Neural network ## Support vector machine (SVM) ## Quadratic Discriminant Analysis (qda), KNN ## Regularization Regularization is a process of introducing additional information in order to solve an ill-posed problem or to prevent overfitting Regularization: Ridge, Lasso and Elastic Net from datacamp.com. Bias and variance trade-off in parameter estimates was used to lead to the discussion. ### Ridge regression Since L2 norm is used in the regularization, ridge regression is also called L2 regularization. Hoerl and Kennard (1970a, 1970b) introduced ridge regression, which minimizes RSS subject to a constraint ${\displaystyle \sum |\beta _{j}|^{2}\leq t}$. Note that though ridge regression shrinks the OLS estimator toward 0 and yields a biased estimator ${\displaystyle {\hat {\beta }}=(X^{T}X+\lambda X)^{-1}X^{T}y}$ where ${\displaystyle \lambda =\lambda (t)}$, a function of t, the variance is smaller than that of the OLS estimator. The solution exists if ${\displaystyle \lambda >0}$ even if ${\displaystyle n. Ridge regression (L2 penalty) only shrinks the coefficients. In contrast, Lasso method (L1 penalty) tries to shrink some coefficient estimators to exactly zeros. This can be seen from comparing the coefficient path plot from both methods. Geometrically (contour plot of the cost function), the L1 penalty (the sum of absolute values of coefficients) will incur a probability of some zero coefficients (i.e. some coefficient hitting the corner of a diamond shape in the 2D case). For example, in the 2D case (X-axis=${\displaystyle \beta _{0}}$, Y-axis=${\displaystyle \beta _{1}}$), the shape of the L1 penalty ${\displaystyle |\beta _{0}|+|\beta _{1}|}$ is a diamond shape whereas the shape of the L2 penalty (${\displaystyle \beta _{0}^{2}+\beta _{1}^{2}}$) is a circle. ### Lasso/glmnet, adaptive lasso and FAQs ### Lasso logistic regression ### Lagrange Multipliers ### How to solve lasso/convex optimization • Convex Optimization by Boyd S, Vandenberghe L, Cambridge 2004. It is cited by Zhang & Lu (2007). The interior point algorithm can be used to solve the optimization problem in adaptive lasso. • Review of gradient descent: • Finding maximum: ${\displaystyle w^{(t+1)}=w^{(t)}+\eta {\frac {dg(w)}{dw}}}$, where ${\displaystyle \eta }$ is stepsize. • Finding minimum: ${\displaystyle w^{(t+1)}=w^{(t)}-\eta {\frac {dg(w)}{dw}}}$. • What is the difference between Gradient Descent and Newton's Gradient Descent? Newton's method requires ${\displaystyle g''(w)}$, more smoothness of g(.). • Finding minimum for multiple variables (gradient descent): ${\displaystyle w^{(t+1)}=w^{(t)}-\eta \Delta g(w^{(t)})}$. For the least squares problem, ${\displaystyle g(w)=RSS(w)}$. • Finding minimum for multiple variables in the least squares problem (minimize ${\displaystyle RSS(w)}$): ${\displaystyle {\text{partial}}(j)=-2\sum h_{j}(x_{i})(y_{i}-{\hat {y}}_{i}(w^{(t)}),w_{j}^{(t+1)}=w_{j}^{(t)}-\eta \;{\text{partial}}(j)}$ • Finding minimum for multiple variables in the ridge regression problem (minimize ${\displaystyle RSS(w)+\lambda \|w\|_{2}^{2}=(y-Hw)'(y-Hw)+\lambda w'w}$): ${\displaystyle {\text{partial}}(j)=-2\sum h_{j}(x_{i})(y_{i}-{\hat {y}}_{i}(w^{(t)}),w_{j}^{(t+1)}=(1-2\eta \lambda )w_{j}^{(t)}-\eta \;{\text{partial}}(j)}$. Compared to the closed form approach: ${\displaystyle {\hat {w}}=(H'H+\lambda I)^{-1}H'y}$ where 1. the inverse exists even N<D as long as ${\displaystyle \lambda >0}$ and 2. the complexity of inverse is ${\displaystyle O(D^{3})}$, D is the dimension of the covariates. • Cyclical coordinate descent was used (vignette) in the glmnet package. See also coordinate descent. The reason we call it 'descent' is because we want to 'minimize' an objective function. • ${\displaystyle {\hat {w}}_{j}=\min _{w}g({\hat {w}}_{1},\cdots ,{\hat {w}}_{j-1},w,{\hat {w}}_{j+1},\cdots ,{\hat {w}}_{D})}$ • See paper on JSS 2010. The Cox PHM case also uses the cyclical coordinate descent method; see the paper on JSS 2011. • Coursera's Machine learning course 2: Regression at 1:42. Soft-thresholding the coefficients is the key for the L1 penalty. The range for the thresholding is controlled by ${\displaystyle \lambda }$. Note to view the videos and all materials in coursera we can enroll to audit the course without starting a trial. • No step size is required as in gradient descent. • Implementing LASSO Regression with Coordinate Descent, Sub-Gradient of the L1 Penalty and Soft Thresholding in Python • Coordinate descent in the least squares problem: ${\displaystyle {\frac {\partial }{\partial w_{j}}}RSS(w)=-2\rho _{j}+2w_{j}}$; i.e. ${\displaystyle {\hat {w}}_{j}=\rho _{j}}$. • Coordinate descent in the Lasso problem (for normalized features): ${\displaystyle {\hat {w}}_{j}={\begin{cases}\rho _{j}+\lambda /2,&{\text{if }}\rho _{j}<-\lambda /2\\0,&{\text{if }}-\lambda /2\leq \rho _{j}\leq \lambda /2\\\rho _{j}-\lambda /2,&{\text{if }}\rho _{j}>\lambda /2\end{cases}}}$ • Choosing ${\displaystyle \lambda }$ via cross validation tends to favor less sparse solutions and thus smaller ${\displaystyle \lambda }$ then optimal choice for feature selection. See "Machine learning: a probabilistic perspective", Murphy 2012. • Classical: Least angle regression (LARS) Efron et al 2004. • Alternating Direction Method of Multipliers (ADMM). Boyd, 2011. “Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.” Foundations and Trends in Machine Learning. Vol. 3, No. 1, 2010, pp. 1–122. • If some variables in design matrix are correlated, then LASSO is convex or not? • Tibshirani. Regression shrinkage and selection via the lasso (free). JRSS B 1996. • Convex Optimization in R by Koenker & Mizera 2014. • Pathwise coordinate optimization by Friedman et al 2007. • Statistical learning with sparsity: the Lasso and generalizations T. Hastie, R. Tibshirani, and M. Wainwright, 2015 (book) • Element of Statistical Learning (book) • https://youtu.be/A5I1G1MfUmA StatsLearning Lect8h 110913 • Fu's (1998) shooting algorithm for Lasso (mentioned in the history of coordinate descent) and Zhang & Lu's (2007) modified shooting algorithm for adaptive Lasso. • Machine Learning: a Probabilistic Perspective Choosing ${\displaystyle \lambda }$ via cross validation tends to favor less sparse solutions and thus smaller ${\displaystyle \lambda }$ than optimal choice for feature selection. ### Quadratic programming ### Constrained optimization Jaya Package. Jaya Algorithm is a gradient-free optimization algorithm. It can be used for Maximization or Minimization of a function for solving both constrained and unconstrained optimization problems. It does not contain any hyperparameters. ### Highly correlated covariates 1. Elastic net 2. Group lasso ### Other Lasso ## Comparison by plotting If we are running simulation, we can use the DALEX package to visualize the fitting result from different machine learning methods and the true model. See http://smarterpoland.pl/index.php/2018/05/ml-models-what-they-cant-learn. ## UMAP # Imbalanced Classification # Deep Learning ## Tensor Flow (tensorflow package) ## Biological applications ## Machine learning resources # Randomization inference # Bootstrap ## Nonparametric bootstrap This is the most common bootstrap method The upstrap Crainiceanu & Crainiceanu, Biostatistics 2018 ## Parametric bootstrap # Cross Validation R packages: ## Difference between CV & bootstrapping • CV tends to be less biased but K-fold CV has fairly large variance. • Bootstrapping tends to drastically reduce the variance but gives more biased results (they tend to be pessimistic). • The 632 and 632+ rules methods have been adapted to deal with the bootstrap bias • Repeated CV does K-fold several times and averages the results similar to regular K-fold ## .632 and .632+ bootstrap ${\displaystyle Err_{.632}=0.368{\overline {err}}+0.632Err_{boot(1)}}$ • Bootstrap, 0.632 Bootstrap, 0.632+ Bootstrap from Encyclopedia of Systems Biology by Springer. • bootpred() from bootstrap function. • The .632 bootstrap estimate can be extended to statistics other than prediction error. See the paper Issues in developing multivariable molecular signatures for guiding clinical care decisions by Sachs. Source code. Let ${\displaystyle \phi }$ be a performance metric, ${\displaystyle S_{b}}$ a sample of size n from a bootstrap, ${\displaystyle S_{-b}}$ subset of ${\displaystyle S}$ that is disjoint from ${\displaystyle S_{b}}$; test set. ${\displaystyle {\hat {E}}^{*}[\phi _{\mathcal {F}}(S)]=.368{\hat {E}}[\phi _{f}(S)]+0.632{\hat {E}}[\phi _{f_{b}}(S_{-b})]}$ where ${\displaystyle {\hat {E}}[\phi _{f}(S)]}$ is the naive estimate of ${\displaystyle \phi _{f}}$ using the entire dataset. ## Create partitions set.seed(), sample.split(),createDataPartition(), and createFolds() functions from the caret package. n <- 42; nfold <- 5 # unequal partition folds <- split(sample(1:n), rep(1:nfold, length = n)) sapply(folds, length)  ## Nested resampling Nested resampling is need when we want to tuning a model by using a grid search. The default settings of a model are likely not optimal for each data set out. So an inner CV has to be performed with the aim to find the best parameter set of a learner for each fold. See a diagram at https://i.stack.imgur.com/vh1sZ.png In BRB-ArrayTools -> class prediction with multiple methods, the alpha (significant level of threshold used for gene selection, 2nd option in individual genes) can be viewed as a tuning parameter for the development of a classifier. ## Pre-validation • Pre-validation and inference in microarrays Tibshirani and Efron, Statistical Applications in Genetics and Molecular Biology, 2002. • http://www.stat.columbia.edu/~tzheng/teaching/genetics/papers/tib_efron.pdf#page=5. In each CV, we compute the estimate of the response. This estimate of the response will serve as a new predictor (pre-validated predictor) in the final fitting model. • P1101 of Sachs 2016. With pre-validation, instead of computing the statistic ${\displaystyle \phi }$ for each of the held-out subsets (${\displaystyle S_{-b}}$ for the bootstrap or ${\displaystyle S_{k}}$ for cross-validation), the fitted signature ${\displaystyle {\hat {f}}(X_{i})}$ is estimated for ${\displaystyle X_{i}\in S_{-b}}$ where ${\displaystyle {\hat {f}}}$ is estimated using ${\displaystyle S_{b}}$. This process is repeated to obtain a set of pre-validated signature estimates ${\displaystyle {\hat {f}}}$. Then an association measure ${\displaystyle \phi }$ can be calculated using the pre-validated signature estimates and the true outcomes ${\displaystyle Y_{i},i=1,\ldots ,n}$. • In CV, left-out samples = hold-out cases = test set ## Custom cross validation ## Cross validation vs regularization ## Cross-validation with confidence (CVC) JASA 2019 by Jing Lei, pdf, code # Clustering See Clustering. # Mixed Effect Model # Model selection criteria ## Akaike information criterion/AIC ${\displaystyle \mathrm {AIC} \,=\,2k-2\ln({\hat {L}})}$, where k be the number of estimated parameters in the model. • Smaller is better • Akaike proposed to approximate the expectation of the cross-validated log likelihood ${\displaystyle E_{test}E_{train}[logL(x_{test}|{\hat {\beta }}_{train})]}$ by ${\displaystyle logL(x_{train}|{\hat {\beta }}_{train})-k}$. • Leave-one-out cross-validation is asymptotically equivalent to AIC, for ordinary linear regression models. • AIC can be used to compare two models even if they are not hierarchically nested. • AIC() from the stats package. ## BIC ${\displaystyle \mathrm {BIC} \,=\,\ln(n)\cdot 2k-2\ln({\hat {L}})}$, where k be the number of estimated parameters in the model. ## Overfitting ## AIC vs AUC Roughly speaking: • AIC is telling you how good your model fits for a specific mis-classification cost. • AUC is telling you how good your model would work, on average, across all mis-classification costs. Frank Harrell: AUC (C-index) has the advantage of measuring the concordance probability as you stated, aside from cost/utility considerations. To me the bottom line is the AUC should be used to describe discrimination of one model, not to compare 2 models. For comparison we need to use the most powerful measure: deviance and those things derived from deviance: generalized 𝑅2 and AIC. ## Variable selection and model estimation • training observations to perform all aspects of model-fitting—including variable selection • make use of the full data set in order to obtain more accurate coefficient estimates (This statement is arguable) # Entropy ## Definition Entropy is defined by -log2(p) where p is a probability. Higher entropy represents higher unpredictable of an event. Some examples: • Fair 2-side die: Entropy = -.5*log2(.5) - .5*log2(.5) = 1. • Fair 6-side die: Entropy = -6*1/6*log2(1/6) = 2.58 • Weighted 6-side die: Consider pi=.1 for i=1,..,5 and p6=.5. Entropy = -5*.1*log2(.1) - .5*log2(.5) = 2.16 (less unpredictable than a fair 6-side die). ## Use When entropy was applied to the variable selection, we want to select a class variable which gives a largest entropy difference between without any class variable (compute entropy using response only) and with that class variable (entropy is computed by adding entropy in each class level) because this variable is most discriminative and it gives most information gain. For example, • entropy (without any class)=.94, • entropy(var 1) = .69, • entropy(var 2)=.91, • entropy(var 3)=.725. We will choose variable 1 since it gives the largest gain (.94 - .69) compared to the other variables (.94 -.91, .94 -.725). Why is picking the attribute with the most information gain beneficial? It reduces entropy, which increases predictability. A decrease in entropy signifies an decrease in unpredictability, which also means an increase in predictability. Consider a split of a continuous variable. Where should we cut the continuous variable to create a binary partition with the highest gain? Suppose cut point c1 creates an entropy .9 and another cut point c2 creates an entropy .1. We should choose c2. ## Related In addition to information gain, gini (dʒiːni) index is another metric used in decision tree. See wikipedia page about decision tree learning. # Ensembles ## Bagging Draw N bootstrap samples and summary the results (averaging for regression problem, majority vote for classification problem). Decrease variance without changing bias. Not help much with underfit or high bias models. ### Random forest Variance importance: if you scramble the values of a variable, and the accuracy of your tree does not change much, then the variable is not very important. Why is it useful to compute variance importance? So the model's predictions are easier to interpret (not improve the prediction performance). Random forest has advantages of easier to run in parallel and suitable for small n large p problems. Random forest versus logistic regression: a large-scale benchmark experiment by Raphael Couronné, BMC Bioinformatics 2018 Arborist: Parallelized, Extensible Random Forests ## Boosting Instead of selecting data points randomly with the boostrap, it favors the misclassified points. Algorithm: • Initialize the weights • Repeat • resample with respect to weights • retrain the model • recompute weights Since boosting requires computation in iterative and bagging can be run in parallel, bagging has an advantage over boosting when the data is very large. ## Time series # p-values ## p-values ## Distribution of p values in medical abstracts ## nominal p-value and Empirical p-values • Nominal p-values are based on asymptotic null distributions • Empirical p-values are computed from simulations/permutations ## (nominal) alpha level Conventional methodology for statistical testing is, in advance of undertaking the test, to set a NOMINAL ALPHA CRITERION LEVEL (often 0.05). The outcome is classified as showing STATISTICAL SIGNIFICANCE if the actual ALPHA (probability of the outcome under the null hypothesis) is no greater than this NOMINAL ALPHA CRITERION LEVEL. ## Normality assumption # T-statistic See T-statistic. # ANOVA See ANOVA. # Goodness of fit ## Chi-square tests ## Fitting distribution # Contingency Tables ## Odds ratio and Risk ratio The ratio of the odds of an event occurring in one group to the odds of it occurring in another group  drawn | not drawn | ------------------------------------- white | A | B | Wh ------------------------------------- black | C | D | Bk  • Odds Ratio = (A / C) / (B / D) = (AD) / (BC) • Risk Ratio = (A / Wh) / (C / Bk) ## Hypergeometric, One-tailed Fisher exact test  drawn | not drawn | ------------------------------------- white | x | | m ------------------------------------- black | k-x | | n ------------------------------------- | k | | m+n  For example, k=100, m=100, m+n=1000, > 1 - phyper(10, 100, 10^3-100, 100, log.p=F) [1] 0.4160339 > a <- dhyper(0:10, 100, 10^3-100, 100) > cumsum(rev(a)) [1] 1.566158e-140 1.409558e-135 3.136408e-131 3.067025e-127 1.668004e-123 5.739613e-120 1.355765e-116 [8] 2.325536e-113 3.018276e-110 3.058586e-107 2.480543e-104 1.642534e-101 9.027724e-99 4.175767e-96 [15] 1.644702e-93 5.572070e-91 1.638079e-88 4.210963e-86 9.530281e-84 1.910424e-81 3.410345e-79 [22] 5.447786e-77 7.821658e-75 1.013356e-72 1.189000e-70 1.267638e-68 1.231736e-66 1.093852e-64 [29] 8.900857e-63 6.652193e-61 4.576232e-59 2.903632e-57 1.702481e-55 9.240350e-54 4.650130e-52 [36] 2.173043e-50 9.442985e-49 3.820823e-47 1.441257e-45 5.074077e-44 1.669028e-42 5.134399e-41 [43] 1.478542e-39 3.989016e-38 1.009089e-36 2.395206e-35 5.338260e-34 1.117816e-32 2.200410e-31 [50] 4.074043e-30 7.098105e-29 1.164233e-27 1.798390e-26 2.617103e-25 3.589044e-24 4.639451e-23 [57] 5.654244e-22 6.497925e-21 7.042397e-20 7.198582e-19 6.940175e-18 6.310859e-17 5.412268e-16 [64] 4.377256e-15 3.338067e-14 2.399811e-13 1.626091e-12 1.038184e-11 6.243346e-11 3.535115e-10 [71] 1.883810e-09 9.442711e-09 4.449741e-08 1.970041e-07 8.188671e-07 3.193112e-06 1.167109e-05 [78] 3.994913e-05 1.279299e-04 3.828641e-04 1.069633e-03 2.786293e-03 6.759071e-03 1.525017e-02 [85] 3.196401e-02 6.216690e-02 1.120899e-01 1.872547e-01 2.898395e-01 4.160339e-01 5.550192e-01 [92] 6.909666e-01 8.079129e-01 8.953150e-01 9.511926e-01 9.811343e-01 9.942110e-01 9.986807e-01 [99] 9.998018e-01 9.999853e-01 1.000000e+00 # Density plot plot(0:100, dhyper(0:100, 100, 10^3-100, 100), type='h')  Moreover,  1 - phyper(q=10, m, n, k) = 1 - sum_{x=0}^{x=10} phyper(x, m, n, k) = 1 - sum(a[1:11]) # R's index starts from 1.  Another example is the data from the functional annotation tool in DAVID.  | gene list | not gene list | ------------------------------------------------------- pathway | 3 (q) | | 40 (m) ------------------------------------------------------- not in pathway | 297 | | 29960 (n) ------------------------------------------------------- | 300 (k) | | 30000  The one-tailed p-value from the hypergeometric test is calculated as 1 - phyper(3-1, 40, 29960, 300) = 0.0074. ## Fisher's exact test Following the above example from the DAVID website, the following R command calculates the Fisher exact test for independence in 2x2 contingency tables. > fisher.test(matrix(c(3, 40, 297, 29960), nr=2)) # alternative = "two.sided" by default Fisher's Exact Test for Count Data data: matrix(c(3, 40, 297, 29960), nr = 2) p-value = 0.008853 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 1.488738 23.966741 sample estimates: odds ratio 7.564602 > fisher.test(matrix(c(3, 40, 297, 29960), nr=2), alternative="greater") Fisher's Exact Test for Count Data data: matrix(c(3, 40, 297, 29960), nr = 2) p-value = 0.008853 alternative hypothesis: true odds ratio is greater than 1 95 percent confidence interval: 1.973 Inf sample estimates: odds ratio 7.564602 > fisher.test(matrix(c(3, 40, 297, 29960), nr=2), alternative="less") Fisher's Exact Test for Count Data data: matrix(c(3, 40, 297, 29960), nr = 2) p-value = 0.9991 alternative hypothesis: true odds ratio is less than 1 95 percent confidence interval: 0.00000 20.90259 sample estimates: odds ratio 7.564602  From the documentation of fisher.test Usage: fisher.test(x, y = NULL, workspace = 200000, hybrid = FALSE, control = list(), or = 1, alternative = "two.sided", conf.int = TRUE, conf.level = 0.95, simulate.p.value = FALSE, B = 2000)  • For 2 by 2 cases, p-values are obtained directly using the (central or non-central) hypergeometric distribution. • For 2 by 2 tables, the null of conditional independence is equivalent to the hypothesis that the odds ratio equals one. • The alternative for a one-sided test is based on the odds ratio, so ‘alternative = "greater"’ is a test of the odds ratio being bigger than ‘or’. • Two-sided tests are based on the probabilities of the tables, and take as ‘more extreme’ all tables with probabilities less than or equal to that of the observed table, the p-value being the sum of such probabilities. ## Chi-square independence test The result of Fisher exact test and chi-square test can be quite different. # https://myweb.uiowa.edu/pbreheny/7210/f15/notes/9-24.pdf#page=4 R> Job <- matrix(c(16,48,67,21,0,19,53,88), nr=2, byrow=T) R> dimnames(Job) <- list(A=letters[1:2],B=letters[1:4]) R> fisher.test(Job) Fisher's Exact Test for Count Data data: Job p-value < 2.2e-16 alternative hypothesis: two.sided R> chisq.test(c(16,48,67,21), c(0,19,53,88)) Pearson's Chi-squared test data: c(16, 48, 67, 21) and c(0, 19, 53, 88) X-squared = 12, df = 9, p-value = 0.2133 Warning message: In chisq.test(c(16, 48, 67, 21), c(0, 19, 53, 88)) : Chi-squared approximation may be incorrect  ## GSEA Determines whether an a priori defined set of genes shows statistically significant, concordant differences between two biological states Two categories of GSEA procedures: • Competitive: compare genes in the test set relative to all other genes. • Self-contained: whether the gene-set is more DE than one were to expect under the null of no association between two phenotype conditions (without reference to other genes in the genome). For example the method by Jiang & Gentleman Bioinformatics 2007 ## McNemar’s test on paired nominal data # Case control study # Confidence vs Credibility Intervals # Power analysis/Sample Size determination See Power. # Common covariance/correlation structures See psu.edu. Assume covariance ${\displaystyle \Sigma =(\sigma _{ij})_{p\times p}}$ • Diagonal structure: ${\displaystyle \sigma _{ij}=0}$ if ${\displaystyle i\neq j}$. • Compound symmetry: ${\displaystyle \sigma _{ij}=\rho }$ if ${\displaystyle i\neq j}$. • First-order autoregressive AR(1) structure: ${\displaystyle \sigma _{ij}=\rho ^{|i-j|}}$. rho <- .8 p <- 5 blockMat <- rho ^ abs(matrix(1:p, p, p, byrow=T) - matrix(1:p, p, p))  • Banded matrix: ${\displaystyle \sigma _{ii}=1,\sigma _{i,i+1}=\sigma _{i+1,i}\neq 0,\sigma _{i,i+2}=\sigma _{i+2,i}\neq 0}$ and ${\displaystyle \sigma _{ij}=0}$ for ${\displaystyle |i-j|\geq 3}$. • Spatial Power • Unstructured Covariance • Toeplitz structure To create blocks of correlation matrix, use the "%x%" operator. See kronecker(). covMat <- diag(n.blocks) %x% blockMat  # Counter/Special Examples ## Correlated does not imply independence Suppose X is a normally-distributed random variable with zero mean. Let Y = X^2. Clearly X and Y are not independent: if you know X, you also know Y. And if you know Y, you know the absolute value of X. The covariance of X and Y is  Cov(X,Y) = E(XY) - E(X)E(Y) = E(X^3) - 0*E(Y) = E(X^3) = 0,  because the distribution of X is symmetric around zero. Thus the correlation r(X,Y) = Cov(X,Y)/Sqrt[Var(X)Var(Y)] = 0, and we have a situation where the variables are not independent, yet have (linear) correlation r(X,Y) = 0. This example shows how a linear correlation coefficient does not encapsulate anything about the quadratic dependence of Y upon X. ## Spearman vs Pearson correlation Pearson benchmarks linear relationship, Spearman benchmarks monotonic relationship. https://stats.stackexchange.com/questions/8071/how-to-choose-between-pearson-and-spearman-correlation x=(1:100); y=exp(x); cor(x,y, method='spearman') # 1 cor(x,y, method='pearson') # .25  ## Spearman vs Wilcoxon • Wilcoxon used to compare categorical versus non-normal continuous variable • Spearman's rho used to compare two continuous (including ordinal) variables that one or both aren't normally distributed ## Spearman vs Kendall correlation • Kendall's tau coefficient (after the Greek letter τ), is a statistic used to measure the ordinal association between two measured quantities. • Kendall Tau or Spearman's rho? ## Anscombe quartet Four datasets have almost same properties: same mean in X, same mean in Y, same variance in X, (almost) same variance in Y, same correlation in X and Y, same linear regression. ## The real meaning of spurious correlations library(ggplot2) set.seed(123) spurious_data <- data.frame(x = rnorm(500, 10, 1), y = rnorm(500, 10, 1), z = rnorm(500, 30, 3)) cor(spurious_datax, spurious_data$y) # [1] -0.05943856 spurious_data %>% ggplot(aes(x, y)) + geom_point(alpha = 0.3) + theme_bw() + labs(title = "Plot of y versus x for 500 observations with N(10, 1)") cor(spurious_data$x / spurious_data$z, spurious_data$y / spurious_data$z) # [1] 0.4517972 spurious_data %>% ggplot(aes(x/z, y/z)) + geom_point(aes(color = z), alpha = 0.5) + theme_bw() + geom_smooth(method = "lm") + scale_color_gradientn(colours = c("red", "white", "blue")) + labs(title = "Plot of y/z versus x/z for 500 observations with x,y N(10, 1); z N(30, 3)") spurious_data$z <- rnorm(500, 30, 6)
cor(spurious_data$x / spurious_data$z, spurious_data$y / spurious_data$z)
# [1] 0.8424597
spurious_data %>% ggplot(aes(x/z, y/z)) + geom_point(aes(color = z), alpha = 0.5) +
theme_bw() + geom_smooth(method = "lm") +
scale_color_gradientn(colours = c("red", "white", "blue")) +
labs(title = "Plot of y/z versus x/z for 500 observations with x,y N(10, 1); z N(30, 6)")