*! version 1.1.1 April 2026

Bivariate mixed Poisson regression (using Maximum Simulated Likelihood)

bimpoisson -- Bivariate mixed Poisson regression (MSL)

Syntax

bimpoisson
depvar1 depvar2
[indepvars]
[if] [in]
[weight],
[options]

can be used when indepvars is a common set of predictors
for each dependent variable. For this case, and also situations where the sets
of predictors for each dependent variable differ, the following syntax can be
used:

bimpoisson
equation1 equation2
[if] [in]
[weight],
[options]

where equation1 and equation2 are each specified as

( [eqname:] depvar [=] [indepvars] )

bimpoisson shares the features of all estimation commands; see help
est.

bimpoisson typed without arguments redisplays the last estimates.

Main
nsim(#) number of simulations (default value is 10)

antithetic use antithetic acceleration

halton use Halton draws (default is standard uniform pseudorandom draws)

burn(#) drop # initial elements when creating the Halton sequences (default is 5)

bias use first-order bias correction

seed(#) set the seed (relevant if using standard uniform pseudorandom draws)

rho_sf(#) set the sampling function correlation (default value is 0)

alt use alternative MSL estimation method (default uses sampling function)

rng0 set the rng state back to its initial value

Model
noconstant suppress constant term

offset1(varname) offset variable for first equation

offset2(varname) offset variable for second equation

constraints(constraints) apply specified linear constraints

SE/Robust
vce(vcetype) vcetype may be oim,
robust, cluster clustvar, opg,
bootstrap, or jackknife

Reporting
irr report incidence rate ratios

eform report exponentiated coefficients

level(#) set confidence level; default is level(95)

nocnsreport do not display constraints

display_options control
Maximization
maximize_options control the maximization process; seldom used

collinear keep collinear variables

depvar1, depvar2, indepvars, and depvar may
contain time-series operators; see tsvarlist.

by and collect are allowed; see prefix.

pweights, fweights, aweights, and iweights
are allowed; see weight.

See bimpoisson postestimation for features available after
estimation.

Description

bimpoisson fits bivariate mixed Poisson regression models by the method
of maximum simulated likelihood (MSL), following the approaches proposed by
Munkin and Trivedi (1999). To quote Munkin and Trivedi (`MT'), p. 35,
"[w]e wish to to specify a general bivariate count model in which
overdispersion comes through bivariate unobserved heterogeneity which affects
the conditional mean functions. Assume that the two count variables
y1 and y2 have a Poisson distribution conditional on
the σ-fields generated by the random variables of unobserved heterogeneity
ε1, ε2 and the set of independent explanatory variables x1
and x2:

yj | xj, εj ~ P(μj), μj > O, j = l, 2,

where we parameterize the mean parameters μ1 and μ2 as
exp(x1'ß1 + ε1) and exp(x2'ß2 + ε2),
respectively. The unobservable variables ε1 and ε2 are
assumed jointly normally distributed, i.e.

(ε1, ε2) ~ N{(0, 0), (σ1^2, ρσ1σ2, σ2^2)}, 0 <= |ρ| <= 1."

MT recommend using a robust variance estimator.

For an overview of alternative approaches to modelling bivariate count
data, and software, see Xu and Hardin (2016).

The default MSL estimator implemented in bimpoisson is that described
by MT, p. 36, in which an importance sampling function is used to address
the issue that maximization of the MSL function is problematic for values
of |ρ| close to one. (Symptoms may include continuous "(backed up)"
messages in the iteration log or "r(430) convergence not achieved".) MT
use a bivariate standard normal distribution with ρ = 0 as their
sampling function. Like Jumamyradov and Munkin (2022), we use a sampling
distribution that is bivariate normal with means zero, but with
standard deviations unrestricted. With this approach, drawn values provide more
effective coverage and improve MSL's estimation of the double integral
that specifies the bivariate density of the bivariate mixed Poisson model.
bimpoisson builds on this idea and also allows users to optionally set
a non-zero value for the sampling function correlation using the
rho_sf(#) option.

MT's MSL estimator not taking account of the maximization feasibility issue,
described in their eq. (11), p. 35, is implemented by selecting
bimpoisson's alt option.

By default and following MT, bimpoisson's two MSL estimators use
pseudorandom numbers (draws from the standard uniform distribution) when
simulating probabilities. Optionally, Halton draws can be used instead.
(Jumamyradov and Munkin (2022) use Halton draws.) According to Train,
"[t]he superior coverage and the negative correlation over observations that
are obtained with Halton draws combine to make Halton draws far more
effective than random [uniform] draws for simulation" (2003, p. 231). Train
also states that "[d]iscarding the initial part of the sequence eliminates the
correlation [between the initial terms of two Halton sequences]. The number of
initial elements to discard needs to be at least as large as the largest
prime that is used in creating the sequences (2003, p. 230)". By default,
burn(#) uses # = 5, noting that our sequences are created using the
first two primes (2, 3).

The default number of simulations is 10. This is a very low value and should
only be used for checking syntax. For serious modelling, the number of
simulations is likely to be much larger. We recommend that users, having
fitted a model with some specific number of simulations, increase the number
systematically to check that estimates of the mixture distribution parameters
in particular are robust.

Antithetic acceleration is a variance reduction technique according to
which negatively correlated draws are used instead of independently
distributed random numbers. We follow MT (section 4.3, pp. 38-39) and
create one antithetic draw for each draw by reversing the sign of both
of the two elements. To implement antithetic acceleration, use the
antithetic option. Antithetic acceleration is commonly used in
conjunction with standard uniform draws. (It is almost never used together
with Halton draws.)

MT (section 4.2, pp. 37-38) also implement a first-order method for correcting
asymptotic bias in the MSL estimator. To implement this method, use the
bias option.

On the basis of their Monte-Carlo analysis of the MSL estimator using the
sampling function approach, MT (p. 40) report that "[t]he first-order
correction of the asymptotic bias significantly reduces the number of
simulations needed for the estimation up a given degree of bias.
First-order corrected SML together with antithetic sampling produces
most efficient estimates." Jumamyradov and Munkin (2022) provide more
extensive analysis of potential bias, and how it varies with the true values
of the mixture distribution parameters.

Options

Main

nsim(#) sets number of simulations (default value is 10).
Unless checking syntax, users should normally specify a much
larger number

antithetic uses antithetic acceleration, as described above.
Given nsims(#), the total number of simulations is 2# if
this option is specified

halton uses Halton draws for simulation rather than the default which is
standard uniform random draws

burn(#) drop # initial elements when creating the Halton sequences
(default is 5)

bias uses the first-order bias correction, as described above, and is
recommended

seed(str) sets the seed. Relevant only when using uniform draws. In this
case, users are recommended to set the seed for reproducibility. If no seed
is set, the current rng state is used

rho_sf(#) sets the sampling function correlation (default # = 0)

alt uses the alternative MSL estimation method (default uses sampling
function

rng0 sets the random number generator (`rng') state back to its
pre-estimation value: see rngstate

Model

noconstant; see
[R] Estimation options.

offset1(varname), offset2(varname),
constraints(constraints); see
[R] Estimation options.

SE/Robust

Reporting

level(#), nocnsreport; see
[R] Estimation options.

Maximization

maximize_options: difficult,
technique(maximize##algorithm_spec:algorithm_spec),
iterate(#), [no]log, trace,
gradient, showstep, hessian,
showtolerance, tolerance(#), ltolerance(#),
nrtolerance(#), nonrtolerance, and
from(init_specs); see [R] Maximize.
These options are seldom used.

Setting the optimization type to technique(bhhh) resets the default
vcetype to vce(opg).

collinear, coeflegend; see [R] Estimation options.

Examples

The examples below are based on the data about numbers of doctor and hospital
visits and their predictors accompanying the bivcnto package
by Xu and Hardin (2016). The data are downloadable via the commands net sj 16-2 st0433 and
net get st0433. MT's empirical application is based on data previously
used by Deb and Trivedi (1997): see "debtrivedi.dta". (We downloaded these
data from the Journal of Applied Econometrics archive website
and converted them to Stata format.)

use "rwm1984.dta", clear

correlate docvis hospvis

generate byte postHS = edlevel1 == 0

generate byte MarM = (married == 1 & female == 0)

generate byte MarF = (married == 1 & female == 1)

generate byte SinM = (married == 0 & female == 0)

generate byte SinF = (married == 0 & female == 1)

global xvars MarM SinM SinF kids outwork postHS

* 400 uniform draws & antithetics, bias correction

bimpoisson docvis hospvis $xvars, robust ///

seed(13159) nsim(400) bias antithetic

matrix b1 = e(b)

* increase # sims (800)

bimpoisson docvis hospvis $xvars, robust ///

seed(4321) from(b1) nsim(800) bias antithetic

* same again, but with Halton draws

bimpoisson docvis hospvis $xvars, robust halton ///

from(b1) nsim(800) bias

* alternative approach (not using sampling function)

bimpoisson docvis hospvis $xvars, robust halton alt ///

nsim(400) bias

Stored results

In addition to the usual results saved after ml,
bimpoisson also saves the following in e():

Scalars

e(nsims) number of simulations used to calculate the
simulated log-likelihood

e(sigma1) estimate of sigma1

e(sigma1_se) estimate of the standard error of sigma1

e(sigma2) estimate of sigma1

e(sigma2_se) estimate of the standard error of sigma2

e(rho) estimate of rho

e(rho_se) estimate of the standard error of rho

e(antithetic) 1 if antithetic option used, 0 otherwise

e(halton) 1 if halton option used, 0 otherwise

e(burn) value of # set by burn(#) option

e(bias) 1 if bias option used, 0 otherwise

e(rho_sf) value of # set by rho_sf(#) option

e(rng0) 1 if rng0 option used; 0 otherwise

Macros

e(seed) contains the seed set by the seed option if
specified. If seed not set, contains "seed not set"

Authors

Stephen P. Jenkins <s.jenkins@lse.ac.uk>

London School of Economics and Political Science (LSE)

Fernando Rios-Avila <f.rios.a@gmail.com>

London School of Economics and Political Science (LSE)

References

Deb, P. and Trivedi, P. K. 1997. Demand for medical care by the elderly:
a finite mixture approach. Journal of Applied Econometrics 12, 313-336.
https://doi.org/10.1002/(SICI)1099-1255(199705)12:3%3C313::AID-JAE440%3E3.0.CO;2-G

Jumamyradov, M. and Munkin, M. K. 2022. Biases in maximum simulated
likelihood estimation of bivariate models. Journal of Econometric Methods
11, 55-70. https://doi.org/10.1515/jem-2021-0003

Munkin, M. K. and Trivedi, P. K. 1999. Simulated maximum likelihood estimation of multivariate mixed-Poisson regression
models, with application. The Econometrics Journal 2, 29-48.
https://doi.org/10.1111/1368-423X.00019

Train, K. E. 2003. Discrete Choice Methods with Simulation.
Cambridge: Cambridge University Press.
https://doi.org/10.1017/CBO9780511753930

Xu, X. and Hardin, J. W. 2016. Regression models for bivariate count outcomes.
The Stata Journal 2, 301-315.
https://doi.org/10.1177/1536867X1601600203

Acknowledgements

The syntax parsing code for bimpoisson borrows heavily from that
accompanying biprobit, and this help file is modelled on the help
file for biprobit.

Also see

poisson, bivcnto if installed (net sj 16-2 st0433),
ppml if installed (net sj 11-2 st0225),
ppmlhdfe if installed (net sj 20-1 st0589)