-
Notifications
You must be signed in to change notification settings - Fork 0
/
102_logit.Rmd
65 lines (39 loc) · 868 Bytes
/
102_logit.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
# Logit
<!-- Discussion of the model goes here
Follow up with R implementation -->
```{r}
inv_logit <- function(x) {
p <- 1 / (1 + exp(-x))
return(p)
}
logit <- function(data, formula) {
mf <- model.frame(formula, data, na.action=na.exclude)
Y <- model.response(mf, type="double")
X <- model.matrix(formula, mf)
b <- numeric(ncol(X))
logitll <- function(b) {
p <- inv_logit((X %*% b))
lli <- Y*log(p) + (1-Y)*log(1-p)
return(-sum(lli))
}
est <- optim(b, logitll, method="BFGS")
b <- est$par
names(b) <- colnames(X)
return(b)
}
```
Testing the function:
```{r}
set.seed(42)
N <- 1000
X1 <- rnorm(N)
X2 <- rnorm(N)
ystar <- 1.5 + 2*X1 + 3*X2 + rnorm(N)
y <- rbinom(N, 1, inv_logit(ystar))
df <- data.frame(
y, X1, X2
)
f <- y ~ X1 + X2
coef(glm(f,binomial,df))
logit(df, f)
```