# Maximize *custom* likelihood function for logistic regression

My dataset contains n observations $X_i$ of n individuals and I want to predict a binary outcome $Y_i$.

1. ### logistic regression model

It is fair to assume that this $Y_i$ is the realization of a Bernoulli experiment and is 1 with probability $\pi_i$. Under the logistic regression model, I model $\log\big(\frac{\pi_i}{1-\pi_i}\big) = \mathbf{X_i} \mathbf{\beta}$ which is fine by me. $\beta$ is solved for by optimizing the log-likelihood of the data under the model. The individual log-likelihood is

$$\log L_\beta(i)=\pi_i^{y_i}(1-\pi_i)^{1-y_i}$$

The total log likelihood is the sum of the individual likelihood, expressed in terms of $\beta$.

2. ### Changes I want to make

The trick is, I also have historical data about previous outcome ($Y$) for some individuals! I do not have the associated previous $X$ though. In my case, it is fair to assume that the probability $\pi_i$ does not change in the life of an individual so I want to take into account the previous outcome in the likelihood!

If I call $n(i)$ the total number of observations ($\geq1$) for individual $i$ and $Y_{i,j}$ the $j^{th}$ outcome observed for individual $i$, my new individual likelihood is:

$$\log L_\beta(i)= \sum_j^{n(i)} \pi_i^{y_{i,j}}(1-\pi_i)^{1-y_{i,j}}$$

And the total likelihood to optimize is the sum of this over $i$.

3. ### My questions

• Does it make sense to do this or is there any obvious flaw I missed? I am open to alternatives.

• Can I implement this using glm (R)? I would need to define my own family, I don't know if it is possible. I guess the likelihood above would have to be in the exponential family and I am not sure it is since the number of observations for a given individual is not fixed.

• Do I need to optimize manually (using optim for R for ex) to fit my model? I don't have a problem with this approach per se, but in the end I would like to use this likelihood in a LASSO or elastic net framework and it would be so much more convenient and effective to be able to use glmnet!

• How do the predictors in your data (your $X$'s) vary over the life of an individual? – Matthew Drury Sep 28 '15 at 16:18
• Some of them have relatively low variance, like Height, Weight, genetic information, etc. Others are specific to the current events leading to the outcome $Y_i$ and may have been different during the previous experiments. But anyway, the fact that $\pi_i$ is the same for all $Y_{i,j}$ is really fair for this data. – asac - Reinstate Monica Sep 28 '15 at 17:26
• I wonder if a Bayesian approach would be appropriate here. You could give each person a prior. Everyone w/ no prior observations would be set at the mean, whereas those w/ prior observations would be shifted from the mean according to their observed previous outcomes. – gung - Reinstate Monica Sep 29 '15 at 1:38

I realized something:

\begin{align} \log L_\beta(i) & = \sum_j^{n(i)} \pi_i^{y_{i,j}}(1-\pi_i)^{1-y_{i,j}} \\ & = \pi_i^{y_{i,1}}(1-\pi_i)^{1-y_{i,1}} + \sum_j^{n(i)-1} \pi_i^{y_{i,j+1}}(1-\pi_i)^{1-y_{i,j+1}} \\ & = \pi_i^{y_{i}}(1-\pi_i)^{1-y_{i}} + \pi_i^{y_{i,2}}(1-\pi_i)^{1-y_{i,2}} & \mbox{assuming } n(i)=2 \mbox{ to simplify} \end{align}

Since $\pi_i = f(X_i)$, the individual likelihood for the $i^{th}$ observation is equal to the likelihood of $n(i)$ observations, all with same $X=X_i$ but with outcomes resp. $Y_{i,1}, ..., Y_{i,n(i)}$ .

I implemented the likelihood described in my question and optimized it manually with optim. I also tried the approach described in my answer which consists of transforming the dataset to create "fake" observations corresponding to previous history, and using glm on the transformed data. The coefficients estimates match!

This great because it allows me to use glmnet in a plug-and-play manner.

Here is what I would do: Notice that I am taking the mean value for the covariates at the individual level and the log likelihood will be: I think you need to hard code this in R.

• It isn't apparent how the reader is intended to interpret these formulas as analyses or operations that you would "do." With your last parenthetical remark even you seem unsure! I think much more explanation may be needed if you want people to understand what you are recommending. – whuber Sep 28 '15 at 19:50
• Indeed I am not sure how to interpret your first two formulas. I do not have any $x_{i,j}$. I have only one vector $x_i$ per individual, corresponding to the outcome $Y_{i,1}$. I also happen to have some $Y_{i,j}$ with $j>1$ from previous outcome for some individuals, but I don't have any more data about these. I edited my text a bit in case that was not clear. – asac - Reinstate Monica Sep 28 '15 at 19:58
• @ antonie-sac: The 'Pi' is the probability of event, which is what you have in your first bullet point in the original question. The problem you have is you observe each individual multiple times and this is Logit in a panel setup. However you are assuming the probability of event (Pi) only changes at the individual level. So my first equation shows the probability of event (Pi) that does not change with time. Since you have only one covariate, you take the average of that covariate across time at the individual. This will provide a time invariant ‘Pi’. – subra Sep 28 '15 at 20:12
• @subra: I don't get your equation: $\bar{x_{i}} = \sum_j x_{i,j}/j$ because I don't have several measurements of $x$ for a given individual. I observe several outcomes $Y$ only. – asac - Reinstate Monica Sep 28 '15 at 20:22
• @ Antonio: The value of predictor variable does not change with time for individuals? Is that rite? In that case you can just use 'x' instead of the average. That should give a unique 'Pi' (probability of event) at the individual level – subra Sep 28 '15 at 20:44