By the end of these workshops and your own reading, you should be able to

• Describe what the likelihood of a model is
• Using unexplained variance to assess model fit
• Describe the limitations of \(R^2\) as model fit statistic
• Use likelihood ratio test (F-test) to compare nested models

Download exercises from Week 3 folder on NOW and move them into your R-project directory.

• Precedence to simplicity: of two competing theories, the simpler explanation of an entity is to be preferred.
• We don’t assume more complex theories, if there is not sufficient evidence / reason to do so.

• Model comparisons are based on how well a model can capture data.
• In any statistical analysis, we assume our data are drawn from some probability distribution with a set of parameters.
• This is what we mean by statistical model, a probabilistic generative model of the statistical population like

\[ y_i \sim \mathcal{N}(\mu_i, \sigma^2), \mu_i = \beta_0 + \beta_1 \cdot x_i \]

• We aim to find the best model of the population in our set of candidate models.

Probability is the chance to obtain a hypothetical outcome given a particular model (hypothesis) such as a fair coin, i.e. \(\mathcal{P}(\text{data}|\theta)\); like p-values, i.e. the probability of the data given the null hypothesis.

Likelihood is the chance to observe a parameter value, model or hypothesis after having seen the data, i.e. \(\mathcal{L}(\theta|\text{data})\)

• Likelihood tells us the compatibility of model and data.
• If data are more likely under model \(\mathcal{M}_0\) than under model \(\mathcal{M}_1\), model \(\mathcal{M}_0\) is more compatible with the data.

A possible model of each observed outcome \(y_i\) is

\[ y_i \sim \mathcal{N}(\mu, \sigma^2)\\ \text{for } i \in 1\dots n, \]

We are modelling \(y\) as normally distributed with two unknowns, the mean \(\mu\) and standard deviation \(\sigma^2\).

Probability of the observed values, \(y_1, y_2, y_3\dots y_n\), assuming values for \(\mu\) and \(\sigma^2\).

\[ \mathcal{P}(y_1\dots y_n\mid\mu, \sigma^2) \]

The joint probability is their product (assuming all \(y\)s are conditionally independent of one another).

\[ \prod_\text{i=1}^n \mathcal{P}(y_i\mid \mu, \sigma^2) \]

Replace \(\mu\) and \(\sigma^2\) with their maximum likelihood estimates \(\hat\mu\) and \(\hat\sigma^2\):

\[ \prod_\text{i=1}^n \mathcal{P}(y_i\mid \hat\mu, \hat\sigma^2) \]

We use the logarithm, cause for large ns, the product will be a very small number cause probabilities are \(0 > p < 1\):

\[ \sum_\text{i=1}^n \text{log } \mathcal{P}(y_i\mid \hat\mu, \hat\sigma) \]

Logarithm is turning product \(\prod\) into sum \(\sum\) cause log(.5) + log(.2) is the same as log(.5 * .2).

model_0 <- lm(log_rt ~ 1, data = blomkvist)
mu_hat <- coef(model_0) # MLE of mu
sigma_hat <- sigma(model_0) # MLE of sigma

# log probability of an observed outcome for given MLEs.
logp <- dnorm(blomkvist$log_rt, mean = mu_hat, sd = sigma_hat, log = TRUE)
sum(logp) # sum of log probabilities is the log likelihood

[1] -17.6

logLik(model_0)

'log Lik.' -17.6 (df=2)

Work through exercise script: log_likelihood_1.R

Work through exercise script: residuals.R

The RSS when using the maximum likelihood estimators is

\[ \begin{aligned} \text{RSS}&=\sum_\text{i=1}^n \mid y_i - (\hat\beta_0+\hat\beta_1\cdot x_i)\mid^2,\\ &=\sum_\text{i=1}^n\mid y_i-\hat{\mu}_i\mid^2 \end{aligned} \]

RSS is a measure of how much variance in the data is not explained by the model.

# RSS: total of unexplained variance
(rss <- sum(residuals(model)^2))

[1] 9.2

calculated on the basis of TSS which is the sum of ESS and RSS

\[ \underbrace{\sum_{i=1}^n(y_i-\bar{y})^2}_\text{TSS} = \underbrace{\sum_{i=1}^n(\hat\mu_i-\bar{y})^2}_\text{ESS} + \underbrace{\sum_{i=1}^n(y_i-\hat\mu_i)^2}_\text{RSS} \]

(tss <- ess + rss)

[1] 17.8

ess / tss

[1] 0.482

summary(model)$r.sq

[1] 0.482

Work through exercise script: rsquared.R

To overcome this spurious increase in \(R^2\), the following adjustment is applied.

\[ \begin{aligned} R^2_\text{Adj} &= 1 - (1 - R^2) \cdot \underbrace{\frac{n-1}{n-K-1}}_{\text{penalty}}\\ \end{aligned} \]

• Large number of parameters \(K\) reduces \(R^2\).
• \(K\) has a small effect for large \(n\)s.

Under null hypothesis, the following two models are identical:

\[ \begin{aligned} \mathcal{M}_0:& \text{logrt}_i \sim N(\hat\mu, \sigma^2), \hat\mu_i = \beta_0\\ \mathcal{M}_1:& \text{logrt}_i \sim N(\hat\mu, \sigma^2), \hat\mu_i = \beta_0 + \beta_\text{sex} \cdot \text{sex}_{i} \end{aligned} \]

• Hence, \(R^2=0\) is only true if \(\beta_\text{sex}=0\).
• Adding \(\text{sex}\) in \(\mathcal{M}_1\) does not explain any variance in the data that model \(\mathcal{M}_0\) isn’t already accounting for.
• \(\mathcal{M}_0\) is nested in \(\mathcal{M}_1\): all the \(K_0\) predictors in \(\mathcal{M}_0\) are also presented in \(\mathcal{M}_1\).

\[ \text{DIC} = -2 \cdot \text{log} \hat{\mathcal{L}} \]

where \(\hat{\mathcal{L}}\) is the likelihood of the model using the MLEs.

Lower deviance indicates better model fit

deviance(model_2)

[1] 17.2

Equivalent to RSS for generalised linear models.

sum(residuals(model_2)^2)

[1] 17.2

… including the varying-intercepts (model_3) and the varying-intercepts and varying-slopes model (model_4): comparing_multiple_nested_models.R

# Specify models
model_1 <- lm(log_rt ~ 1, data = blomkvist)
model_2 <- lm(log_rt ~ sex, data = blomkvist)
model_3 <- lm(log_rt ~ sex + age, data = blomkvist)
model_4 <- lm(log_rt ~ sex * age, data = blomkvist)

# Compare models
anova(model_1, model_2, model_3, model_4)

Analysis of Variance Table

Model 1: log_rt ~ 1
Model 2: log_rt ~ sex
Model 3: log_rt ~ sex + age
Model 4: log_rt ~ sex * age
  Res.Df   RSS Df Sum of Sq      F  Pr(>F)
1    265 17.78                            
2    264 17.22  1      0.56  15.97 8.4e-05
3    263  9.25  1      7.97 226.81 < 2e-16
4    262  9.20  1      0.05   1.37    0.24

• Model comparisons don’t tell us which model is the best but which model is worse (i.e. “all models are wrong”).
• Failing to find evidence for an alternative model is not mean that the simpler model is the best
• It also doesn’t mean that the assumption of the more complex model is not theoretically relevant,
• or that both models are equally good.

• Comparing models with nested sets of predictors.
• Evaluation is based on likelihood of data under the model.
• Penalising more complex models by considering number of parameters.
• Other comparison techniques for models that are not nested within each other: e.g. different probability models (e.g. Poisson vs negative binomial).
• Cross-validation (e.g. “LOO”: Leave-one-out CV, k-fold CV, Monte Carlo CV).

Next week

• Thom Baguley’s session on generalized regression models for binomial and count data.

Formative assessment

• Data-analysis report using linear regression to address a psychologically relevant question.
• One reproducible data analysis hence 4 pages max (written in RMarkdown).
• Submit single zip archive via Dropbox until 24th October 2025, 2pm.