7
$\begingroup$

It is said that Least Squares estimates would differ from Maximum Likelihood estimates if the underlying data were non-normal. This should be the reason why LS estimates can be used in linear regression, while for generalized linear models, the parameters are derived from ML estimates. Even though I read this many times, I never saw this in action. Thus, I wonder if one could provide an example in R where the estimates using both methods clearly deviate from each other.

$\endgroup$
1
  • $\begingroup$ Just to clear things up: a least square fit IS a maximum likelihood fit! It's one where you assume for every point a Gaussian model that has the prediction as mean, 1 as std and the data as "x". You can try that: take the negative log of a Gaussian and the least square pops out (which is equivalent when minimized to minimizing the negative log-likelihood). Therefore, the least square works if a gaussian with mean "predictions" is the actual model the data was drawn from. Your just maximum likelihood-fit many gaussians in the least square. No magic. $\endgroup$
    – Mayou36
    Apr 20 at 17:07

3 Answers 3

4
$\begingroup$

Let's simulate data where the outcome y is gamma-distributed, with its scale parameter modeled as an exponential function of x:

> set.seed(42)
> n <- 5e3
> x <- rnorm(n)
> y <- rgamma(n, shape=2, scale=exp(.9 + .7*x))
> 
> ## Shapiro test
> shapiro.test(y)

    Shapiro-Wilk normality test

data:  y
W = 0.6585, p-value < 2.2e-16

The Shapiro-Wilk test confirms non-normality. Histogram and QQ-plot look - surprise - gamma-distributed.

enter image description here

Maximum likelihood

Using maximum likelihood, we derive the expected (modeled) estimates for the intercept and slope.

> ## ML
> L <- function(theta, ...) {
+   shape <- 2  ## a priori shape parameter for demonstration
+   scale <- exp(theta[1] + theta[2]*x)
+   loglik <- sum(dgamma(y, shape, scale=scale, log=TRUE))
+   return(-loglik)
+ }
> optim(c(0, 0), L, x=x, y=y)$par
[1] 0.9019497 0.7098924

OLS

However, OLS does not properly model the data due to assumption of normally distributed residuals:

> OLS <- lm(y ~ x)
> as.numeric(OLS$coe)
[1] 6.337453 4.388758

GLM

Using glm with a Gamma family and log link accurately models the data with maximum likelihood:

> GLM <- glm(y ~ x, family=Gamma(link="log"))
> as.numeric(GLM$coe)
[1] 1.5949384 0.7098734

After subtracting log(2) from the intercept to compensate for the fixed shape parameter, the GLM results are consistent with those obtained from the manual approach.

> as.numeric(GLM$coe) - c(log(2), 0)
[1] 0.9017912 0.7098734

Conversely, we could adapt the function used in optim to match the GLM results.

> L_glm <- function(theta, y, x) {
+   shape <- 2
+   log_mean <- theta[1] + theta[2]*x
+   scale <- exp(log_mean)/shape  ## scale as mean by (fixed) shape parameter
+   loglik <- -sum(dgamma(y, shape, scale=scale, log=TRUE))
+   return(loglik)
+ }
> optim(c(0, 0), L_glm, x=x, y=y)$par
[1] 1.5951432 0.7098887

Residuals

Examining the residuals shows, while the residuals of the GLM (blue) adequately cluster around zero, those of the OLS model (red) clearly violate the OLS-assumption of normal distributed residuals.

> quantile(GLM$resid) |> zapsmall()
       0%       25%       50%       75%      100% 
-0.978066 -0.517014 -0.159916  0.338285  4.952284 
> quantile(OLS$resid) |> zapsmall()
       0%       25%       50%       75%      100% 
-15.11840  -3.37093  -0.70010   1.97892  87.12154 

enter image description here

$\endgroup$
3
$\begingroup$

The ordinary least-squares (OLS) estimator is the maximum likelihood estimator (MLE) assuming Gaussian errors. That is, the OLS estimator is the generalized linear model (GLM) for the Gaussian distribution (and linear link function).

Below is an example comparing OLS with the GLM with binomial distribution and logit link function (that is, logistic regression):

> data = iris
> data$outcome <- data$Species == "virginica"
> lm(outcome ~ Sepal.Length, data=data)

Call:
lm(formula = outcome ~ Sepal.Length, data = data)

Coefficients:
 (Intercept)  Sepal.Length  
     -1.7962        0.3644  

> glm(outcome ~ Sepal.Length, data=data, family="gaussian")

Call:  glm(formula = outcome ~ Sepal.Length, family = "gaussian", data = data)

Coefficients:
 (Intercept)  Sepal.Length  
     -1.7962        0.3644  

Degrees of Freedom: 149 Total (i.e. Null);  148 Residual
Null Deviance:      33.33 
Residual Deviance: 19.76        AIC: 127.7
> glm(outcome ~ Sepal.Length, data=data, family="binomial")

Call:  glm(formula = outcome ~ Sepal.Length, family = "binomial", data = data)

Coefficients:
 (Intercept)  Sepal.Length  
     -16.320         2.592  

Degrees of Freedom: 149 Total (i.e. Null);  148 Residual
Null Deviance:      191 
Residual Deviance: 117.3        AIC: 121.3
$\endgroup$
2
$\begingroup$

Here is a non-regression example distinguishing between maximum likelihood and least squares, which might also illustrate the bias-variance trade-off.

Suppose you have a sample sized $n$ from $\mathcal N(\mu,\sigma^2)$ where $\mu$ and $\sigma^2$ are unknown. It is natural to estimate $\mu$ using $\bar x$ and this is indeed the least squares and maximum likelihood estimate. For $\sigma^2$ the question is harder, though natural to use some fraction of $\sum (x_i-\bar x)^2$. We can find that $\frac1{\sigma^2}\sum (X_i-\bar X)^2 \sim \chi^2_{n-1}$.

Let's set $S^2_k = \frac1k\sum (X_i-\bar X)^2$ and consider different values of $k$.

  • The estimator $S^2_{n-1}$ has an expectation of $\sigma^2$ since $\chi^2_{n-1}$ has expectation $n-1$, and so is unbiased (and is the best unbiased estimator in various senses). This is associated with Bessel's correction when estimating variances.

  • The estimator $S^2_{n}$ is the maximum likelihood estimator of $\sigma^2$ though, being smaller, its expectation is slightly lower and so is biased; its variance is also smaller and, after combining these, its expected squared error is lower than the unbiased estimator.

  • The estimator $S^2_{n+1}$ minimises the expected squared error here. Being even smaller, it is more biased but has a smaller variance, and is (when sampling from a normal distribution) the point where the increase in the square of the bias matches the decrease in the variance.

All this can be shown analytically, but to illustrate the points with a simulation in R, this takes a million cases of a sample sized $n=10$, looking at $s^2_{n-1}, s^2_n, s^2_{n+1}$, and showing the average log-likelihood, the average bias, the average variance and the average sum-of-squares error:

set.seed(2024)
n <- 10
mu <- 1234
sigmasq <- 56
cases <- 10^6
xmat <- matrix(rnorm(n*cases, mu, sqrt(sigmasq)), ncol=n)
xbar <- rowMeans(xmat)
sumsquarexminusxbar <- rowSums((xmat-xbar)^2)

results <- matrix(nrow=4, ncol=5)
colnames(results) <- c("S^2_n-2", "S^2_n-1", "S^2_n", "S^2_n+1", "S^2_n+2")
rownames(results) <- c("meanloglike", "bias", "var", "sqerror")

for (k in (n-2):(n+2)){
  ssq_k <- sumsquarexminusxbar / k
  results[, k-n+3] <- c(
     mean(rowSums(log(dnorm(xmat, xbar, sqrt(ssq_k))))), 
     mean(sigmasq - ssq_k),
     mean((ssq_k - mean(ssq_k))^2),
     mean((sigmasq - ssq_k)^2))
  }

giving the simulation results

> print(results)
               S^2_n-2      S^2_n-1      S^2_n   S^2_n+1   S^2_n+2
meanloglike -33.330535 -33.24162001 -33.214817 -33.23827 -33.30321
bias         -7.013697  -0.01217492   5.589043  10.17186  13.99087
var         882.764300 697.49278015 564.969152 466.91665 392.33969
sqerror     931.956242 697.49292837 596.206549 570.38333 588.08410

You can see that the middle $S^2_n$ column has the highest (least negative) log-likelihood, that the $S^2_{n-1}$ column has zero bias (allowing for a little simulation noise), and that the $S^2_{n+1}$ column has the lowest expected squared error (essentially the variance plus the square of the bias), illustrating the assertions above.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.