Estimate a Random Effects Negative Binomial regression model

Usage

renb(
  formula,
  group_var,
  data,
  method = "NM",
  max.iters = 1000,
  print.level = 0,
  bootstraps = NULL,
  offset = NULL
)

Arguments

formula: an R formula.
group_var: the grouping variable(s) for the random effects (e.g., individual ID or other panel ID variables).
data: a dataframe that has all of the variables in the formula.
method: a method to use for optimization in the maximum likelihood estimation. For options, see maxLik. Note that "BHHH" is not available for this function due to the implementation for the random effects.
max.iters: the maximum number of iterations to allow the optimization method to perform.
print.level: Integer specifying the verbosity of output during optimization.
bootstraps: Optional integer specifying the number of bootstrap samples to be used for estimating standard errors. If not specified, no bootstrapping is performed.
offset: an optional offset term provided as a string.

Value

An object of class countreg which is a list with the following components:

model: the fitted model object.
data: the data frame used to fit the model.
call: the matched call.
formula: the formula used to fit the model.

Details

This function estimates a random effects negative binomial (RENB) regression model. This model is based on the NB-1 model. The PDF for the RENB is: $$f(y_{it}|\lambda_{it}, a, b) = \frac{\Gamma(a+b) \Gamma(a + \sum_{t = 1}^{n_i} \\lambda_{it}) \Gamma(b + \sum_{t=1}^{n_i}y_{it})} {\Gamma(a) \Gamma(b) \Gamma(a + b + \sum_{t=1}^{n_i}\lambda_{it} + \sum_{t=1}^{n_i}y_{it})} \prod_{t=1}^{n_i} \frac{\Gamma(\lambda_{it}+y_{it})}{\Gamma(\lambda_{it})\Gamma(y_{it})}$$

Where $y_{it}$ is the count outcome for individual $i$ at time $t$, and $\lambda_{it}$ is the latent Poisson mean parameter for individual $i$ at time $t$. The parameters $a$ and $b$ are the shape parameters for the beta distribution that is used to model the random effects. The RENB model allows for overdispersion in the count data and accounts for unobserved heterogeneity across individuals by including random effects in the model. This formulation follows the approach described in the paper by Hausman, Hall, and Griliches (1984) for modeling panel data with random effects.

The marginal mean and marginal variance of the RENB model are given by: $$E[y_{it}] = \lambda_{it}\frac{b}{a-1}=\mu_it=exp(X_{it}\beta)$$

$$Var[y_{it}] = \frac{a+b-1}{a-2}\mu_{it} + \frac{a+b-1}{b(a-2)} \mu_{it}^2$$

Thus, the formulation of the model estimated here allows the use of the estimated coefficients to directly compute the marginal mean.

Note that the RENB model is a panel data model, and the group_var argument must be specified to indicate the grouping variable(s) for the random effects. The model is estimated using maximum likelihood estimation, and the optimization is performed using the maxLik package. The user can specify the optimization method and maximum iterations.

References

Hausman, Jerry A., Bronwyn H. Hall, and Zvi Griliches. "Econometric models for count data with an application to the patents–R&D relationship." Econometrica: Journal of the Econometric Society (1984): 909-938.

Examples

# \donttest{
## RENB Model
data("washington_roads")
washington_roads$AADTover10k <- 
  ifelse(washington_roads$AADT > 10000, 1, 0) # create a dummy variable
renb.mod <- renb(Animal ~ lnaadt + speed50 + ShouldWidth04 + AADTover10k,
                                data=washington_roads,
                                offset = "lnlength",
                                group_var="ID",
                                method="nm",
                                max.iters = 1000)
#> Warning: NaNs produced
summary(renb.mod)
#> Call:
#>  Animal ~ lnaadt + speed50 + ShouldWidth04 + AADTover10k 
#> 
#>  Method:  RENB 
#> Iterations:  736 
#> Convergence:  successful convergence  
#> Log-likelihood:  -263.7953 
#> 
#> Parameter Estimates:
#> # A tibble: 8 × 7
#>   parameter           coeff `Std. Err.` `t-stat` `p-value` `lower CI` `upper CI`
#>   <chr>               <dbl>       <dbl>    <dbl>     <dbl>      <dbl>      <dbl>
#> 1 (Intercept)        -9.35        1.22     -7.64     0        -11.7       -6.95 
#> 2 lnaadt              0.965       0.141     6.82     0          0.688      1.24 
#> 3 speed50            -0.988       0.34     -2.91     0.004     -1.65      -0.322
#> 4 ShouldWidth04      -0.415       0.283    -1.47     0.143     -0.97       0.14 
#> 5 AADTover10k        -0.868       0.511    -1.70     0.089     -1.87       0.133
#> 6 ln(a-1)             2.91       NA        NA       NA         NA         NA    
#> 7 ln(b)               0.619       0.309     2.00     0.045      0.013      1.23 
#> 8 lnlength (Offset …  1          NA        NA       NA         NA         NA    
# }