MT5762 Lecture 14

Here we look to unify several of the models seen previously (which generated CIs and \( p \)-values)

The class of models are Linear Models

They have a specific structure for the signal and noise

My favourite equation

Many models are effectively regressions - they are broadly of the form:

\[ y = f(x_1,..., x_p) + noise \]

So we need to estimate the signal \( f \), which is usually contingent on an idea of the distribution of the error (because we have to tease the signal and noise apart)

As we'll see, our estimates for the signal are optimised in light of our model for noise

The term linear can be used differently in different contexts

• It might refer to some shape of relationship e.g. a straight line relating some variable \( x \) to \( y \)
• Here we mean the linear form. In short we model the signal as \( y = f(x1,...x_p) \) where it is

\[ y = \mathbf{X}\boldsymbol{\beta} \]

• Recalling from last week, this matrix equation means things like:

\[ y_i = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p \]

This is veery wide. Naturally it includes a straight line relating \( x \) and \( y \):

\[ y = \beta_0 + \beta_1 x_1 \]

or indeed just modelling \( y \) as a constant (a mean would be best)

\[ y = \beta_0 \]

Let's look at some examples

We frequently need to put categorical variables into an equation - demanding some numeric representation

We do this using dummy variable encoding/representation (comp. sci. focussed people use the ridiculous term one-hot encoding)

We cover this here - linear models encompass this

We just need to get clever with the design matrix i.e. the bit on the LHS of the regression equation expressing our covariates

• We create a binary matrix
• The 0/1 codings indicate which level of the factor variable (our categorical covariate) each observation belongs to
• Each level therefore gets a parameter in the parameter vector

So - say \( X \) has 3 levels, with the first 6 values being \( \mathbf{x} = [A, A, B, B, C, C] \), we might encode:

\[ \mathbf{x} \rightarrow \begin{bmatrix} 1 & 0 & 0\\ 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ 0 & 0 & 1\\ \end{bmatrix} \]

A “1” in the first column means level A, a “1” in the second column means level B and a “1” in the third column means level C

Through the magic of matrix multiplication:

\[ \hat{\mathbf{y}} = \mathbf{X}\boldsymbol{\hat{\beta}} = \begin{bmatrix} 1 & 0 & 0\\ 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ 0 & 0 & 1\\ \end{bmatrix} \begin{bmatrix} \hat{\beta_0}\\ \hat{\beta_1}\\ \hat{\beta_2}\\ \end{bmatrix} = \begin{bmatrix} \hat{\beta_0}\\ \hat{\beta_0}\\ \hat{\beta_1}\\ \hat{\beta_1}\\ \hat{\beta_2}\\ \hat{\beta_2}\\ \end{bmatrix} \]

So, with an appropriate fitting objective, the estimated \( \hat{\beta_j} \) could be means of each factor level (A, B, C).

Armed with this idea of constructing design matrices to give complex models, we can combine any number of categorical and numeric covariates

Through the magic of matrix multiplication:

\[ \hat{\mathbf{y}} = \mathbf{X}\boldsymbol{\hat{\beta}} = \begin{bmatrix} 1 & 0 & 0 & x_{1,2}\\ 1 & 0 & 0 & x_{2,2}\\ 1 & 1 & 0 & x_{3,2}\\ 1 & 1 & 0 & x_{4,2}\\ 1 & 0 & 1 & x_{5,2}\\ 1 & 0 & 1 & x_{6,2}\\ \end{bmatrix} \begin{bmatrix} \hat{\beta_0}\\ \hat{\beta_1}\\ \hat{\beta_2}\\ \hat{\beta_3}\\ \end{bmatrix} = \begin{bmatrix} \hat{\beta_0} + \hat{\beta_3}x_{1,2}\\ \hat{\beta_0} + \hat{\beta_3}x_{2,2}\\ \hat{\beta_0} + \hat{\beta_1} + \hat{\beta_3}x_{3,2}\\ \hat{\beta_0} + \hat{\beta_1} + \hat{\beta_3}x_{4,2}\\ \hat{\beta_0} + \hat{\beta_2} + \hat{\beta_3}x_{5,2}\\ \hat{\beta_0} + \hat{\beta_2} + \hat{\beta_3}x_{6,2}\\ \end{bmatrix} \]

So, with an appropriate fitting objective, we now estimate a strightline relationship between \( x_2 \) and some adjustment for the factor level of \( x_1 \)

We construct the design matrix i.e. create dummy variables

  x_Design <- cbind(rep(1, 6), rep(c(0,1,0), c(2,2,2)), 
                    rep(c(0,0,1), c(2,2,2)), x_2)

  lmData <- data.frame(y, x_Design)  

  lmData

          y V1 V2 V3      x_2
1 11.491858  1  0  0 5.303632
2  5.528969  1  0  0 1.356557
3 17.623767  1  1  0 4.235338
4 30.908142  1  1  0 9.813859
5 34.429791  1  0  1 7.826814
6 24.410992  1  0  1 2.542402

What is this saying? Let's just take the estimated coefficients.

  coef(exampleLM)

(Intercept)        x_1B        x_1C         x_2 
   1.790517    8.300291   17.167700    2.017930 

• For starters, level A isn't listed
• It is the baseline (i.e. “contained” in the intercept)
• When x_2 is zero and you are in the baseline class(es), the predicted mean value of y is the “(Intercept)” parameter (here 1.79)

Now we know a bit more, we can look at our previous models through the prism of a linear model

The ANOVA provides:

  weedANOVA <- aov(Ti ~ Group, data = TiData)

  summary(weedANOVA)

            Df Sum Sq Mean Sq F value   Pr(>F)    
Group        2 118.60   59.30   28.31 1.43e-08 ***
Residuals   43  90.06    2.09                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  TukeyHSD(weedANOVA)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Ti ~ Group, data = TiData)

$Group
            diff        lwr      upr     p adj
nth-bhb 2.477778  0.9544691 4.001086 0.0008236
pm-bhb  3.750000  2.5402571 4.959743 0.0000000
pm-nth  1.272222 -0.1008697 2.645314 0.0743237

Using lm (which dummy-codes Group):

  coef(weedLM)

(Intercept)    Groupnth     Grouppm 
   5.100000    2.477778    3.750000 

• Observe estimate of baseline mean
• Estimates of differences between baseline and nth and pm
• Overall \( F \)-stat, df and \( p \)-values match

A nice thing about doing things as linear models, is our processes for getting estimates, errors, assumption checking etc use the same tools each time:

  qqnorm(weedLM$residuals)
  qqline(weedLM$residuals)

So, \( t \)-test, ANOVA, regression, whatever, will conduct examinations of residuals. They're all linear models.

  shapiro.test(weedLM$residuals)

    Shapiro-Wilk normality test

data:  weedLM$residuals
W = 0.98038, p-value = 0.6216

An ANalysis of COVAriance (ANCOVA) is sometimes referred to and is simply a mix of a numeric and categorical covariate for predicting numeric \( y \) (assuming Normal noise):

• It is of course, a linear model
• It may involve an interaction, which we cover later
• You should be familar with the term, even though it is just a simple case of a linear model

So now we'll turn the handle to make more complex models of signal

The model for noise (independent draws from a Normal distribution) remains the same throughout the linear models

• BAD - The target variable/response, code as 1 for a loan default and 0 if paid back.
• REASON - What the loan was for: HomeImp for home improvement and DebtCon for debt consolidation.
• JOB - a set of six occupational categories.
• MORTDUE - the amount remaining on the existing mortgage.
• VALUE - value of the current property
• DEBTINC the debt-to-income ratio i.e. if this is high, then the person owes a lot relative to their income (they're highly indebted).
• YOJ - the number of years at the present job.
• DEROG - the number of major derogatory reports
• CLNO - the number of trade lines
• DELINQ - the number of delinquent trade lines
• CLAGE - the age of the oldest trade line (months)
• NINQ - the number of recent credit checks.

Specify the linear model: \( y = f(Reason, Value) \), where the function/signal is of the linear form \( \mathbf{y} = \mathbf{X}\boldsymbol{\beta} \) and our model for noise is independent draws from a Normal distribution.

  loanData <- loanData %>% select(DEBTINC, REASON, VALUE) %>% na.omit()

  loanLM <- lm(DEBTINC ~ REASON + VALUE, data = loanData)  

  summary(loanLM)

• All estimates are statistically significant
• The intercept parameter test is not interesting (why?)
• The REASONHomeImp says Debt Consolidation loans are associated with higher Debt-to-income ratios (1 unit in fact)
• As VALUE increases by 1 unit, average DEBTINC is predicted to increase by 1.9e-5 units