1A: Overview

1B: Assumptions for linear regression

The following conditions should be true in a scatterplot for a line to be considered a reasonable approximation to the relationship in the plot and for the application of the methods of inference (discussed later):

• Linearity: the data show a linear trend.
• Constant variability: the variability of the response variable about the line remains roughly constant as the predictor variable changes.
• Independent observations: the (x,y) pairs are independent; i.e., values of one pair provide no information about values of other pairs.
• Approximate normality of residuals: definition coming later…

• The variability in volume is noticeably less for smaller values of body size than for larger values.

1C: Residuals in linear regression

Our overarching goal is to fit a ‘good’ line to this data cloud. Suppose we know what the ‘best’ line that fits this data is. This means, for any given value of predictor \(x\) we obtain a predicted \(\hat{y}.\) In particular, for each observed data \((x_i, y_i)\) with \(y_i\) being the observed response, we get a predicted response \(\hat{y}_i.\)

• The vertical distance between a point in the scatterplot and the predicted value on the regression line is the residual for the point. For an observation \((x_i, y_i)\), where \(\hat{y}_i\) is the predicted value given by \[\hat{y} = b_0 + b_1 x,\] the residual is the value \[e_i = y_i - \hat{y}_i.\]

1D: Least squares regression

The least squares regression line is the line which minimizes the sum of the squared residuals for all the points in the plot.

In other words, the least squares line is the line with coefficients \(b_0\) and \(b_1\) such that the quantity \[e_1^2 + e_2^2 + \ldots + e_n^2\] is the smallest posssible value it can take.

1E: Statistical model for least squares regression

For a general population of ordered pairs \((x, y)\), the population regression model is \[Y = \beta_0 + \beta_1X + \epsilon,\] where \(\epsilon \sim N(0, \sigma)\) is the error term.

Since \(\mathrm{E}(\epsilon) = 0\), the population regression model may be rewritten in terms of conditional behaviour of \(Y\) given \(X = x\), i.e., \[\mathrm{E}(Y|X=x) = \beta_0 + \beta_1x.\]

The terms \(\beta_0\) and \(\beta_1\) are parameters - with \(b_0\) and \(b_1\) serving as estimates.

We will use R to (a) obtain estimates \(b_0\) and \(b_1\), (b) check validity of assumptions of linear regression.

1F: lm() for categorical predictors with two levels

Although the response variable in linear regression is necessarily numerical, the predictor may be either numerical or categorical. Simple linear regression only allows for categorical predictor variables with two levels. Examining categorical predictors with more than two levels requires multiple linear regression.

Fitting a simple linear regression model with a two-level categorical predictor is analogous to comparing the means of two groups, where the groups are defined by the categorical variable.

Here, we examine if there are any gender-based differences in RFFT scores in the PREVEND dataset. First, we compare the gender-stratified distribution of RFFT scores by means of violin plots.

We compare the group means as follows

prevend.samp %>% 
  rowwise() %>% 
  mutate(Gender = factor(ifelse(Gender == 1, "Male", "Female"))) %>% 
  group_by(Gender) %>% ## stratifying by gender
  summarise(m = mean(RFFT), ## strata-specific means
            q1 = quantile(RFFT, 0.25), ## strata-specific quantiles 
            q2 = quantile(RFFT, 0.5),
            q3 = quantile(RFFT, 0.75))

## # A tibble: 2 × 5
##   Gender     m    q1    q2    q3
##   <fct>  <dbl> <dbl> <dbl> <dbl>
## 1 Male    69.2  48.5    66    87
## 2 Female  67.7  44      68    88

and fit a linear model of RFFT by gender

lm(RFFT ~ Gender, data = data)

## 
## Call:
## lm(formula = RFFT ~ Gender, data = data)
## 
## Coefficients:
##  (Intercept)  GenderFemale  
##       69.238        -1.578

Note that the estimated intercept is the group mean for one category (called the baseline) - here it is the Gender = Male group. The estimated slope is the difference of the group means.

1G: The quantity \(R^2\)

1H: Hypothesis testing in regression

Inference in a regression context is usually about the slope parameter \(\beta_1\). The null hypothesis is most commonly a hypothesis of ‘no association’, which may be formulated mathematically as \[H_0: \beta_1 = 0 \text{ [denoting X and Y are NOT associated]}.\] The alternative hypothesis is given by \[H_0: \beta_1 \neq 0 \text{ [denoting X and Y are associated]}.\] We use estimator \(b_1\) in the statistic used to test whether hypothesis \(H_0\) is true. The test statistic is given by \[t = \frac{b_1 - \beta_1^0}{\text{s.e.}(b_1)} = \frac{b_1}{\text{s.e.}(b_1)},\] where \(\beta_1^0\) equals zero when the null hypothesis is true.

The \(t-\)statistic defined above follows a \(t-\)distribution with degrees of freedom \(n − 2\), where \(n\) is the number of ordered pairs \((x_i, y_i)\) in the dataset.

In order to test whether \(H_0\) holds or not, we construct the \(95\%\) confidence interval associated with \(\beta_1\) and investigate if the resultant confidence interval contains zero or not. The \(95\%\) C.I. is given by \[b_1 \pm \left(t^* \times \text{s.e.}(b_1) \right),\] where \(t^*\) is the 97.5-th percentile of a t-distribution with (n-2) degrees of freedom. (may be computed in R as qt(p = 0.975, df = n-2)).

Lab 1: R for simple linear regression

Notes for review may be found here.

Exercise 1 of 1: Using lm() to analyse data from the PREVEND study

install.packages("librarian")
librarian::shelf(tidyverse, patchwork, ggsci)

nhanes.samp.adult.500 <- read_csv("https://raw.githubusercontent.com/soumikp/bdsi_2022/main/data/nhanes.samp.adult.500.csv")
prevend.samp <- read_csv("https://raw.githubusercontent.com/soumikp/bdsi_2022/main/data/prevend.samp.csv")

2A: Introduction to multiple linear regression

In most practical settings, more than one explanatory variable is likely to be associated with a response.

Multiple linear regression is used to estimate the linear relationship between a response variable \(Y\) and several predictors \(x_1, x_2,\ldots, x_p\), where \(p\) is the number of predictors.

We extend the statistical model from just one predictor to \(p-\)many predictors as follows \[\mathrm{E}(Y|x_1, x_2, \ldots, x_p) = \beta_0 +\beta_1x_1 +\beta_2x_2 +···+\beta_px_p.\]

There are several applications of multiple regression. We will focus on two areas:

• Estimating an association between a response variable and primary predictor of interest while adjusting for possible confounding variables.
• Constructing a model that effectively explains the observed variation in the response variable.

We will again make use of the prevend.sample dataset in this module, but this time we’ll consider a more complicated example concerning the effect of statin use on cognitive function (through RFFT scores).

Statin use and cognitive function.

Statins are a class of drugs widely used to lower cholesterol. If followed, recent guidelines for prescribing statins would lead to statin use in almost half of Americans between 40 - 75 years of age and nearly all men over 60.

A few small studies have suggested that statins may be associated with lower cognitive ability.

The PREVEND study collected data on statin use as well as other demographic factors.

2B: Age, statin use and cognitive function

2C: Interpreting coefficients in multiple linear regression

The statistical model for multiple regression is based on \[\mathrm{E}(Y|x_1, x_2, \ldots, x_p) = \beta_0 + \beta_1 x_1 + \ldots + \beta_p x_p,\] where \(p\) is the number of predictors.

The coefficient \(\beta_j\) of a predictor \(X_j\) is estimated by \(b_j\), say. It iss the predicted mean change in the response corresponding to a one unit change in \(x_j\) , when the values of all other predictors remain constant.

The practical interpretation is that a coefficient in multiple regression estimates the association between a response and that predictor, after adjusting for the other predictors in the model.

Assumptions for multiple regression:

Similar to those of simple linear regression…

1. Linearity: For each predictor variable \(X_j\), change in the predictor is linearly related to change in the response variable when the value of all other predictors is held constant.
2. Constant variability: The residuals have approximately constant variance.
3. Independent observations: Each set of observations \((y, x_1, x_2, \ldots, x_p)\) is independent.
4. Normality of residuals: The residuals are approximately normally distributed.

It is not possible to make a scatterplot of a response against several simultaneous predictors :(

2D: Model diagnostics and assessment

It is not possible to make a scatterplot of a response against several simultaneous predictors.

• To assess linearity: use a modified residual plot.
  • For each (numerical) predictor, plot the residuals on the y-axis and the predictor values on the x-axis.
  • Patterns/curvature are indicative of non-linearity.
• To assess constant variablity: use the same approach as for simple regression
  • For each (numerical) predictor, plot the residuals on the y-axis and the predicted values on the x-axis.
• To assess normality of residuals: use the same approach as for simple regression
  • construct and study Q-Q plot for the residuals.

Adjusted \(R^2\) as a tool for model assessment

More information is always nice: As variables are added, \(R^2\) always increases.

But at what cost?: As variables are added, model complexity increases.

We consider using adjusted \(R^2\) (\(R^2_{\text{adj}}\)) as a tool for model assesment.

\[R^2_{\text{adj}} = 1 - (1 - R^2)\times \frac{n-1}{n-1-p}.\]

It is often used to balance predictive ability with model complexity. \(R^2_{\text{adj}}\) incorporates a penalty for including predictors that do not contribute much towards explaining observed variation in the response variable.

Unlike \(R^2\), \(R^2_{\text{adj}}\) does not have an inherent interpretation.

2E: Testing hypotheses for multiple linear regression models

Testing hypothesss about slope coefficients: \(t-\)tests.

Typically, the hypotheses of interest are \(H^k_{0}: \beta_k = 0\) (\(X_k\) and \(Y\) are NOT associated) with alternative \(H^k_A: \beta_k \neq 0\) (\(X_k\) and \(Y\) are associated).

The test statistic is again a t-statistic under the null, with degrees of freedom \(= n - p - 1\), given by \[t_k = \frac{b_k - \beta_k^0}{\text{s.e.}(b_k)},\] with a \(95\%\) confidence interval for \(\beta_k\) given by \[b_k \pm \left(t^* \times \text{s.e.}(b_k) \right),\] where \(t^*\) is the \(97.5-\)th percentile for a \(t-\)distribution with \(n-p-1\) degrees of freedom.

Testing overall goodness of model: \(F-\)tests.

The \(F-\)statistic is used in an overall test of the model to assess whether the predictors in the model, considered as a group, are associated with the response.

\[H_0 :\beta_1 = \beta_2 =\ldots=\beta_p = 0 \text{ vs } H_A: H_0 \text{ is false,}\] where \(H_0\) is false if and only if at least one of the slope coefficients \(\beta_k\) is not 0

There is sufficient evidence to reject \(H_0\) if the \(p-\)value of the \(F-\)statistic is smaller than or equal to \(\alpha.\)

Again, R does the hard work for us: the F-statistic and its associated p-value are displayed in the output.

2F: Interaction in regression

The multiple regression model assumes that when one of the predictors \(x_j\) is changed by 1 unit and the values of the other variables remain constant, the predicted response changes by \(\beta_j\), regardless of the values of the other variables.

A statistical interaction occurs when this assumption is not true, such that the effect of one explanatory variable \(x_j\) on the response depends on the particular value(s) of one or more other explanatory variables.

As an example, we specifically examine interaction in a two-variable setting, where one of the predictors is categorical and the other is numerical.

Consider data on total cholesterol level (mmol/L) from age (yrs.) and diabetes status in a dataset of size \(n = 473\).

We fit the linear model \[\mathrm{E}(\text{Total cholesterol|Age, Diabetes status}) = \beta_0 + \beta_1 \text{Age} + \beta_2 \mathrm{I}(\text{has diabetes}).\] This model yields two fitted lines, one for individuals with diabetes and one for individuals without diabetes - we overlay these fitted lines on the scatterplot.

Next, we consider two separate models for the relationship between total cholesterol and age; one in diabetic individuals and one in non-diabetic individuals.

\[\mathrm{E}(\text{Total cholesterol|Age, Diabetes status}) = \beta_0 + \beta_1 \text{ Age } + \beta_2 \ \mathrm{I}(\text{has diabetes}) + \beta_3 \text{ Age } \times \ \mathrm{I}(\text{has diabetes}).\] The estimated model coefficients are given by

##     (Intercept)             Age     DiabetesYes Age:DiabetesYes 
##     4.695702513     0.009638183     1.718704342    -0.033451562

\[\widehat{TotChol} = 4.70 + 0.0096 \times \text{ Age} + 1.72 \times \ \mathrm{I}(\text{has diabetes}) - 0.033 \times \text{ Age }\times \mathrm{I}(\text{has diabetes}).\]

Hence the fitted model equation for diabetics is given by \[\widehat{TotChol} = 6.42 - 0.023 \times \text{Age},\] and that for non-diabetics is given by \[\widehat{TotChol} = 4.70 + 0.0096 \times \text{Age}.\]

Lab 2: R for multiple linear regression

Notes for review may be found here.

Exercise 1 of 2: Verify validity of linear model assumptions for PREVEND study

library(oibiostat)
data(prevend.samp)

prevend.samp <- prevend.samp %>% 
  rowwise() %>% 
  mutate(Statin = factor(ifelse(Statin == 1, "Yes", "No"))) %>% 
  select(c(RFFT, Age, Statin))

#fit a multiple regression model
model1 = lm(RFFT ~ Statin + Age, data = prevend.samp)

Exercise 2 of 2: Model interaction effects for PREVEND study

#load the data
library(oibiostat)
data("prevend.samp")
#convert Statin to a factor
prevend.samp$Statin = factor(prevend.samp$Statin, levels = c(0, 1),
                             labels = c("NonUser", "User"))
#fit the model
model.RFFT.interact = lm(RFFT ~ Age*Statin, data = prevend.samp)
coef(model.RFFT.interact)

##    (Intercept)            Age     StatinUser Age:StatinUser 
##    140.2031114     -1.3149119    -13.9720216      0.2474466

\(\infty\)