TU BASICS: Regression modeling

linear regression is also called linear model:

• simple linear model: Y and X
• multiple linear model: Y and X1, X2, …, Xp
• generalized linear models (and in particular logistic regression): Y and X1, X2, …, Xp when Y is not Gaussian

Study on risks factor for low birth weights from Baystate medical center in Springfield (MA, USA) during 1986 [Hosmer & Lemeshow (1989), Applied Logistic Regression]

- low: low birth weight indicator coded as: 1 (low birth weight, i.e. ≤2.5kg) and 0 (normal birth weight)
- age: mother’s age (in years) at the time of birth
- lwt: mother’s weight (in pounds) at the time of the last menstrual period
- race: mother’s race coded as: 1 (white), 2 (black) and 3 (other)
- smoke: mother’s smoking status coded as 1 (smoker) and 0 (non-smoker)
- ptl: number of previous premature labours
- ht: hypertension history indicator coded as: 1(yes) and 0 (no)
- ui: uterine irritability indicator coded as: 1 (yes) and 0 (no)
- ftv: number of physician visits during the first trimester
- bwt: infant’s birth weight (in grams)

Are X and Y values linked ?

A few definitions:

• cor(X,Y)=Cov(X,Y)√Var(X)Var(Y)
• Var(X)=E[(X−E(X))2]
• Cov(X,Y)=E[(X−E(X))(Y−E(Y))]=E(XY)−E(X)E(Y)
• E(X)=∑xxP(X=x)

Those quantities are estimated over a sample (x1,…,xn) by:

• ˉx=1n∑ni=1xi
• s2x=∑ni=1(xi−ˉx)2n−1
• sxy=∑ni=1(xi−ˉx)(yi−ˉy)n−1
• r=sxy√s2xs2y

NB: the correlation is a scaled version of the covariance,
it is thus an association measure between −1 and 1
(the closer to 0, the weaker the link)

Here, the correlation between the mother’s age and the birth weight is 0.09

This is the linear correlation, also called Pearson correlation coefficient.

Many other measures have been proposed to explore and quantify the link between two variables:

• Spearman correlation
• Mutual Information
• Maximum Information Criterion
• …

It all started with a scatter plot…

⇒ Regression, a.k.a. linear model:
what is the optimal line going through this scatter plot ?

y=a+bx

We can write such an equation for relating yi values to xi values:

yi=β0+β1xi+εi

• β0 is interpreted as the average value of Y when xi is 0
• β1 represents the slope of the regression line
• εi is the random error separating the regression line from the observed value

The regression problem can be formalized as:

How do we characterize the conditional distribution of the random variable Y knowing the values x of the random variable X ?

Y|x∼?

• its expected value E(Y|x)
• its variance Var(Y|x)
• its probability distribution

The simple linear model can be written as:

Y|x∼N(β0+β1x,σ2)

The regression equation then writes itself as: E(Y|x)=β0+β1x

1. linearity of the relationship between Y and X
2. homoskedasticity of the observations yi (same variance across observations)
3. normality of the observations yi
4. independence of the yi across observations

When X increases by 1 unit, Y increases on average by β1 units
regardless of the inital value of X

E(Yi|Xi=a+1)−E(Yi|Xi=a)=[β0+β1(a+1)]−[β0+β1a]=β1

Var(Yi|xi)=σ2 for all xi
usually not enough data to assess numerically – instead a graphical assesment can be performed

Ex: for any mother’s age (x), the conditional birth weight (Y|x) must have the same variance

Yi|xi∼N for all xi
usually not enough data to assess – instead a graphical assessment can be performed on the marginal histogram of Y (instead of Y|x)

Ex: for any mother’s age (x), the conditional birth weight (Y|x) must have a Gaussian distribution

Cov(Yi,Yj)=0 for all i≠j
Such a hypothesis cannot be checked directly on the data, but must be assessed from what we know from the data collection process

Ex: birth weights of the 189 infants are independent, assuming there is no twins (or even siblings) among them

• E(Yi|xi)=β0+β1xi(i=1,…,n)
• Yi=β0+β1xi+εi(i=1,…,n)

where β0+β1xi is the determinist part
and εi is the random part: εi=Yi−E(Yi|xi)=Yi−(β0+β1xi)

εi are called the random errors, or the residuals

The 4 hypothesis from the linear model can be expressed with respect to the random errors:

• independence of the εi ⇒ independence of Yi
• E(εi)=0 ⇒ linearity hypothesis
• Var(εi)=σ2 ⇒ Var(Yi|xi)=σ2 (homoskedasticity)
• normality of the εi ⇒ normality of the Yi|xi.

• Linearity: E(Yi|Xi=xi)=E(Yi|xi)=β0+β1xi or E(εi)=0 (i=1,…,n)
  ✅ check: X×Y scatter plot + post-fit residuals

• Homoskedasticity: Var(Yi|xi)=σ2 or Var(εi)=σ2 for all i
  ✅ check: X×Y scatter plot + post-fit residuals

• Normality: Yi|xi∼N for any xi or εi∼N for all i
  ✅ check: histogram of Y before estimation + post-fit residuals

• Independence: cov(Yi,Yj)=0 or cov(εi,εj)=0 for all pairs i≠j
  ✅ check: context

One way to define the “Best” line is the line that will minimize the sum of squarred errors.

The Least Squares criteria minimizes: Φ(β0,β1)=n∑i=1(yi−β0−β1xi)2

Minimizing Φ(β0,β1) gives ˆβ0 and ˆβ1 (solved setting the partial derivatives with respect to β1 and β0 respectively to 0):

ˆβ1=∑ni=1(xi−ˉx)(yi−ˉy)∑ni=1(xi−ˉx)2etˆβ0=ˉy−^β1ˉx

with ˉy=1n∑ni=1yi and ˉx=1n∑ni=1xi

Another possible formula is: ˆβ1=¯xy−ˉxˉy¯x2−ˉx2
with ¯xy=1n∑ni=1xiyi and ¯x2=1n∑ni=1x2i

It can be shown that, as n (the number of observations) gets larger:

• ^β0∼N(β0,σ2^β0) ⇒ ^β0−β0s(^β0)∼Student(n−2)

• ^β1∼N(β1,σ2^β1) ⇒ ^β1−β1s(^β1)∼Student(n−2)

It follows that

• Confidence interval at 100×(1−α)% for β0 : [^β0±tn−21−α/2×s(ˆβ0)]

• Confidence interval at 100×(1−α)% for β1 : [^β1±tn−21−α/2×s(ˆβ1)]

The most important test in the simple linear model is whether ˆβ1 is significantly different from 0,
i.e. does X significantly informs the prediction of Y:

• null hypothesis H0: β1=0
  (no association between Y & X)
• alternative hypothesis H1: β1≠0
  (linear association between Y & X)

⇒ test statistics: T=ˆβ1s(ˆβ1)∼H0Student(n−2)

• CI1−α%(β1) contains 0 ⇔ non-rejection of H0 at the significance level α
• CI1−α%(β1) does NOT contain 0 ⇔ rejection of H0 at the significance level α

CI gives a idea of the precision around the estimation of β1, but not its significance level (p-value)…
One can only know if the p-value is below or above α from the CI

R2 coefficient is the square of the correlation and is related to β1

R2=β21σ2Xσ2Y

The correlation estimator can then be written as: r=^β1sXsY

📝

• 👉 write the likelihood for the simple linear model

+

• 👉 load the data from birthweight_data.txt
  
• 👉 using dplyr, compute a new data.frame with both the mean and the standard deviation (look at the group_by and summarize_all() function) of all original variables for byh avaier’s age. Do value observedes it makes sense ?
  
• 👉 add the number of observations summarized with the add_tally() function
  
• 👉 draw a scatter plot with ggplot2 of average birth weight as a function of the mother’s age

• 👉 using the lm() function (look at the help and the examples), fit a simple linear model explaining the average birth weight from the mother’s age
• 👉 explore the results with the ls() function, the $ operator and the summary() function
• 👉 assess the linear model assumpions
• 👉 interpret your results (take particular care for β0)
• 👉 fit a second simple linear model explaining the individual birth weight from the mother’s age using the non-aggregated original data
• 👉 compare the results. How can you explain the differences ?

Regardless of which regression model you are using:

1. Model specification according to the research question, and the data available
2. Estimation with point estimates and Confidence Intervals for the model parameters
3. Significance testing for each model parameters of interest
4. Model adequation diagnoses
5. Results presentation (often as a table)
6. Results interpretation and discussion

yi=β0+β1x1i+β2x2i+⋯+βpx2p+εi=β0+(p∑i=1βixi)+εi

Y=Xβ+ϵ

ˆβ=(XTX)−1XTY

❗️Requires XTX to be invertible !

This is no the case if:

• p≥n (high-dimensional case)
• XTX is not of full rank (e.g. one covariate can be computed as a perfect linear combination of the other)
  numerical instabilities when covaraites are extremely correlated

For categorical predictors, the design matrix X uses modality indicator variables:

Xj=[m1mKm1⋮mK] → Xjm1=[101⋮0], …, Xjmn=[010⋮1]

❗️Need to set a reference class to ensure identifiability

?model.matrix()

Xj×Xm

⇒ the value of Xm modifies the effect of Xj on Y

❗️if an interaction is included in a linear model, the main effects involved must always be included as well
otherwise the interaction is not interpretable

Due to numerical complexity, it is very rare to use higher order (>2) interactions…

potential confounders must be included in the model !
X2 is a counfounding factor for the relationship between X1 and Y if:

• X2 is associated with both X1 and Y
• X2 is a direct or indirect cause of Y, but not a consequence
• X2 is not on the causal path between X1 and Y

Assessthe relative variation RV=|βmodel11−βmodel21|βmodel11
⇒ arbitrary thresholds of 10% or 20% to declare potential confounding factors

The Likelihood Ratio Test (LRT) can be used to compare nested models, in particular for testing several coefficients at once.

AIC & BIC criteria:

• based on −2log-likelihood
• penalize for the number of parameters (the more parameters, the better the fit, mechanically)
• can be used to compare non-nested models on the same data
• the lower the better

AIC=2p−2log-likelihood
BIC=plog(n)−2log-likelihood
R2adj=1−(1−R2)n−1n−p−1

❗️when p becomes large, the number of coefficients tested in the model becomes large as well, and one must be careful and deal with the multiple testing issue

To be continued…

One of the recurring question in multiple regression is:

Which variables should be included in the model ?

Often dealt with stepwise method inclusion

❗️this is not statistically sound and should be avoided

Modern methods, such as the LASSO or Random Forest encompass a framework for suitable variable selection in multivariate regression models… To be continued

+

• 👉 load the data from birthweight_data.txt
• 👉 how can you best explain the individual birth weight from the other variables ?

Be able to:

• identify wether a linear or logistic regression would let you answer the scientific question of interest given the available data
• propose an adequate model for the question (model, outcome, predictors)
• fit the model in R
• interpret and discuss the results

Linear regression model:

E(Y|x1,x2,…,xp)=β0+β1x1+β2x2+⋯+βpxp

If Y is a binary variable: ∈ [0,1]:

E(Y|x1,x2,…,xp)=P(Y=1|x1,x2,…,xp)

logit(p)=log(p1−p) ∈ ]−∞,∞[

expit(x)=ex1+ex ∈ [0,1]

NB: expit(logit(p))=p et logit(expit(x))=x

logit(P(Y=1|x1,x2,…,xp))=β0+β1x1+β2x2+⋯+βpxp

⇒P(Y=1|x1,x2,…,xp)=eβ0+β1x1+β2x2+⋯+βpxp1+eβ0+β1x1+β2x2+⋯+βpxp

• the logit is called the link function:
  makes the link between the expectation of the outcome and the linear model

NB: No error term, randomness is directly encompassed into the model

• p1−p

If logit(P(Y=1|x)=β0+β1x: log(P(Y=1|X=a)1−P(Y=1|X=a))−log(P(Y=1|X=a)1−P(Y=1|X=a+1))=β1

CI(OR)1−α=[e^β1±u1−αˆs(β1)]

❗️Normal asymptotic approximation ⇒ Wald test:

^β1ˆs(β1)∼H0N(0,1)whenn→∞

One way to assess the linearity hypothesis is to split-up the continuous covariates into categories and check that the coefficients associated with each categories are coherent with a linear relationship.

This is the only hypothesis to be check in the logistic regression

+

• 👉 load the data from birthweight_data.txt
• 👉 how can you best explain the binary variable low from the other variables ? Use the glm() function with the "binomial link"
• 👉 change the reference category for race to obtain the CI for the OR for smoking in black women

• Poisson regression for discrete quantitative outcome (e.g. counts)
  link function is the log
• Probit & Tobit models as alternatives to logistic regression
  (different link functions)
• …

GLM assume that all the observations are iid

Linear mixed-effect models can be used to model grouped data

Ex: longitudinal data, batch effects, etc…

Regardless of which regression model you are using:

1. Model specification according to the research question, and the data available
2. Estimation with point estimates and Confidence Intervals for the model parameters
3. Significance testing for each model parameters of interest
4. Model fit diagnoses
5. Results presentation (often as a table)
6. Results interpretation and discussion

Points of significance | Statistics for Biologists series:

• Altman & Krzywinski (2015) Simple linear regression Nature Methods 12:999–1000
• Altman & Krzywinski (2015) Multiple linear regression Nature Methods 12:1103–1104
• Lever, Altman & Krzywinski (2015), Logistic regression Nature Methods 13:541–542
• Altman & Krzywinski (2016) Analyzing outliers: influential or nuisance ? Nature Methods 13:281–282
• Altman & Krzywinski (2016) Regression diagnosis Nature Methods 13:385–386