Why does log likelihood function for a model use SSE/n and not SSE/df? I'm trying to find out how log-likelihood function works for linear regression. I found the formula here and here. Making some experiments with it (see code below), I was quite surprised that the likelihood uses SSE/n instead of MSE (SSE/df). MSE was used everywhere up to now! I thought MSE is much better estimator of $\sigma^2$ mentioned in the formula in the 1st resource (page 6) - the actual residual variance. But the 2nd resource and my experiment clearly states that $\sigma^2$ is defined as SSE/n (where n is length of the outcome variable vector).
Here is the code to play with:
set.seed(128)
y = c(rnorm(200, 20, 4), rnorm(300, 30, 4), rnorm(400, 40, 4), rnorm(500, 50, 4))
cat1 = as.factor(c(rep(1, 200), rep(2, 300), rep(3, 400), rep(4, 500)))
rand_order = sample(1:length(cat1))
cat2 = cat1[rand_order]
cat2y = c(rep(1, 200), rep(-2, 300), rep(3, 400), rep(-4, 500))
y = y + cat2y[rand_order]
m1 = lm(y ~ 0 + cat1 + cat2)

# logLik using residual degrees of freedom (-3941.94):
-length(m1$model$y)/2*log(2*pi) - length(m1$model$y)/2*log(sum((m1$residual)^2)/m1$df.residual) - 1/2*m1$df.residual

# logLik using N (-3941.931)
-length(m1$model$y)/2*log(2*pi) - length(m1$model$y)/2*log(sum((m1$residual)^2)/length(m1$model$y)) - 1/2*length(m1$model$y)

# real logLik (-3941.931)
logLik(m1)

 A: The same situation happens in the simplest Normal means model: $Y=\mu + \epsilon$, with $\epsilon$ ~ N(0,$\sigma^2$).
The MLE of $\sigma^2$ is the sum of squares about the mean divided by n: $$S^2=\sum \frac{(Y_i-\bar{Y})^2}{n}$$.
However, this quantity is a biased estimator: $E(S^2) \neq \sigma^2$. Dividing the sum of squares by $n-1$ instead of $n$ gives an unbiased estimator of the variance.
Furthermore, $\sum \frac{(Y_i-\bar{Y})^2}{\sigma^2(n-1)}$ has a $\chi^2$ distribution with $n-1$ degrees of freedom. And ratios of independent mean squares will have the F distribution.
Statistical inference offers several criteria for "bestness" in estimators: these include unbiasedness, minimum variance, minimizer of a loss function and (increasingly) predictive accuracy. Being the maximum likelihood estimator is also considered desirable, because the likelihood (supposedly) contains all relevant information about the model (this is debatable and has been debated). The MLE usually manages to be asymptotically unbiased and efficient. What this means for a finite sample depends on the model.
Linear regression with Normal errors works much like the simple normal means model that I gave here: the unbiased estimator for the variance is not the MLE. 
The unbiased estimator is preferred because of its nice distributional properties. Geometrically, you are partitioning $R^n$ into two linear subspaces: one to contain the model and one to contain the residuals. The dimension of the residual space is the degrees of freedom.
A: To expand on a very good answer that Placidia gave:


*

*Unbiasedness is not necessarily the best possible property for an estimator. Shrinkage estimators applied in situations of multiple collinearity or with many possible regressors (lasso) are intentionally biased, and this is done to improve their other properties (easier interpretation of the results). Any Bayesian posterior mean or posterior mode estimator with an informative prior is biased; it does not mean that we want to dismiss this whole area of statistics.

*As far as other criteria of performance of a statistical estimator go, the mean squared error (MSE) is a popular criterion: by how far the estimator is off, on average, disregarding whether it is biased or not. The best estimator of the population $\sigma^2$ is then the one that has not $n-1$, not $n$, but $n+1$ in the denominator. So if your inferential target is the population variance, you might want to use the estimator that divides the sum of squared errors by $n+1$.

*The observation that the MLE $S^2$ does not appear to make much sense in regression context has been made earlier, of course, and corrections have been proposed to force it to use the "right" degrees of freedom. This is the idea of restricted maximum likelihood (REML), where the estimators are defined conditioning on the residual subspace that has the "right" dimension.

*Yet another useful property of the ML is transformation invariance. If $S^2$ is the MLE of $\sigma^2$, then automatically $S$ is the MLE of $\sigma$, and $\ln S$ is the estimator of $\ln(\sigma)$. This comes handy in software code: maximizing with respect to $S$ or $S^2$ is being complicated by the quantity being non-negative, while maximizing with respect to $\ln S$ does not involve any constraints. (You'd observe that by Jensen's inequality, unbiasedness of an estimator $$s^2 = \frac1{n-1} \sum (Y_i - \bar Y)^2$$ is easily destroyed by any transformation: $s$ is not an unbiased estimate of $\sigma$. In fact, the unbiased estimator of $\sigma$ is pretty difficult to construct, and I won't be too much ashamed to admit that I don't know one off the top of my head.)

