# Linear Regression

We start our demonstrations with a standard regression model via maximum likelihood or least squares loss. Also included are examples for QR decomposition and normal equations. This can serve as an entry point for those starting out in the wider world of computational statistics, as maximum likelihood is the fundamental approach used in most applied statistics, but which is also a key aspect of the Bayesian approach. Least squares loss is not confined to the standard regression setting, but is widely used in more predictive/‘algorithmic’ approaches e.g. in machine learning and elsewhere. You can find a Python version of the code in the supplemental section.

## Data Setup

We will simulate some data to make our results known and easier to manipulate.

library(tidyverse)

set.seed(123)  # ensures replication

# predictors and target
N = 100 # sample size
k = 2   # number of desired predictors
X = matrix(rnorm(N * k), ncol = k)
y = -.5 + .2*X[, 1] + .1*X[, 2] + rnorm(N, sd = .5)  # increasing N will get estimated values closer to these

dfXy = data.frame(X, y)

## Functions

A maximum likelihood approach.

lm_ml <- function(par, X, y) {
# par: parameters to be estimated
# X: predictor matrix with intercept column
# y: target

# setup
beta   = par[-1]                             # coefficients
sigma2 = par[1]                              # error variance
sigma  = sqrt(sigma2)
N = nrow(X)

# linear predictor
LP = X %*% beta                              # linear predictor
mu = LP                                      # identity link in the glm sense

# calculate likelihood
L = dnorm(y, mean = mu, sd = sigma, log = TRUE) # log likelihood
# L =  -.5*N*log(sigma2) - .5*(1/sigma2)*crossprod(y-mu)    # alternate log likelihood form

-sum(L)                                      # optim by default is minimization, and we want to maximize the likelihood
}

An approach via least squares loss function.

lm_ls <- function(par, X, y) {
# arguments-
# par: parameters to be estimated
# X: predictor matrix with intercept column
# y: target

# setup
beta = par                                   # coefficients

# linear predictor
LP = X %*% beta                              # linear predictor
mu = LP                                      # identity link

# calculate least squares loss function
L = crossprod(y - mu)
}

## Estimation

Setup for use with optim.

X = cbind(1, X)

Initial values. Note we’d normally want to handle the sigma differently as it’s bounded by zero, but we’ll ignore for demonstration. Also sigma2 is not required for the LS approach as it is the objective function.

init = c(1, rep(0, ncol(X)))
names(init) = c('sigma2', 'intercept', 'b1', 'b2')

fit_ML = optim(
par = init,
fn  = lm_ml,
X   = X,
y   = y,
control = list(reltol = 1e-8)
)

fit_LS = optim(
par = init[-1],
fn  = lm_ls,
X   = X,
y   = y,
control = list(reltol = 1e-8)
)

pars_ML = fit_ML$par pars_LS = c(sigma2 = fit_LS$value / (N-k-1), fit_LS$par) # calculate sigma2 and add ## Comparison Compare to lm which uses QR decomposition. fit_lm = lm(y ~ ., dfXy) Example of QR. Not shown. # QRX = qr(X) # Q = qr.Q(QRX) # R = qr.R(QRX) # Bhat = solve(R) %*% crossprod(Q, y) # alternate: qr.coef(QRX, y) sigma2 intercept b1 b2 fit_ML 0.219 -0.432 0.133 0.112 fit_LS 0.226 -0.432 0.133 0.112 fit_lm 0.226 -0.432 0.133 0.112 The slight difference in sigma is roughly dividing by N vs. N-k-1 in the traditional least squares approach. It diminishes with increasing N as both tend toward whatever sd^2 you specify when creating the y target above. Compare to glm, which by default assumes gaussian family with identity link and uses lm.fit. fit_glm = glm(y ~ ., data = dfXy) summary(fit_glm)  Call: glm(formula = y ~ ., data = dfXy) Deviance Residuals: Min 1Q Median 3Q Max -0.93651 -0.33037 -0.06222 0.31068 1.03991 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.43247 0.04807 -8.997 1.97e-14 *** X1 0.13341 0.05243 2.544 0.0125 * X2 0.11191 0.04950 2.261 0.0260 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for gaussian family taken to be 0.2262419) Null deviance: 24.444 on 99 degrees of freedom Residual deviance: 21.945 on 97 degrees of freedom AIC: 140.13 Number of Fisher Scoring iterations: 2 Via normal equations. coefs = solve(t(X) %*% X) %*% t(X) %*% y # coefficients Compare. sqrt(crossprod(y - X %*% coefs) / (N - k - 1)) summary(fit_lm)$sigma
sqrt(fit_glm$deviance / fit_glm$df.residual)
c(sqrt(pars_ML[1]), sqrt(pars_LS[1]))

# rerun by adding 3-4 zeros to the N
          [,1]
[1,] 0.4756489
[1] 0.4756489
[1] 0.4756489
sigma2    sigma2
0.4684616 0.4756490 

## Python

The above is available as a Python demo in the supplemental section.

## Source

Original code available at https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/standard_lm.R