---
title: "Panel generalized linear models with panglm"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Panel generalized linear models with panglm}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4,
  eval = requireNamespace("plm", quietly = TRUE)
)
```

## Introduction

Longitudinal and panel data contain repeated observations for each sampling
unit. Dependence within units must therefore be represented either through
unit-specific effects or through an appropriate covariance estimator. The
`panglm` package provides a common interface for pooled, fixed-effects, and
random-effects generalized linear models for continuous, binary, and count
outcomes.

This vignette defines the models fitted by the package, describes their
estimating equations, and illustrates the associated inferential methods.
The emphasis is on the statistical specification of each estimator. The
numerical routines use Rcpp, RcppArmadillo, and RcppParallel.

```{r setup}
library(panglm)
data(Grunfeld, package = "plm")
```

## Notation and data organization

Let $y_{it}$ denote the response for unit $i$ at occasion $t$, and let
$x_{it}$ be a vector of observed covariates. The conditional mean is linked
to a linear predictor by

$$
g\{E(y_{it} \mid x_{it}, \alpha_i, \gamma_t)\}
= x_{it}^{\mathsf T}\beta + \alpha_i + \gamma_t,
$$

where $\alpha_i$ is an individual effect and $\gamma_t$ is an optional time
effect. The `index` argument identifies the individual and time variables.
The `model` argument determines how the individual effect is treated:

* `"pooling"` omits $\alpha_i$ from the mean model;
* `"within"` treats $\alpha_i$ as a fixed nuisance parameter;
* `"random"` assigns a distribution to $\alpha_i$ and integrates it out.

The optional `effect = "twoways"` includes both individual and time effects
for supported fixed-effects models.

Model matrices follow the usual R formula rules, including factors and
interactions. Before estimation, `panglm` checks rank after applying the
transformation associated with the requested estimator. Unidentified
columns, including time-invariant regressors in a within model, are reported
and omitted. Missing responses or covariates follow `na.action`; missing or
duplicate panel index pairs are errors. Nontrivial observation weights and
nonzero offsets are outside the current estimator definitions and are
rejected explicitly.

## Pooled models

For Gaussian, Poisson, and binomial responses, pooled models maximize the
ordinary generalized linear model likelihood by iteratively reweighted
least squares. The negative binomial model uses the NB2 variance function

$$
\operatorname{Var}(y_{it} \mid x_{it}) = \mu_{it} +
\frac{\mu_{it}^2}{\theta},
$$

and alternates coefficient updates with estimation of the shape parameter
$\theta$.

```{r pooled}
fit_pool <- panglm(
  inv ~ value + capital,
  data = Grunfeld,
  index = c("firm", "year"),
  model = "pooling",
  family = "gaussian"
)
coef(fit_pool)
```

Because pooling does not model within-unit dependence, inference should
usually use a covariance estimator that is robust to clustering by unit.

```{r pooled-cluster}
sqrt(diag(vcov(fit_pool, type = "cluster")))
```

## Fixed-effects models

### Gaussian response

The one-way Gaussian estimator applies the within transformation

$$
\widetilde y_{it} = y_{it} - \overline y_i, \qquad
\widetilde x_{it} = x_{it} - \overline x_i,
$$

and estimates $\beta$ by least squares using the transformed variables.
Time-invariant regressors are removed by this transformation and cannot be
identified. With `effect = "twoways"`, alternating projections remove both
individual and time means. This procedure also accommodates unbalanced
panels.

```{r within-gaussian}
fit_fe <- panglm(
  inv ~ value + capital,
  data = Grunfeld,
  index = c("firm", "year"),
  model = "within",
  family = "gaussian"
)
coef(fit_fe)
```

### Poisson response

For a one-way Poisson model, the individual intercept has a closed-form
profile solution. If $Y_i = \sum_t y_{it}$ and
$L_i(\beta) = \sum_t \exp(x_{it}^{\mathsf T}\beta)$, then

$$
\widehat\alpha_i(\beta) = \log\{Y_i/L_i(\beta)\}.
$$

Substitution into the likelihood gives the conditional Poisson estimator.
For two-way Poisson models, `panglm` uses an outer iteratively reweighted
least-squares loop and weighted alternating projections to absorb the two
sets of fixed effects.

### Negative binomial response

The one-way fixed-effects NB2 estimator jointly estimates the covariate
coefficients, one intercept for each individual, and the shape parameter.
Groups for which every count is zero have an unbounded intercept and are
excluded with an informative message. The two-way model absorbs individual
and time effects through weighted alternating projections while updating
the NB2 shape parameter.

The one-way specification is an unconditional nonlinear fixed-effects
estimator. It is distinct from the Hausman-Hall-Griliches conditional
negative binomial model and avoids interpreting that model's scale
restriction as an individual intercept. Explicit dummy intercepts do not,
however, remove the general finite-T incidental-parameter concern. Results
from short panels should therefore be interpreted with appropriate caution.

### Binomial response

The one-way binomial estimator is the exact conditional logistic model. For
unit $i$, inference is conditional on the observed number of successes
$S_i = \sum_t y_{it}$. The contribution to the conditional likelihood is

$$
L_i(\beta \mid S_i) =
\frac{\exp\left(\sum_t y_{it}x_{it}^{\mathsf T}\beta\right)}
{\sum_{a:\,\sum_t a_t=S_i}
 \exp\left(\sum_t a_tx_{it}^{\mathsf T}\beta\right)}.
$$

The denominator and score are evaluated by a forward and backward dynamic
program. Units with all-zero or all-one responses do not contribute to this
conditional likelihood and are excluded.

## Random-effects models

Random-effects models retain time-invariant regressors and account for
within-unit dependence by integrating over a latent individual effect. The
implemented specifications are:

* Gaussian responses use the Swamy-Arora feasible generalized least-squares
  estimator;
* Poisson responses use a Poisson-Gamma mixture, which yields a closed-form
  marginal likelihood;
* negative binomial responses use a two-parameter beta-negative-binomial
  marginal likelihood;
* binomial responses use a Gaussian random intercept integrated by
  Gauss-Hermite quadrature.

The number of quadrature nodes for the binomial model is controlled by `R`.
Random-effects covariance matrices are based on the observed information.

```{r random-gaussian}
fit_re <- panglm(
  inv ~ value + capital,
  data = Grunfeld,
  index = c("firm", "year"),
  model = "random",
  family = "gaussian"
)
coef(fit_re)
```

## Covariance estimation and confidence intervals

The default covariance matrix is model based. Pooled and fixed-effects
models also support `type = "HC1"` and `type = "cluster"`. Let $s_j(\beta)$
be a score contribution and let $B$ be the inverse information matrix. The
sandwich estimator has the form

$$
\widehat V(\widehat\beta) = B
\left(\sum_j s_j(\widehat\beta)s_j(\widehat\beta)^{\mathsf T}\right)B.
$$

For cluster-robust inference, observation scores are first summed within
clusters. The panel individual is used as the default cluster. For exact
conditional binomial models, the independent likelihood contribution is a
complete panel stratum. HC1 therefore uses stratum-level scores, and a
custom clustering variable may combine but cannot split strata. Gaussian,
Poisson, negative binomial, and conditional binomial fixed-effects fits all
provide the score information required by these estimators.

```{r inference}
vcov(fit_fe, type = "cluster")
confint(fit_fe)
```

## Model assessment

For models that store fitted values, `plot()` provides coefficient,
residual, and observed-versus-fitted displays. Coefficient plots are
available for every `panglm` fit.

`fitted()` and in-sample `predict()` return values in the original data row
order. Random-effects response predictions are marginal over the fitted
latent-effect distribution. Fixed-effects predictions for new data include
estimated effects for known panel levels; unseen levels are rejected unless
`allow.new.levels = TRUE`. For exact conditional logistic regression,
in-sample fitted values are probabilities conditional on each unit's
observed success total. Arbitrary new-data probabilities are not identified
because the individual intercept has been conditioned out.

```{r plots, fig.alt = "Coefficient estimates with confidence intervals"}
plot(fit_fe, which = "coefficients")
```

`panglm_dispersiontest()` reports the Pearson statistic divided by its
residual degrees of freedom for Poisson and negative binomial fits.
`panglm_hausman()` compares compatible fixed-effects and random-effects
estimates using their covariance difference.

For likelihood-based estimators, `logLik()`, `AIC()`, and `BIC()` include
estimated shape, dispersion, variance, and absorbed fixed-effect parameters
when these parameters belong to the reported likelihood. Conditional
Poisson and conditional logistic likelihoods count only the regression
parameters that remain after conditioning. Gaussian within and
Swamy-Arora random-effects estimators do not report a likelihood, so their
information criteria are undefined.

```{r hausman}
panglm_hausman(fit_fe, fit_re)
```

## Hurdle models for panel counts

For outcomes with a distinct zero-generating process, `panglm_hurdle()`
fits two fixed-effects components. The first is an exact conditional
logistic model for $I(y_{it}>0)$. The second is a zero-truncated Poisson or
NB2 model for positive counts. Its density is

$$
f_+(y \mid \mu, \theta) =
\frac{f(y \mid \mu, \theta)}{1-f(0 \mid \mu, \theta)}, \qquad y>0.
$$

The count distribution is selected with `count_family`. Both positive-count
models jointly estimate covariate coefficients and individual intercepts;
the NB2 model also estimates $\theta$.

```{r hurdle}
data(copd)
fit_hurdle <- panglm_hurdle(
  exacerbations ~ crp,
  data = copd,
  index = c("id", "visit"),
  count_family = "negbin"
)
plot(fit_hurdle)
```

## Computational validation

The test suite compares estimators with established implementations where
the statistical models coincide. Comparisons include `stats::glm()` for
pooled generalized linear models, `MASS::glm.nb()` for pooled NB2,
`plm::plm()` for Gaussian panel models, `survival::clogit()` for exact
conditional logistic regression, `fixest` for absorbed fixed-effects count
models, `pglm` for compatible random-effects likelihoods, and
`pscl::hurdle()` for zero-truncated count likelihoods.

## Model scope

The package focuses on Gaussian, Poisson, binomial, and NB2 panel models.
It does not provide ordinal, Tobit, or between estimators. A general
two-way exact conditional logistic likelihood is not available, so the
binomial fixed-effects model is restricted to individual effects. These
restrictions define the current statistical scope and prevent unlike
estimators from being presented under a common option.

## References

Allison, P. D., and Waterman, R. P. (2002). Fixed-effects negative binomial
regression models. *Sociological Methodology*, 32, 247-265.

Chamberlain, G. (1980). Analysis of covariance with qualitative data.
*Review of Economic Studies*, 47, 225-238.

Guimaraes, P., and Portugal, P. (2010). A simple feasible procedure to fit
models with high-dimensional fixed effects. *Stata Journal*, 10, 628-649.

Hausman, J., Hall, B. H., and Griliches, Z. (1984). Econometric models for
count data with an application to the patents-R&D relationship.
*Econometrica*, 52, 909-938.