-
Notifications
You must be signed in to change notification settings - Fork 17
/
logistic-regression.Rmd
146 lines (99 loc) · 3.33 KB
/
logistic-regression.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
# 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][Linear Regression] 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`.
```{r data-setup-logreg}
library(tidyverse)
set.seed(1235) # ensures replication
N = 10000 # sample size
k = 2 # number of desired predictors
X = matrix(rnorm(N * k), ncol = k)
# the linear predictor
lp = -.5 + .2*X[, 1] + .5*X[, 2] # increasing N will get estimated values closer to these
y = rbinom(N, size = 1, prob = plogis(lp))
dfXy = data.frame(X, y)
```
## Functions
A maximum likelihood approach.
```{r logreg_ml}
logreg_ml <- function(par, X, y) {
# Arguments
# par: parameters to be estimated
# X: predictor matrix with intercept column
# y: target
# setup
beta = par # coefficients
N = nrow(X)
# linear predictor
LP = X %*% beta # linear predictor
mu = plogis(LP) # logit link
# calculate likelihood
L = dbinom(y, size = 1, prob = mu, log = TRUE) # log likelihood
# 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.
```{r logreg_exp}
logreg_exp <- 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
# calculate exponential loss function (convert y to -1:1 from 0:1)
L = sum(exp(-ifelse(y, 1, -1) * .5 * LP))
}
```
## Estimation
Setup for use with <span class="func" style = "">optim</span>.
```{r logreg-est}
X = cbind(1, X)
# initial values
init = rep(0, ncol(X))
names(init) = c('intercept', 'b1', 'b2')
fit_ml = optim(
par = init,
fn = logreg_ml,
X = X,
y = y,
control = list(reltol = 1e-8)
)
fit_exp = optim(
par = init,
fn = logreg_exp,
X = X,
y = y,
control = list(reltol = 1e-15)
)
pars_ml = fit_ml$par
pars_exp = fit_exp$par
```
## Comparison
Compare to `glm`.
```{r logreg-compare}
fit_glm = glm(y ~ ., dfXy, family = binomial)
```
```{r logreg-compare-show, echo=FALSE}
rbind(
fit_ml = pars_ml,
fit_exp = pars_exp,
fit_glm = coef(fit_glm)
) %>%
kable_df()
```
## Python
The above is available as a Python demo in the [supplemental section](#python-logreg).
## Source
Original code available at https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/standard_logistic.R