How to estimate logit model I am trying to understand how to fit a logit model using maximum likelihood described in a paper:
$$p_{it}=\frac{exp(\alpha+\beta q_{it})}{1+exp(\alpha+\beta q_{it})} $$
where
$$q_{it}=\frac{q_0d^{t-1}+x_{i1}d^{t-2}+...+x_{it-2}d+x_{it-1}}{d^{t-1}+d^{t-2}+...+1}$$
My question is how to estimate $\alpha$, $\beta$, $q_0$ and $d$ and also their standard errors.
I am thinking to write out the expressions:
$$q_{it}=a+b_1x_{i1}+b_2x_{i2}+...+b_{t-1}x_{it-1}$$ such that
$$a=\alpha+\frac{\beta q_0d^{t-1}}{d^{t-1}+...+1}$$
$$b_1=\frac{\beta}{d^{t-1}+...+1}d^{t-2}$$
$$...$$
$$b_{t-1}=\frac{\beta}{d^{t-1}+...+1}$$
I kind of know how to solve the equations to get $\beta$ and $d$, but I don't know how to get $\alpha$ and $q_0$. And the model is over parameterized ($t=75$), so I do not know how to get the standard errors. So far I think of it in terms of regressions, and perhaps it is not a right way to get the parameters, but I have no clue how to do it.
 A: It is unlikely that there is a closed form solution for the ML-estimators of the parameters. You must thus find the values numerically with an optimization (minimization/maximization) algorithm.
The outline is as follows

*

*Define the log-likelihood function as $$l = \sum_{i,t} \log p_{it}$$


*Use a software library routine, e.g. optim in standard R, to find the parameters that maximize l. Caveat: if the function minimizes the given function (like optim does), you must replace $l$ with $-l$.
To estimate standard errors for the parameters, there is a complicated method based on the inverse Hessian matrix (which I occasionally observe to be numerically unstable). A more robust and very simple method is instead to use the jacknife method, which I have described in section 5.2 of this article:

Dalitz: "Construction of confidence intervals." Technical Report No. 2017-01, pp. 15-28, Hochschule Niederrhein, Fachbereich Elektrotechnik und Informatik, 2017

A: If you are open to a Bayesian approach and willing to specify at least weakly informative proper prior distributions for the model parameters, then using the brms R package should make this relatively straightforward. This is basically a logistic regression - I assume this is what you meant by a logit-model, if you are modeling logit-transformed probabilities that you observe instead of Bernoulli trials, not much changes below except that you'd e.g. use beta-regression - with a non-linear function of the predictors, which is a feature that package implements easily.
Taking a Bayesian approach also helps with over-parameterized models: if you specify not-too-weak priors, these will to some extent resolve issues with model parameters being not identified or only poorly identified.
Something along the lines below should let you fit this model with a bit of experimentation. I've put in some placeholder prior distributions, but I have no idea whether those make any sense in your context (who knows, maybe some model parameters should even have joint priors).
library(brms)

brmfit1 = brm( family="bernoulli", 
     bf( y ~ alpha + beta * ( q0 * d^(t-1) + ...  ) / ( d^(t-1) + d^(t-2) + ... + 1 ),
         nl = TRUE,
         alpha ~ 1, beta ~ 1, q0 ~ 1, d ~ 1), # assumes the t and x_it are observed covariates
     prior=c(prior(normal(0,1), nlpar="alpha"),
             prior(normal(0,1), nlpar="beta"),
             prior(normal(0,1), nlpar="q0"),
             prior(lognormal(0,1), nlpar="d")),
     control=list(adapt_delta=0.99))
prior_summary(brmfit1
summary(brmfit1)
stancode(brmfit1) # if you want to see the Stan code generated by brms, that way you could modify it yourself and then use it with rstan

