Logistic Regression
The following demo regards a standard logistic regression model via maximum likelihood or exponential loss. This can serve as an entry point for those starting out to 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. Exponential loss is not confined to the standard GLM setting, but is widely used in more predictive/‘algorithmic’ approaches e.g. in machine learning and elsewhere.
This follows the linear regression model approach.
Data Setup
Predictors and target. This follows the same approach as the linear regression example, but now draws the target variable from the binomial distribution with size = 1
.
library(tidyverse)
set.seed(1235) # ensures replication
= 10000 # sample size
N = 2 # number of desired predictors
k = matrix(rnorm(N * k), ncol = k)
X
# the linear predictor
= -.5 + .2*X[, 1] + .5*X[, 2] # increasing N will get estimated values closer to these
lp
= rbinom(N, size = 1, prob = plogis(lp))
y
= data.frame(X, y) dfXy
Functions
A maximum likelihood approach.
<- function(par, X, y) {
logreg_ml # Arguments
# par: parameters to be estimated
# X: predictor matrix with intercept column
# y: target
# setup
= par # coefficients
beta = nrow(X)
N
# linear predictor
= X %*% beta # linear predictor
LP = plogis(LP) # logit link
mu
# calculate likelihood
= dbinom(y, size = 1, prob = mu, log = TRUE) # log likelihood
L # L = y*log(mu) + (1 - y)*log(1-mu) # alternate log likelihood form
-sum(L) # optim by default is minimization, and we want to maximize the likelihood
# (see also fnscale in optim.control)
}
Another approach via exponential loss function.
<- function(par, X, y) {
logreg_exp # Arguments
# par: parameters to be estimated
# X: predictor matrix with intercept column
# y: target
# setup
= par # coefficients
beta
# linear predictor
= X %*% beta # linear predictor
LP
# calculate exponential loss function (convert y to -1:1 from 0:1)
= sum(exp(-ifelse(y, 1, -1) * .5 * LP))
L }
Estimation
Setup for use with optim.
= cbind(1, X)
X
# initial values
= rep(0, ncol(X))
init names(init) = c('intercept', 'b1', 'b2')
= optim(
fit_ml par = init,
fn = logreg_ml,
X = X,
y = y,
control = list(reltol = 1e-8)
)
= optim(
fit_exp par = init,
fn = logreg_exp,
X = X,
y = y,
control = list(reltol = 1e-15)
)
= fit_ml$par
pars_ml = fit_exp$par pars_exp
Comparison
Compare to glm
.
= glm(y ~ ., dfXy, family = binomial) fit_glm
intercept | b1 | b2 | |
---|---|---|---|
fit_ml | -0.504 | 0.248 | 0.503 |
fit_exp | -0.503 | 0.245 | 0.501 |
fit_glm | -0.504 | 0.248 | 0.503 |
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_logistic.R