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

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

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

• Are data are compatible with the model: what is the probability of the data according to the model.
• If the probability of observing the data is higher under model \(\mathcal{M}_0\) than under \(\mathcal{M}_1\), which model is more compatible with the data? Therefore which model is better?
• The probability of the data according to the model is the model’s likelihood.

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

\[ y_i \sim 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\) (based on a particular model).

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

The joint probability is the product, if all \(y\)s are conditionally independent of one another.

\[ \prod_\text{i=1}^n 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 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 (i.e. numeric underflow) cause probabilities are \(0 > p < 1\):

\[ \sum_\text{i=1}^n \text{log } 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: exercises/log_likelihood_1.R

The sum of squared residuals 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; or the model’s lack of fit.

# 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

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} \]

\(\text{ for } i \in 1\dots n\)

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

# 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

Check exercises/comparing_multiple_nested_models.R

• Model comparisons don’t tell us which model is the best (i.e. “all models are wrong”).
• Failing to find evidence for an alternative model is not mean that the simpler model is the best, nor does it mean that the assumption of the more complex model is not theoretically relevant.

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