Hypothesis testing: Why is a null model that fits the data well better than one that doesn't? Let's say we have two models: a null model, $M_0$, and an alternative model $M_1$.    The only difference between them is that, in $M_0$ one parameter is fixed at $0$ and in $M_1$, that parameter is fixed at the value that maximizes the likelihood of model $M_1$.  This is a typical setup for a likelihood ratio test.
My intuition is that the better $M_0$ describes the data-generating process, and thus the less the residual variation in the fitted model, the better.  By "better", I mean for a given sample size, effect size, and false positive rate, I will have more power to reject the null.
That's a bit hand-wavey.  I'll make a simulation with a linear regression model.
set.seed(27599)

lm_lrs_no_covar <- function(y, x) {
    2*(logLik(lm(formula = y ~ x)) - logLik(lm(formula = y ~ 1)))
}

lm_lrs_yes_covar <- function(y, x, z) {
  2*(logLik(lm(formula = y ~ x + z)) - logLik(lm(formula = y ~ z)))
}

n <- 1e2
num_sims <- 1e4

no_covar <- yes_covar <- rep(NA, num_sims)

for (sim_idx in 1:num_sims) {

  x <- runif(n = n)
  z <- runif(n = n)
  y <- rnorm(n = n, mean = 0.2*x + z, sd = 0.2)

  yes_covar[sim_idx] <- lm_lrs_yes_covar(y = y, x = x, z = z)
  no_covar[sim_idx] <- lm_lrs_no_covar(y = y, x = x)
}


plot(x = sort(no_covar),
     y = sort(yes_covar),
     type = 'l')
abline(a = 0, b = 1)

This plot shows that the LR statistic from the model with the covariate is pointwise greater than the LR statistic from the model without the covariate.

But, why, from a likelihood perspective is this so?
 A: The models are dominated by the large residual variability, implicitly set to the default value sd=1 of the rnorm() function. Differences in the covariate structure contributes little to the quality of the model, and the small difference cannot be detected with your simulations. Use sd=0.1 and you get different results that match your intuition.
    set.seed(27599)

lm_lrs_no_covar <- function(y, x) {
    2*(logLik(lm(formula = y ~ x)) - logLik(lm(formula = y ~ 1)))
}

lm_lrs_yes_covar <- function(y, x, z) {
  2*(logLik(lm(formula = y ~ x + z)) - logLik(lm(formula = y ~ z)))
}

n <- 1e2
num_sims <- 1e2

no_covar <- yes_covar <- rep(NA, num_sims)

for (sim_idx in 1:num_sims) {

  x <- runif(n = n)
  z <- runif(n = n)
  y <- rnorm(n = n, mean = 0.2*x + z,sd=0.1)

  yes_covar[sim_idx] <- lm_lrs_yes_covar(y = y, x = x, z = z)
  no_covar[sim_idx] <- lm_lrs_no_covar(y = y, x = x)
}


plot(x = sort(no_covar),
     y = sort(yes_covar),
     type = 'l')
abline(a = 0, b = 1)

A: I haven't been able to thoroughly crack this one, and it's only gotten one partial answer on SE, but here's what progress I've made.

It's not too surprising, but both the null and the alt models fit better when the covariate is included (blue and purple are above red and green).  But, as reflected in the LR, there's more improvement from null to alt when the covar is modeled (blue->purple has more improvement than red->green).
I guess my read on this is that the more correct the null model is, the more improvement we get from adding in the tested parameter.  In the absence of the covariate, the MLE of the parameter that we add to go from null to alt may be some compromise between explaining the unmodeled effect of the covariate and explaining what that parameter really does in the data-generating process.
