David's Wiki

Testmathj

Revision as of 22:34, 7 September 2019 by David (talk | contribs)
\( \newcommand{\P}[]{\unicode{xB6}} \newcommand{\AA}[]{\unicode{x212B}} \newcommand{\empty}[]{\emptyset} \newcommand{\O}[]{\emptyset} \newcommand{\Alpha}[]{Α} \newcommand{\Beta}[]{Β} \newcommand{\Epsilon}[]{Ε} \newcommand{\Iota}[]{Ι} \newcommand{\Kappa}[]{Κ} \newcommand{\Rho}[]{Ρ} \newcommand{\Tau}[]{Τ} \newcommand{\Zeta}[]{Ζ} \newcommand{\Mu}[]{\unicode{x039C}} \newcommand{\Chi}[]{Χ} \newcommand{\Eta}[]{\unicode{x0397}} \newcommand{\Nu}[]{\unicode{x039D}} \newcommand{\Omicron}[]{\unicode{x039F}} \DeclareMathOperator{\sgn}{sgn} \def\oiint{\mathop{\vcenter{\mathchoice{\huge\unicode{x222F}\,}{\unicode{x222F}}{\unicode{x222F}}{\unicode{x222F}}}\,}\nolimits} \def\oiiint{\mathop{\vcenter{\mathchoice{\huge\unicode{x2230}\,}{\unicode{x2230}}{\unicode{x2230}}{\unicode{x2230}}}\,}\nolimits} \)

Statisticians

Statistics for biologists

http://www.nature.com/collections/qghhqm

Transform sample values to their percentiles

https://stackoverflow.com/questions/21219447/calculating-percentile-of-dataset-column 

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

File:Boxplot.svg

  • 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.

Dot-box plot

File:Boxdot.svg

Other boxplots

File:Lotsboxplot.png

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

Box-Cox transformation

the Holy Trinity (LRT, Wald, Score tests)

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)
d2$u %*% diag(d2$d) %*% t(d2$v) # equal to cmat
d2$u %*% diag(sqrt(d2$d))
# [,1]      [,2]
# [1,] -0.6322816 0.7692937
# [2,]  0.9305953 0.5226872

Linear Regression

Regression Models for Data Science in R by Brian Caffo

Comic https://xkcd.com/1725/

Different models (in R)

http://www.quantide.com/raccoon-ch-1-introduction-to-linear-models-with-r/

dummy.coef.lm() in R

Extracts coefficients in terms of the original levels of the coefficients rather than the coded variables.

model.matrix, design matrix

ExploreModelMatrix: Explore design matrices interactively with R/Shiny

Contrasts in linear regression

  • Page 147 of Modern Applied Statistics with S (4th ed)
  • https://biologyforfun.wordpress.com/2015/01/13/using-and-interpreting-different-contrasts-in-linear-models-in-r/ This explains the meanings of 'treatment', 'helmert' and 'sum' contrasts.
  • A (sort of) Complete Guide to Contrasts in R by Rose Maier 
    mat

##      constant NLvMH  NvL  MvH
## [1,]        1  -0.5  0.5  0.0
## [2,]        1  -0.5 -0.5  0.0
## [3,]        1   0.5  0.0  0.5
## [4,]        1   0.5  0.0 -0.5
mat <- mat[ , -1]

model7 <- lm(y ~ dose, data=data, contrasts=list(dose=mat) )
summary(model7)

## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  <br />
## (Intercept)  118.578      1.076 110.187  < 2e-16 ***
## doseNLvMH      3.179      2.152   1.477  0.14215  <br />
## doseNvL       -8.723      3.044  -2.866  0.00489 ** 
## doseMvH       13.232      3.044   4.347 2.84e-05 ***

# double check your contrasts
attributes(model7$qr$qr)$contrasts
## $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

https://en.mimi.hu/mathematics/exposure_variable.html

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

Dealing with heteroskedasticity; regression with robust standard errors using R

Linear regression with Map Reduce

https://freakonometrics.hypotheses.org/53269

Relationship between multiple variables

Visualizing the relationship between multiple variables

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 = T P^\mathrm{T} + E\)
\(\displaystyle Y = U Q^\mathrm{T} + F\)

where Template:Mvar is an \(\displaystyle n \times m\) matrix of predictors, Template:Mvar is an \(\displaystyle n \times p\) matrix of responses; Template:Mvar and Template:Mvar are \(\displaystyle n \times l\) matrices that are, respectively, projections of Template:Mvar (the X score, component or factor matrix) and projections of Template:Mvar (the Y scores); Template:Mvar and Template:Mvar are, respectively, \(\displaystyle m \times l\) and \(\displaystyle p \times l\) orthogonal loading matrices; and matrices Template:Mvar and Template:Mvar are the error terms, assumed to be independent and identically distributed random normal variables. The decompositions of Template:Mvar and Template:Mvar are made so as to maximise the covariance between Template:Mvar and Template:Mvar (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

https://en.wikipedia.org/wiki/Feature_selection

Independent component analysis

ICA is another dimensionality reduction method.

ICA vs PCA

ICS vs FA

Correspondence analysis

https://francoishusson.wordpress.com/2017/07/18/multiple-correspondence-analysis-with-factominer/ and the book Exploratory Multivariate Analysis by Example Using R

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.

Visualize the random effects

http://www.quantumforest.com/2012/11/more-sense-of-random-effects/

Calibration

  • How to determine calibration accuracy/uncertainty of a linear regression?
  • Linear Regression and Calibration Curves
  • Regression and calibration Shaun Burke
  • calibrate package
  • 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)         <br />
z1 = 1 + 0*x      <br />
pr1 = 1/(1+exp(-z1))
y1 = rbinom(1000,1,pr1)<br />
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)         <br />
z2 =  -2 + 0*x      <br />
pr2 = 1/(1+exp(-z2))<br />
y2 = rbinom(1000,1,pr2)<br />
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.

ROC curve

See ROC.

NRI (Net reclassification improvement)

Maximum likelihood

Difference of partial likelihood, profile likelihood and marginal likelihood

Generalized Linear Model

Lectures from a course in Simon Fraser University Statistics.

Doing magic and analyzing seasonal time series with GAM (Generalized Additive Model) in R

Link function

Link Functions versus Data Transforms

Quasi Likelihood

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

Plot

https://strengejacke.wordpress.com/2015/02/05/sjplot-package-and-related-online-manuals-updated-rstats-ggplot/

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<br />
# -14.1083   -4.2773   -0.5484    5.4838   15.2922<br />
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)  <br />
# (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 *<br />
#   TreatFT       4.5631     2.1333   2.139 0.036035 *<br />
#   ---
#   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

Saturated model

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

Simulate data from a specified density

Signal to noise ratio

\(\displaystyle \frac{\sigma^2_{signal}}{\sigma^2_{noise}} = \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).

File:Hist bh.svg

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

File:Scatterhist.svg

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.

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).

String Permutations Algorithm

https://youtu.be/nYFd7VHKyWQ

Empirical Bayes Normal Means Problem with Correlated Noise

Solving the Empirical Bayes Normal Means Problem with Correlated Noise Sun 2018

The package cashr and the source code of the paper

Bayes

Bayes factor

Empirical Bayes method

Naive Bayes classifier

Understanding Naïve Bayes Classifier Using R

MCMC

Speeding up Metropolis-Hastings with Rcpp

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

https://en.wikipedia.org/wiki/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

https://stats.stackexchange.com/questions/66586/is-there-a-test-to-determine-whether-glm-overdispersion-is-significant

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

See Survival data analysis

Logistic regression

Simulate binary data from the logistic model

https://stats.stackexchange.com/questions/46523/how-to-simulate-artificial-data-for-logistic-regression 

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

https://www.analyticsvidhya.com/blog/2015/10/basics-logistic-regression

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)))<br />

AUC

A small introduction to the ROCR package 

       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

Retrieved from "https://wiki.davidl.me/index.php?title=Testmathj&oldid=51"