Computing Log-likelihood Model Manually for Logit Model

Question

A log-likelihood function takes on the form of: logL$=\Sigma^{n}_{i=1}(y_ilog(\frac{e^{\beta' x_i}}{1+e^{\beta'x_i}} )+(1-y_i)(\frac{e^{\beta' x_i}}{1+e^{\beta'x_i}})$

My logit model is estimated as follows:

glm.logit=glm(model,binomial(),data)

Estimating $y_i$

yi=data$y

x's as a matrix of the dependent variable

 xi=cbind(data$x1+data$x2+...)

Taking the estimates of $\beta$ from my model

betai=coef(glm.logit)

Putting the these together:

xibetai<-xi%*%t(betai)

Estimating the logistic form:

logiti<-exp(xibetai)/(1+exp(xibetai))

Putting everything together in the form of a log-likelihood model:

LogLi<-yi%*%log(logiti)+(1-yi)%*%log(1-logiti)

Issue with the fact that you can't take take the log of a negative value, so I only get about half of my coefficients in the output.

Answer 1

Both your log-likelihood equation and your R code have some errors:

$$ \begin{split} \eta_i & = \mathbf{X}_i \boldsymbol{\beta} \quad \textrm{(so far so good)} \\ p_i & = \exp(\eta_i)/(1+\exp(\eta_i)) = 1/(1+\exp(-\eta_i)) \\ L_i & = y_i \log p_i + (1-y) \log(1-p_i) \\ L & = \sum L_i \end{split} $$

In R you can use plogis() for the logistic ($1/(1+\exp(-\eta_i))$):

p <- plogis(X %*% beta)
sum(y*log(p) + (1-y)*log(1-p))

(or dbinom(y, prob=p, size=1, log=TRUE))

If your coefficients include an intercept, $\mathbf{X}$ should include a column of 1s at the corresponding position (usually the first column). You can check various intermediate results with predict(fitted_model, type="link") ($\eta$), predict(fitted_model, type="response") ($p$).