### Load Libraries and Set Working Directory in R

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.md

Initial setup for R scripts, including loading necessary libraries and setting the working directory. Ensure packages are installed before loading.

```r
rm(list=ls())
#install.packages("glm2")
#install.packages("alpaca")
#install.packages("FENmlm")
library(alpaca)
library(glm2)
library(FENmlm)
setwd("C:/Git/ppml_hdfe_demo/guides")
```

--------------------------------

### Install ppmlhdfe from GitHub

Source: https://context7.com/sergiocorreia/ppmlhdfe/llms.txt

Installs the latest version of ppmlhdfe and its dependencies from GitHub. Compiles Mata libraries and verifies the installation.

```stata
* Install from GitHub (latest version)
net install ftools, from("https://raw.githubusercontent.com/sergiocorreia/ftools/master/src/")
net install reghdfe, from("https://raw.githubusercontent.com/sergiocorreia/reghdfe/master/src/")
net install ppmlhdfe, from("https://raw.githubusercontent.com/sergiocorreia/ppmlhdfe/master/src/")

ftools, compile
reghdfe, compile
ppmlhdfe, version
```

--------------------------------

### Install ppmlhdfe from GitHub

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/README.md

Installs the latest versions of ftools, reghdfe, and ppmlhdfe directly from GitHub. This method requires uninstalling existing versions and then installing from the specified GitHub URLs. It also includes commands to compile the packages and check their versions.

```stata
* Install ftools
cap ado uninstall ftools
net install ftools, from("https://raw.githubusercontent.com/sergiocorreia/ftools/master/src/")

* Install reghdfe
cap ado uninstall reghdfe
net install reghdfe, from("https://raw.githubusercontent.com/sergiocorreia/reghdfe/master/src/")

* Install ppmlhdfe
cap ado uninstall ppmlhdfe
net install ppmlhdfe, from("https://raw.githubusercontent.com/sergiocorreia/ppmlhdfe/master/src/")

* Create compiled files
ftools, compile
reghdfe, compile

* Check versions
ppmlhdfe, version

* Clear programs already in memory
program drop _all

* Test program
sysuse auto, clear
reghdfe price weight, a(turn)
ppmlhdfe price weight, a(turn)
```

--------------------------------

### Install ppmlhdfe from SSC

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/README.md

Installs stable versions of ftools, reghdfe, and ppmlhdfe from SSC. Includes commands to uninstall previous versions and compile the installed packages. A test program is provided to verify installation.

```stata
cap ado uninstall ftools
cap ado uninstall reghdfe
cap ado uninstall ppmlhdfe

ssc install ftools
ssc install reghdfe
ssc install ppmlhdfe

clear all
ftools, compile
reghdfe, compile

* Test program
sysuse auto, clear
reghdfe price weight, a(turn)
ppmlhdfe price weight, a(turn)
```

--------------------------------

### Install ppmlhdfe from SSC

Source: https://context7.com/sergiocorreia/ppmlhdfe/llms.txt

Installs the stable version of ppmlhdfe and its dependencies (ftools, reghdfe) from SSC. Compiles Mata libraries and verifies the installation.

```stata
* Install from SSC (stable version)
ssc install ftools
ssc install reghdfe
ssc install ppmlhdfe

* Compile Mata libraries
ftools, compile
reghdfe, compile

* Verify installation
ppmlhdfe, version
```

--------------------------------

### glm2 Package: Poisson Regression on Example 2 Data with Convergence Issues

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.md

Fits a Poisson regression model using the `glm2` package with data from 'example2.csv'. This example may also produce convergence warnings, similar to the base R `glm` function.

```r
rm(list=ls())
data <- read.csv(file="csv/example2.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 + x3 + x4
mod <- glm2(formula, data, family=poisson())
summary(mod)
```

--------------------------------

### Test GLM Convergence with Poisson Regression

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_benchmarks.md

Executes a Poisson regression using the glm package. This example demonstrates a case where the algorithm fails to converge and issues warning messages.

```R
>         data <- read.csv(file="example2.csv", header=TRUE, sep=",")
>         formula <- y ~ x1 + x2 + x3 + x4
>         mod <- glm(formula, data, family=poisson())
Warning messages:
1: glm.fit: algorithm did not converge 
2: glm.fit: fitted rates numerically 0 occurred 
>         summary(mod)

Call:
glm(formula = formula, family = poisson(), data = data)

Deviance Residuals: 
       Min          1Q      Median          3Q         Max  
-3.109e-03  -2.000e-08  -2.000e-08  -2.000e-08   2.000e-08  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -367.83    7448.51  -0.049    0.961
x1            512.42   10345.15   0.050    0.960
x2          -1644.86   33284.96  -0.049    0.961
x3           -105.85    2138.00  -0.050    0.961
x4             20.56     413.81   0.050    0.960

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 4.6719e+01  on 11  degrees of freedom
Residual deviance: 9.6716e-06  on  7  degrees of freedom
AIC: 19.933

Number of Fisher Scoring iterations: 25
```

--------------------------------

### Handle Errors in feglm with try

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.md

This example shows how to use `try` to catch potential errors when fitting a model with `feglm`, such as when the optimization algorithm fails to converge. The `try` function prevents the script from stopping if an error occurs.

```r
rm(list=ls())
data <- read.csv(file="csv/fe2.csv", header=TRUE, sep=",")
formula <- y ~ x1 | i + j
try(mod <- feglm(formula, data, family=poisson()))
```

--------------------------------

### Fit Poisson GLM with glm2

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_benchmarks.md

Examples of fitting a Poisson regression model using the glm2 package on different datasets.

```R
>         data <- read.csv(file="csv/example2.csv", header=TRUE, sep=",")
>         formula <- y ~ x1 + x2 + x3 + x4
>         mod <- glm2(formula, data, family=poisson())
Warning messages:
1: glm.fit2: algorithm did not converge. Try increasing the maximum iterations 
2: glm.fit2: fitted rates numerically 0 occurred 
>         summary(mod)

Call:
glm2(formula = formula, family = poisson(), data = data)

Deviance Residuals: 
       Min          1Q      Median          3Q         Max  
-3.109e-03  -2.000e-08  -2.000e-08  -2.000e-08   2.000e-08  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -367.83    7448.51  -0.049    0.961
x1            512.42   10345.15   0.050    0.960
x2          -1644.86   33284.96  -0.049    0.961
x3           -105.85    2138.00  -0.050    0.961
x4             20.56     413.81   0.050    0.960

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 4.6719e+01  on 11  degrees of freedom
Residual deviance: 9.6716e-06  on  7  degrees of freedom
AIC: 19.933

Number of Fisher Scoring iterations: 25
```

```R
>         data <- read.csv(file="csv/example1.csv", header=TRUE, sep=",")
>         formula <- y ~ x1 + x2 + x3 + x4
>         mod <- glm2(formula, data, family=poisson())
>         summary(mod)

Call:
glm2(formula = formula, family = poisson(), data = data)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.97355  -0.75131  -0.16879   0.07357   2.70708  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept)    0.59095    0.30291   1.951   0.0511 .
x1           -17.78017 3467.85856  -0.005   0.9959  
x2            17.32952 3467.85857   0.005   0.9960  
x3            -0.47085    0.23117  -2.037   0.0417 *
x4            -0.03779    0.04375  -0.864   0.3878  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 31.912  on 11  degrees of freedom
Residual deviance: 15.956  on  7  degrees of freedom
AIC: 46.991

Number of Fisher Scoring iterations: 15
```

--------------------------------

### Use Simplex Solver with ppmlhdfe

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/undocumented.md

Directly call the Simplex solver in ppmlhdfe to verify results. This example passes a matrix of observations and regressors to the simplex option.

```stata
ppmlhdfe, simplex( (1, -1, -1, 0 \ -1, 1, -1, 0 \ -1, -1, 1, 0)' ) ans(0, 0, 0, 1)
```

--------------------------------

### Base R glm: Poisson Regression on Example 2 Data with Convergence Issues

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.md

Fits a Poisson regression model using base R's `glm` function with data from 'example2.csv'. This example may produce convergence warnings, indicating potential numerical instability.

```r
rm(list=ls())
data <- read.csv(file="csv/example2.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 + x3 + x4
mod <- glm(formula, data, family=poisson())
summary(mod)
```

--------------------------------

### Run glm2 GLM

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.nb.html

Executes a Poisson GLM using the glm2 package on different example datasets.

```R
rm(list=ls())
data <- read.csv(file="csv/example1.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 + x3 + x4
mod <- glm2(formula, data, family=poisson())
summary(mod)
```

```R
rm(list=ls())
data <- read.csv(file="csv/example2.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 + x3 + x4
mod <- glm2(formula, data, family=poisson())
summary(mod)
```

--------------------------------

### sep(mu) Fragility with Skewed Data

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_primer.md

Illustrates the fragility of the `sep(mu)` method with skewed data. A lower `mu_tol` can cause it to fail detection. This example shows how a very small `mu_tol` might not detect separation.

```stata
clear
input double(y id1 id2)
0 1 1
1 1 1
0 2 1
0 2 2
1 2 2
1e-6 2 1
1e-6 2 1
1e-6 2 1
end

ppmlhdfe y, a(id1 id2) sep(mu) mu_tol(1e-6)
```

--------------------------------

### glm2 Package: Poisson Regression on Example 1 Data

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.md

Fits a Poisson regression model using the `glm2` package with data from 'example1.csv'. `glm2` is an alternative implementation that may offer performance benefits or different convergence properties.

```r
rm(list=ls())
data <- read.csv(file="csv/example1.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 + x3 + x4
mod <- glm2(formula, data, family=poisson())
summary(mod)
```

--------------------------------

### Fit Poisson GLM with glm2

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.nb.html

Fits a Poisson regression model using the glm2 package. The first example demonstrates a standard fit, while the second shows handling of non-convergence.

```R
rm(list=ls())
data <- read.csv(file="csv/example1.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 + x3 + x4
mod <- glm2(formula, data, family=poisson())
summary(mod)
```

```R
rm(list=ls())
data <- read.csv(file="csv/example2.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 + x3 + x4
mod <- glm2(formula, data, family=poisson())

glm.fit2: algorithm did not converge. Try increasing the maximum iterationsglm.fit2: fitted rates numerically 0 occurred

summary(mod)
```

--------------------------------

### Generate Regression Tables with esttab and margins

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/undocumented.md

Use this code to generate journal-style regression tables by combining ppmlhdfe with estpost margins and esttab. Ensure esttab and margins are installed.

```stata
cls
estimates clear 
sysuse auto, clear

qui ppmlhdfe price weight, a(turn) d
qui estpost margins, dydx(_all)
qui eststo

qui ppmlhdfe price weight length, a(turn trunk) d
qui estpost margins, dydx(_all)
qui eststo

estfe *, labels(turn "Turn FE" trunk "Trunk FE")
esttab, indicate(`r(indicate_fe)', labels("yes" "")) stat(N ll, fmt(%9.0fc %10.1fc)) se starlevels(* 0.1 ** 0.05 *** 0.01) label compress
```

--------------------------------

### Fixed Effects with Individual Slopes (Full Interaction)

Source: https://context7.com/sergiocorreia/ppmlhdfe/llms.txt

This example demonstrates absorbing fixed effects and heterogeneous slopes using the `##` operator, specifically modeling a full journal-year interaction (`issn#year`). This allows for distinct trends for each journal across all years.

```stata
* Load citation data (included with package examples)
use "examples/citations_example", clear

* Model 3: Full journal-year interaction (issn#year)
ppmlhdfe cit nbaut, absorb(issn#year type jel2 pubyear)
```

--------------------------------

### Base R glm: Poisson Regression on Example 1 Data

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.md

Fits a Poisson regression model using base R's `glm` function with data from 'example1.csv'. This is suitable for standard Poisson regression tasks.

```r
rm(list=ls())
data <- read.csv(file="csv/example1.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 + x3 + x4
mod <- glm(formula, data, family=poisson())
summary(mod)
```

--------------------------------

### Aggressive sep(mu) and Incorrect Dropping of Observations

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_primer.md

Shows how an aggressive `mu_tol` in `sep(mu)` can lead to slow convergence and erroneously drop non-separated observations. This example uses a high `mu_tol` which causes convergence issues and incorrect dropping.

```stata
clear
input double(y id1 id2)
0 1 1
1 1 1
0 2 1
0 2 2
1 2 2
1e-6 2 1
1e-6 2 1
1e-6 2 1
end

ppmlhdfe y, a(id1 id2) sep(mu) mu_tol(1e-2)
```

--------------------------------

### Create and Regress Data with PPMLHDfe

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_primer.md

This snippet demonstrates how to create sample data and run a PPMLHDfe regression. It highlights how PPMLHDfe handles cases that might lead to infinite estimates by dropping a variable.

```stata
* Create data
clear
set obs 6
gen y = max(0, _n-3)
gen x1 = 2*(_n==1) - (_n==2) + cond(_n>4, _n, 0)
gen x2 = 2 * (_n==2) - 2 * cond(_n>4, _n, 0) - (_n==1)

ppmlhdfe y x1 x2
```

--------------------------------

### Compare Poisson regression methods for separation

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_benchmarks.md

Demonstrates how different Stata commands handle data separation by comparing ppmlhdfe, poisson, and glm implementations.

```stata
loc url "https://raw.githubusercontent.com/sergiocorreia/ppmlhdfe/master/guides/csv"
insheet using "`url'/example1.csv", clear

ppmlhdfe y x*
poisson y x*, iter(100) nolog
cap noi glm y x*, family(poisson) ml iter(100) nolog
cap noi glm y x*, family(poisson) irls iter(100) nolog
ppml y x*
ppml y x*, strict
```

--------------------------------

### Run GLM regression in Julia

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_benchmarks.md

Demonstrates fitting a Generalized Linear Model in Julia using the GLM package.

```Julia
import Pkg
#Pkg.add("CSV")
#Pkg.add("DataFrames")
#Pkg.add("GLM")
#Pkg.add("CategoricalArrays")

using CSV, DataFrames, GLM
data = CSV.File("./csv/example1.csv") |> DataFrame
ans = fit(GeneralizedLinearModel, @formula(y ~ x1 + x2 + x3 + x4), data, Poisson())
round(deviance(ans), digits=5)
```

```Julia
import Pkg
using CSV, DataFrames, GLM
data = CSV.File("./csv/example2.csv") |> DataFrame
ans = fit(GeneralizedLinearModel, @formula(y ~ x1 + x2 + x3 + x4), data, Poisson())
round(deviance(ans), digits=5)
```

```Julia
import Pkg
using CSV, DataFrames, GLM, CategoricalArrays
data = CSV.File("csv/fe1.csv") |> DataFrame
data.id1 = CategoricalArray(data.i)
data.id2 = CategoricalArray(data.j)
ans = fit(GeneralizedLinearModel, @formula(y ~ x1 + x2 + id1 + id2), data, Poisson())
round(deviance(ans), digits=5)
```

--------------------------------

### Compare Poisson Regression Convergence

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_benchmarks.md

A comparison of various Stata commands for Poisson regression on a dataset prone to convergence issues.

```stata
clear
input int(y x1 x2 x3)
 0   0   1  0
 0   0   0  0
 0   0   0  0
 0   0   0  0
 0   1   9  0
 2  21  21 21
 3   0   0  0
 5   0   0  0
 7   0   0  0
10 -18 -18  0
end


poisson y x*, iter(1000) // does not converge, but deviance and LL are close to correct values (see figure)
glm y x*, family(poisson) ml iter(1000) // does not converge, but deviance and LL are close to correct values
glm y x*, family(poisson) irls iter(1000) // does not converge, but deviance and LL are close to correct values
ppml y x* // does not converge, but deviance and LL are close to correct values
ppml y x*, strict // converges to incorrect solution; with b[x3] = 0.961 and deviance and log-likelihoods far from the correct values
ppmlhdfe y x*, sep(none) maxiter(100) // If we disable separation checks, ppmlhdfe will also fail to converge

ppmlhdfe y x* // only command that converges to correct estimate; b[x3] = 0.066
```

--------------------------------

### Comparison of sep(ir) and sep(mu) for Correct Convergence

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_primer.md

Compares `sep(ir)` and `sep(mu)` for correct convergence. `sep(ir)` converges quickly and to the correct number of observations, while `sep(mu)` with a conservative tolerance also achieves the correct result.

```stata
ppmlhdfe y, a(id1 id2) sep(ir)
ppmlhdfe y, a(id1 id2) sep(mu) mu_tol(1e-6)
```

--------------------------------

### Fixed Effects with Individual Slopes (Journal-Specific Trends)

Source: https://context7.com/sergiocorreia/ppmlhdfe/llms.txt

Use the `##` operator to absorb both fixed effects and heterogeneous slopes, such as individual-specific trends. This example models journal-specific year trends by including `issn##c.year` in the `absorb()` option.

```stata
* Load citation data (included with package examples)
use "examples/citations_example", clear

* Model 2: Journal-specific year trends (issn#c.year)
ppmlhdfe cit nbaut, absorb(issn##c.year type jel2 pubyear)
```

--------------------------------

### Detect Separation with ppml

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_examples.md

Using the ppml command with the check option to identify potential existence issues.

```stata
. ppml y x z, check

note: checking the existence of the estimates

Number of regressors excluded to ensure that the estimates exist: 1
Excluded regressors:  z
Number of observations excluded: 2
```

--------------------------------

### ppmlhdfe Configuration Options

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/undocumented.md

Overview of the configuration parameters for the ppmlhdfe command, including IRLS, separation, and display options.

```APIDOC
## ppmlhdfe Configuration Options

### Description
Configuration options for the ppmlhdfe command, categorized by algorithm type and display settings.

### IRLS Options
- **tolerance** (float) - Convergence threshold (default 1e-8).
- **use_exact_solver** (int) - Whether to always use an exact least-squares solver (default 0).
- **use_exact_partial** (int) - Whether to partial out (z, X) every iteration (default 0).
- **target_inner_tol** (float) - Inner tolerance (default 1e-9).
- **start_inner_tol** (float) - Initial tolerance when partialling-out (default 1e-4).
- **min_ok** (int) - Number of times epsilon < tolerance must be observed (default 1).
- **maxiter** (int) - Maximum number of iterations (default 1000).
- **use_heuristic_tol** (int) - Anticipate convergence to increase inner tolerance (default 1).

### Simplex Separation Options
- **simplex_tol** (float) - Internal tolerance for simplex step (default 1e-12).
- **simplex_maxiter** (int) - Maximum iterations for simplex step (default 1000).

### IR/ReLU Separation Options
- **relu_tol** (float) - Internal tolerance (default 1e-4).
- **relu_zero_tol** (float) - Tolerance for converting values to zero (default 1e-8).
- **relu_maxiter** (int) - Maximum iterations for ReLU step (default 100).
- **relu_strict** (int) - Raise error if maxiter exceeded (default 0).
- **relu_accelerate** (int) - Increase weights of negative observations (default 0).

### Certificate of Separation Options
- **tagsep(varname)** - Save indicator variable for separated observations.
- **zvarname(varname)** - Save the certificate of separation.
- **r2** - Run least-squares regression between certificate z and regressors.

### Mu Separation Options
- **mu_tol** (float) - Criteria for tagging observation as separated (default 1e-6).
```

--------------------------------

### Estimate Poisson Regression with Tagged Separation in Stata

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_examples.md

This Stata command estimates a Poisson regression with fixed effects, similar to the previous example, but additionally tags separated observations using 'sep_relu' and includes 'zvar(z)' for potential further analysis. The output indicates perfect R-squared and near-zero residuals, suggesting a potential issue or specific data characteristic.

```stata
. ppmlhdfe y x*, a(id1 id2) sep(relu) tagsep(sep_relu) r2 zvar(z) 
<some output omitted...>
HDFE Linear regression                            Number of obs   =        100
Absorbing 2 HDFE groups                           F(   2,     79) =          .
                                                  Prob > F        =          .
                                                  R-squared       =     1.0000
                                                  Adj R-squared   =     1.0000
                                                  Within R-sq.    =     1.0000
                                                  Root MSE        =     0.0000

------------------------------------------------------------------------------
           z | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
          x1 |   3.25e-11   4.20e-11     0.77   0.442    -5.11e-11    1.16e-10
          x2 |          1   2.88e-10  3.5e+09   0.000            1           1
       _cons |  -9.37e-13   3.72e-11    -0.03   0.980    -7.49e-11    7.30e-11
<some output omitted...>
```

--------------------------------

### Creating Publication Tables with esttab and ppmlhdfe

Source: https://context7.com/sergiocorreia/ppmlhdfe/llms.txt

Illustrates generating formatted regression tables using `estout`/`esttab` with margins for marginal effects after estimating models with `ppmlhdfe`. This includes creating tables with fixed effect indicators.

```stata
estimates clear
sysuse auto, clear

* Estimation 1: Single FE
qui ppmlhdfe price weight, absorb(turn) d
qui estpost margins, dydx(_all)
qui eststo

* Estimation 2: Two FEs
qui ppmlhdfe price weight length, absorb(turn trunk) d
qui estpost margins, dydx(_all)
qui eststo

* Create table with FE indicators
estfe *, labels(turn "Turn FE" trunk "Trunk FE")
esttab, indicate(`r(indicate_fe)', labels("yes" "")) \
    stat(N ll, fmt(%9.0fc %10.1fc)) se \
    starlevels(* 0.1 ** 0.05 *** 0.01) label compress
/*
--------------------------------------------
                      (1)             (2)
--------------------------------------------
weight              2.419***        1.704***
                  (0.437)         (0.511)
length                              3.508**
                                  (1.614)
--------------------------------------------
Turn FE               yes             yes
Trunk FE               no             yes
--------------------------------------------
N                      74              74
ll                 -534.6          -523.9
--------------------------------------------
*/
```

--------------------------------

### Using Exposure and Offset Variables in ppmlhdfe

Source: https://context7.com/sergiocorreia/ppmlhdfe/llms.txt

Demonstrates how to include exposure variables (with logged coefficient constrained to 1) or offset variables for rate models in `ppmlhdfe`. Ensure the correct variable is specified for exposure or offset.

```stata
* Load ships data
use "https://www.stata-press.com/data/r16/ships", clear

* Exposure: include ln(service) with coefficient = 1
ppmlhdfe accident op_75_79 co_65_69 co_70_74 co_75_79, absorb(ship) exposure(service)

* Offset: include variable directly with coefficient = 1
gen log_service = ln(service)
ppmlhdfe accident op_75_79 co_65_69 co_70_74 co_75_79, absorb(ship) offset(log_service)

* Report incidence rate ratios instead of coefficients
ppmlhdfe accident op_75_79 co_65_69 co_70_74 co_75_79, absorb(ship) exposure(service) irr
```

--------------------------------

### Implement iterative rectifier algorithm manually

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_primer.md

A manual implementation of the IR algorithm using standard Stata code.

```stata
* Create data
clear
input byte(y id1 id2)
0 1 1
1 1 1
0 2 1
0 2 2
1 2 2
end

* Run IR (iterative rectifier) algorithm
loc tol = 1e-5
gen u =  !y
su u, mean
loc K = ceil(r(sum) / `tol' ^ 2)
gen w = cond(y, `K', 1) 

while 1 {
	qui reghdfe u [fw=w], absorb(id1 id2) resid(e)
	predict double xb, xbd
	qui replace xb = 0 if abs(xb) < `tol'

	* Stop once all predicted values become non-negative
	qui cou if xb < 0
	if !r(N) {
		continue, break
	}

	replace u = max(xb, 0)
	drop xb w
}

rename xb z
gen is_sep = z > 0
list y id1 id2 is_sep
```

--------------------------------

### Summarize FENmlm Model Results

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.md

After fitting a model with `FENmlm::femlm`, use `summary(mod)` to view the estimation results. Pay attention to convergence warnings, as they indicate that the results may not be reliable.

```r
summary(mod)
```

--------------------------------

### Tuning Convergence and Tolerance

Source: https://context7.com/sergiocorreia/ppmlhdfe/llms.txt

Adjusts algorithm parameters like tolerance, iteration limits, and initial values for complex models.

```stata
* Stricter tolerance (slower but more precise)
ppmlhdfe trade fta, absorb(imp#year exp#year imp#exp) tolerance(1e-10)

* Looser inner tolerance for faster iteration (less precise intermediate steps)
ppmlhdfe trade fta, absorb(imp#year exp#year imp#exp) itol(1e-6)

* Increase maximum iterations for slow-converging models
ppmlhdfe trade fta, absorb(imp#year exp#year imp#exp) maxiter(2000)

* OLS-based initial values (sometimes helps convergence)
ppmlhdfe trade fta, absorb(imp#year exp#year imp#exp) guess(ols)

* Suppress iteration log
ppmlhdfe trade fta, absorb(imp#year exp#year imp#exp) nolog

* Verbose output for debugging
ppmlhdfe trade fta, absorb(imp#year exp#year imp#exp) verbose(1)
```

--------------------------------

### Handle Errors in FENmlm::femlm with try

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.md

This code uses `try` to gracefully handle errors that may occur during model fitting with `FENmlm::femlm`. This is useful for preventing script termination due to issues like zero outcomes in clusters.

```r
rm(list=ls())
data <- read.csv(file="csv/fe2.csv", header=TRUE, sep=",")
formula <- y ~ x1 | i + j
try(mod <- FENmlm::femlm(formula, data, family="poisson"))
```

--------------------------------

### Sample CSV Data for Poisson Regression

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_examples.md

This is a sample dataset in CSV format, including outcome variable 'y', regressors 'x1' and 'x2', categorical fixed effects 'id1' and 'id2', and a 'separated' indicator. Used for testing Poisson regressions.

```csv
y,x1,x2,id1,id2,separated
0.0000000000,-0.9303550124,1,1,4,1
0.0000000000,0.1835959703,1,2,1,1
0.0000000000,-0.6371972561,0,2,6,0
0.0000000000,-0.4237562418,0,2,7,0
0.1527670026,-1.1799178123,0,8,4,0
0.1553160399,0.8860545158,0,1,7,0
0.1734523475,1.0502026081,0,8,3,0
0.2217264324,-0.2490162849,0,9,1,0
0.2260344625,0.9635434151,0,7,6,0
<90 rows omitted>
```

--------------------------------

### Detect separation with simplex method

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_primer.md

Demonstrates the simplex method for separation detection, noting its limitation with fixed effects.

```stata
clear
input byte(y id1 id2)
0 1 1
1 1 1
0 2 1
0 2 2
1 2 2
end

ppmlhdfe y, a(id1 id2) sep(fe simplex)
```

--------------------------------

### Run alpaca GLM

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.nb.html

Executes a Poisson GLM with fixed effects using the alpaca package.

```R
rm(list=ls())
data <- read.csv(file="csv/fe1.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 | i + j
mod <- feglm(formula, data, family=poisson())
summary(mod)
```

```R
rm(list=ls())
data <- read.csv(file="csv/fe2.csv", header=TRUE, sep=",")
formula <- y ~ x1 | i + j
try(mod <- feglm(formula, data, family=poisson()))
#summary(mod)
```

--------------------------------

### Analyze Endometrial Cancer Data with Stata Logit

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_examples.md

Demonstrates the use of tabulate and logit to identify quasi-complete separation in a dataset where a regressor perfectly predicts the outcome.

```stata
. tab y x1

           |          x1
         y |         0          1 |     Total
-----------+----------------------+----------
         0 |        49          0 |        49 
         1 |        17         13 |        30 
-----------+----------------------+----------
     Total |        66         13 |        79 

. logit y x1 x2 x3

note: x1 != 0 predicts success perfectly;
      x1 omitted and 13 obs not used.

Iteration 0:  Log likelihood = -37.653392  
<some output omitted...>

Logistic regression                                     Number of obs =     66
                                                        LR chi2(2)    =  19.91
                                                        Prob > chi2   = 0.0000
Log likelihood = -27.69663                              Pseudo R2     = 0.2644

------------------------------------------------------------------------------
           y | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
          x1 |          0  (omitted)
          x2 |  -.0421834    .044332    -0.95   0.341    -.1290725    .0447057
          x3 |  -2.902606   .8455516    -3.43   0.001    -4.559856   -1.245355
       _cons |   4.304518   1.637299     2.63   0.009     1.095471    7.513564
------------------------------------------------------------------------------
```

--------------------------------

### Fit Poisson Model with Fixed Effects using alpaca

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.html

Fits a Poisson model with fixed effects using the feglm function from the alpaca package.

```R
rm(list=ls())
data <- read.csv(file="csv/fe1.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 | i + j
mod <- feglm(formula, data, family=poisson())
summary(mod)
```

```R
rm(list=ls())
data <- read.csv(file="csv/fe2.csv", header=TRUE, sep=",")
formula <- y ~ x1 | i + j
try(mod <- feglm(formula, data, family=poisson()))
```

--------------------------------

### Fit Poisson Fixed Effects Model with FENmlm::femlm

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.md

This snippet demonstrates fitting a Poisson fixed effects model using the `femlm` function from the `FENmlm` package. Be aware of potential convergence warnings and check the output for reliability.

```r
rm(list=ls())
data <- read.csv(file="csv/fe1.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 | i + j
mod <- FENmlm::femlm(formula, data, family="poisson")
```

--------------------------------

### Inspecting Directions of Recession with ppmlhdfe

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_examples.md

Use `ppmlhdfe` with `tagsep(sep)` and `zvar(z)` options to inspect the directions of recession when dealing with separated observations. This command helps in recovering estimates when point estimates do not exist.

```stata
ppmlhdfe y x c, absorb(i) vce(cluster i) tagsep(sep) zvar(z) r2
```

--------------------------------

### Run FENmlm GLM

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.nb.html

Executes a Poisson GLM with fixed effects using the FENmlm package.

```R
rm(list=ls())
data <- read.csv(file="csv/fe1.csv", header=TRUE, sep=",")
formula <- y ~ x1 + x2 | i + j
mod <- FENmlm::femlm(formula, data, family="poisson")
summary(mod)
```

```R
rm(list=ls())
data <- read.csv(file="csv/fe2.csv", header=TRUE, sep=",")
formula <- y ~ x1 | i + j
try(mod <- FENmlm::femlm(formula, data, family="poisson"))
#summary(mod)
#diagnostic(mod)
```

--------------------------------

### Run Poisson regression with fixed effects in R using alpaca

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_benchmarks.md

Uses the alpaca package to estimate Poisson models with fixed effects. Note that this may fail to converge or handle separation correctly depending on the dataset.

```R
> # install.packages("alpaca")
> data <- read.csv(file="fe1.csv", header=TRUE, sep=",")
> library(alpaca)
> formula <- y ~ x1 + x2 | i + j
> mod <- feglm(formula, data, family=poisson())
> summary(mod)
```

```R
> data <- read.csv(file="fe2.csv", header=TRUE, sep=",")
> library(alpaca)
> formula <- y ~ x1 | i + j
> mod <- feglm(formula, data, family=poisson())
```

```R
> data <- read.csv(file="fe3.csv", header=TRUE, sep=",")
> library(alpaca)
> formula <- y ~ x1 | i + j
> mod <- feglm(formula, data, family=poisson())
> summary(mod)
```

--------------------------------

### Run Poisson regression in Matlab

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_benchmarks.md

Uses glmfit or fitglm to perform Poisson regression.

```Matlab
data = readtable('example1.csv');
X = [data.x1 data.x2 data.x3 data.x4];
Y = data.y;
[b,dev,stats] = glmfit(X,Y , 'poisson');
b'
dev
```

```Matlab
data = readtable('example2.csv');
X = [data.x1 data.x2 data.x3 data.x4];
Y = data.y;
model = fitglm(X, Y, 'linear', 'Distribution', 'Poisson');
```

--------------------------------

### Customize Regression Output Display

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/undocumented.md

Apply standard Stata estimation options to control confidence levels, formatting, and displayed coefficients in ppmlhdfe output.

```stata
sysuse auto, clear
gen w = weight
ppmlhdfe price weight w i.rep, level(90) cformat(%6.3fc) noci baselevels noomitted vsquish
```

--------------------------------

### Accessing Stored Results from ppmlhdfe

Source: https://context7.com/sergiocorreia/ppmlhdfe/llms.txt

Explains how to access comprehensive results stored in `e()` after running `ppmlhdfe` for postestimation analysis. This includes scalars, coefficient vectors, variance matrices, and macros.

```stata
sysuse auto, clear
ppmlhdfe price weight, absorb(turn)

* Key scalars
display "Observations: " e(N)
display "Pseudo R2: " e(r2_p)
display "Log-likelihood: " e(ll)
display "Deviance: " 2 * (e(ll_0) - e(ll))
display "Singletons dropped: " e(num_singletons)
display "Separated obs dropped: " e(num_separated)
display "Iterations: " e(ic)

* Coefficient vector and variance matrix
matrix list e(b)
matrix list e(V)

* Fixed effects degrees of freedom table
matrix list e(dof_table)

* Key macros
display "Dependent variable: " e(depvar)
display "Absorbed variables: " e(absvars)
display "VCE type: " e(vce)
display "Separation methods: " e(separation)
```

--------------------------------

### Reproduce Separation Issue in Stata

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/nonexistence_primer.md

Creates a synthetic dataset where separation occurs and attempts to run a standard Poisson regression.

```stata
* Create data
clear
set obs 6
gen y = max(0, _n-3)
gen x1 = 2*(_n==1) - (_n==2) + cond(_n>4, _n, 0)
gen x2 = 2 * (_n==2) - 2 * cond(_n>4, _n, 0) - (_n==1)
* Try to run -poisson-
poisson y x1 x2
```

--------------------------------

### Basic PPML Regression Comparison

Source: https://context7.com/sergiocorreia/ppmlhdfe/llms.txt

Compares the results of a standard Stata `poisson` regression with robust standard errors to an equivalent `ppmlhdfe` command without fixed effects, using sample airline data.

```stata
* Load sample data
use http://www.stata-press.com/data/r14/airline, clear

* Compare with standard poisson
poisson injuries XYZowned, vce(robust)
/*
------------------------------------------------------------------------------
     injuries | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
--------------+----------------------------------------------------------------
     XYZowned |   .6840413   .3895877     1.76   0.079    -.0795366    1.447619
        _cons |   3.688879   .0576012    64.04   0.000     3.575983    3.801776
------------------------------------------------------------------------------
*/

* Equivalent ppmlhdfe command
ppmlhdfe injuries XYZowned
/*
HDFE PPML regression                              No. of obs      =         21
                                                  Wald chi2(1)    =       3.08
Deviance             =  37.66127821               Prob > chi2     =     0.0791
Log pseudolikelihood =   -115.5043                Pseudo R2       =     0.0162
------------------------------------------------------------------------------
     injuries | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
--------------+----------------------------------------------------------------
     XYZowned |   .6840413   .3895877     1.76   0.079    -.0795366    1.447619
        _cons |   3.688879   .0576012    64.04   0.000     3.575982    3.801776
------------------------------------------------------------------------------
*/
```

--------------------------------

### Bootstrap Table Styling and MathJax Injection

Source: https://github.com/sergiocorreia/ppmlhdfe/blob/master/guides/r/glm_benchmarks.html

JavaScript functions to apply Bootstrap classes to Pandoc tables and dynamically load MathJax for rendering.

```JavaScript
function bootstrapStylePandocTables() { $('tr.header').parent('thead').parent('table').addClass('table table-condensed'); } $(document).ready(function () { bootstrapStylePandocTables(); }); (function () { var script = document.createElement("script"); script.type = "text/javascript"; script.src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"; document.getElementsByTagName("head")[0].appendChild(script); })();
```

--------------------------------

### Detecting Separation in Poisson Models

Source: https://context7.com/sergiocorreia/ppmlhdfe/llms.txt

Demonstrates how ppmlhdfe automatically identifies and drops separated observations that cause convergence issues in standard poisson models.

```stata
* Create data with separation
clear
set obs 6
gen y = max(_n - 3, 0)  // y = (0, 0, 0, 1, 2, 3)
gen x = _n < 3          // x = (1, 1, 0, 0, 0, 0)

* Standard poisson converges to incorrect estimates
poisson y x
/*
         y | Coefficient
-----------+------------
         x |   -18.8...   <- Arbitrary negative value, should be -infinity
     _cons |    .4054...
*/

* ppmlhdfe correctly detects separation
ppmlhdfe y x
/*
(dropped 2 separated observations identified by the FE method)
note: 1 variable omitted because of collinearity: x

         y | Coefficient
-----------+------------
         x |  (omitted)   <- Correctly omitted due to separation
     _cons |   .4054651
*/
```