pweight with melogit

[Marie Kuehnzwei]

I have a Panel dataset with 260,647 data points consisting of 260,647 Respondents within 43,400 households and between 1 and 11 survey waves. Using this dataset, I am trying to run a logistic random effects regression. So far, I am using xtlogit, which works fine except that it doesn't allow me to use the survey weights that are needed to make statistical inference. Additionally, I can't account for the household level in the model, which is why I am currently using robust standard errors. I would prefer to specify a model which allows me to include my survey weight with (pweight) and to explicitly model both the panel structure and the multilevel structure. The model looks like this:

webuse union.dta, clear
xtset idcode year
xtlogit union, vce(robust) re

(I am using the union data here because I am not allowed to share the actual data that I am using)

The multilevel structure can be added using xtmelogit, but as I understand it, it is no longer part of the official stata which makes me a bit suspicious and also it doesn't allow me to use pweights (and it takes a really long time to run even for the empty model).

A good option seems to be melogit, but for some reason it doesn't converge even when I am using it to fit the empty model without the houdehold level as soon as I include the weights. The model looks like this (of course, there is no weighting variable in the union dataset. To illustrate the problem, I have created a weighting variable pw = 1, but with this, the model runs. So this is just an illustration of what my code looks like and unfortunately can't replicate the problem):

gen pw = 1
melogit union [pw=pw] || idcode:

With my data, the output from that model looks like this:

Fitting fixed-effects model:

Iteration 0: log likelihood = -4.512e+08
Iteration 1: log likelihood = -4.500e+08
Iteration 2: log likelihood = -4.500e+08
Iteration 3: log likelihood = -4.500e+08

Refining starting values:

Grid node 0: log likelihood = .
Grid node 1: log likelihood = .
Grid node 2: log likelihood = .
Grid node 3: log likelihood = .
(note: Grid search failed to find values that will yield a log likelihood value.)

Fitting full model:

initial values not feasible
r(1400);

Does anyone have an Idea what could be the issue here even without looking at the data? Or maybe there is another way to specify xtlogit-like models in Stata that allows both pweights and another level?

[Andrew Musau]

As you allude to, the issue is specific to your data and not with melogit and weights (as you are able to successfully run a regression with weights). If you get results without weights, you can attempt to use these as starting values for the regression with weights. Also, see if you need to rescale some variables in the regression.

melogit yvar xvars || id:
mat b= e(b)
melogit yvar xvars [pw= pw], from(b) || id:

[Marie Kuehnzwei]

Thank you very much for the idea and the code example! That made it very easy for me to implement. Unfortunately, using the starting values from the unweighted model yields the same results.

In the meantime I have also tried using svyset and svy:, which doesn't work because it requires a weight at every specified stage. Is it possible that this is also the issue with the model described above and that melogit expects a weighting variable that is constant within persons?

Or is there anything else I can try?

[Andrew Musau]

I would give gllamm from SSC a shot at estimating the model. It implements an adaptive quadrature maximization algorithm in addition to Newton-Raphson and may be able to find a starting value. Here is the syntax, where you put your probability weight variable in the highlighted.

gllamm depvar [varlist] [if exp] [in range] , i(varlist)
[ noconstant offset(varname) nrf(#,...,#)
eqs(eqnames) frload(#,...,#)^
ip(string) nip(#,...,#) peqs(eqname)
bmatrix(matrix) geqs(eqnames) nocorrel
constraints(clist) weight(varname) pweight(varname)
family(family) fv(varname) denom(varname)
s(eqname) link(link) lv(varname)
expanded(varname varname string) basecategory(#)
composite(varnames) thresh(eqnames) ethresh(eqnames)
from(matrix) copy skip long
lf0(# #) gateaux(# # #) search(#)
noest eval init iterate(#) adoonly adapt
robust cluster(varname)
level(#) eform allc trace nolog nodisplay dots
]

The option for adaptive quadrature is -adapt-. Otherwise, start with a simple model and add one regressor at a time until the point where the estimation fails and look at the issue with that added regressor. It may be that you need to rescale it - divide by some positive factor if the values are too large or multiply if the values are too small.