### Example with Weights

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/reg.md

Demonstrates how to apply observation weights to a regression model.

```APIDOC
## With Weights

```julia
m = reg(df, @formula(Sales ~ NDI + fe(State)), weights = :Pop)
```
```

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

### FixedEffectModel Display Output Example

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/fixedeffectmodel.md

This is an example of the formatted summary output when a FixedEffectModel object is displayed. It includes key statistics like the number of observations, model convergence status, R-squared values, F-statistics, and a coefficient table.

```text
FixedEffectModel
===============================================================
Number of obs:                 1380   Converged:                     true
dof (model):                      1   dof (residuals):                 45
R²:                           0.803   R² adjusted:                  0.798
F-statistic:                13.3382   P-value:                      0.001
R² within:                    0.139   Iterations:                       5
===============================================================
      Estimate  Std. Error    t-stat  Pr(>|t|)   Lower 95%    Upper 95%
───────────────────────────────────────────────────────────────────────
NDI  -0.00526  0.00144  -3.652    0.0007  -0.0081  -0.0024
===============================================================

```

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

### Install FixedEffectModels.jl

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/API-OVERVIEW.md

Use the Julia package manager to add the FixedEffectModels package.

```julia
]
add FixedEffectModels
```

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

### Example with Fixed Effects

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/reg.md

Illustrates how to include fixed effects in a regression model.

```APIDOC
## With Fixed Effects

```julia
m = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)))
coeftable(m)
```
```

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

### Example Usage of fe() and @formula

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/types.md

Demonstrates how to create a `FixedEffectTerm` using `fe(:state)` and how to incorporate it into a formula using the `@formula` macro.

```julia
fe(:state)  # FixedEffectTerm(:state)
@formula(y ~ x + fe(state))  # fe(:state) inside formula
```

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

### Simple OLS Example

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/reg.md

A basic example demonstrating Ordinary Least Squares (OLS) regression without fixed effects.

```APIDOC
## Simple OLS

```julia
using RDatasets, FixedEffectModels
df = dataset("plm", "Cigar")
m = reg(df, @formula(Sales ~ NDI))
coef(m)
```
```

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

### Examples of Covariance Estimator Usage

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/types.md

Illustrates how to specify different covariance estimators (simple, robust, clustered) as the third argument to the `reg()` function.

```julia
reg(df, formula, Vcov.simple())
reg(df, formula, Vcov.robust())
reg(df, formula, Vcov.cluster(:state))
reg(df, formula, Vcov.cluster(:state, :year))
```

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

### Example with Clustered Standard Errors

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/reg.md

Shows how to compute clustered standard errors for robust inference.

```APIDOC
## With Clustered Standard Errors

```julia
using Vcov
m = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), Vcov.cluster(:State))
stderror(m)
```
```

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

### Example with Instrumental Variables

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/reg.md

Illustrates the usage of instrumental variables for Two-Stage Least Squares (2SLS) estimation.

```APIDOC
## With Instrumental Variables

```julia
m = reg(df, @formula(Sales ~ (Price ~ PriceMin)), method = :cpu)
```
```

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

### Formula Syntax for Fixed Effects and IV

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/API-OVERVIEW.md

Examples of using the formula syntax for basic regression, fixed effects, instrumental variables, and interactions.

```julia
# Basic: y ~ x1 + x2
@formula(y ~ x1 + x2)

# With fixed effects: y ~ x1 + fe(id) + fe(time)
@formula(y ~ x1 + fe(id) + fe(time))

# With instruments: y ~ x1 + (endogenous ~ instrument)
@formula(y ~ x1 + (price ~ supplycost))

# FE interactions
@formula(y ~ x1 + fe(state) & fe(year))  # Full state-year cells
@formula(y ~ x1 + fe(state) & year)      # State-specific slopes
@formula(y ~ x1 + fe(state) * year)      # Both effects and slopes

# Programmatic construction
term(:y) ~ term(:x1) + fe(:state)
```

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

### Sequential Model Building and R-squared Comparison

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/ADVANCED-USAGE.md

Build models sequentially, starting with a simple OLS and adding controls. This snippet demonstrates how to extract and compare R-squared values between models.

```julia
# Model 1: Simple OLS
m1 = reg(df, @formula(Sales ~ NDI), Vcov.robust())
r2_1 = r2(m1)

# Model 2: Add controls
m2 = reg(df, @formula(Sales ~ NDI + Price + Population), Vcov.robust())
r2_2 = r2(m2)
```

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

### Canonical Difference-in-Differences Setup

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/fixed-effect-terms.md

Use this snippet for standard difference-in-differences regressions with state and year fixed effects. Requires RDatasets and FixedEffectModels packages.

```julia
using RDatasets, FixedEffectModels

df = dataset("plm", "Cigar")
# Difference-in-differences: state and year fixed effects
m = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)))

```

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

### Plot Coefficients with StatsPlots.jl

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Example of plotting the coefficients and their standard errors of a fitted fixed effects model using StatsPlots.jl.

```julia
using StatsPlots
plot(coef(m), yerror = stderror(m))
```

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

### Basic OLS with Fixed Effects

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/API-OVERVIEW.md

Estimate an OLS model with fixed effects for 'State' and 'Year' using data from the 'Cigar' dataset. Includes examples of inspecting results like coefficients, standard errors, R-squared, and coefficient tables.

```julia
using RDatasets, FixedEffectModels

df = dataset("plm", "Cigar")

# OLS with fixed effects
m = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)))

# Inspect results
coef(m)       # Coefficients
stderror(m)   # Standard errors
r2(m)         # R-squared
coeftable(m)  # Formatted table
```

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

### Specifying Contrasts for Dummy Variables

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/README.md

Illustrates how to use the `contrasts` option to specify that a column should be treated as a set of dummy variables, with an optional base level.

```julia
reg(df, @formula(Sales ~ Price + Year); contrasts = Dict(:Year => DummyCoding()))
```

```julia
reg(df, @formula(Sales ~ Price + Year); contrasts = Dict(:Year => DummyCoding(base = 80)))
```

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

### GPU Acceleration

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/reg.md

Demonstrates how to leverage GPU acceleration for faster model fitting.

```APIDOC
## GPU Acceleration

```julia
using CUDA
m = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), method = :CUDA)
```
```

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

### Run a Simple Regression

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Use the `reg` function to perform a basic linear regression. The `coeftable` function displays the results.

```julia
using FixedEffectModels, RDatasets
df = dataset("plm", "Cigar")

m = reg(df, @formula(Sales ~ NDI))
coeftable(m)
```

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

### Timing Model Estimation with @time and @benchmark

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/ADVANCED-USAGE.md

Time model estimation using the `@time` macro for a quick estimate and the `@benchmark` macro for more accurate benchmarking of the `reg` function.

```julia
using BenchmarkTools

m = @time reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)))

# For accurate benchmarking
@benchmark reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)))
```

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

### Performance Tuning: CPU with Increased Tolerance and Reduced Precision

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/configuration.md

This snippet demonstrates performance tuning on the CPU by increasing the tolerance (`tol`) and reducing precision (`double_precision = false`). This can speed up computations at the cost of some accuracy.

```julia
# Or increase tolerance and reduce precision on CPU
m = reg(df, formula; tol = 1e-5, double_precision = false)
```

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

### Performance Tuning: Using GPU for Large Samples

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/configuration.md

This snippet shows how to leverage the GPU for performance gains when dealing with large samples or many fixed effects. It requires the CUDA package and sets `method` to `:CUDA` and `double_precision` to `false`.

```julia
# Use GPU if available
using CUDA
m = reg(df, formula; method = :CUDA, double_precision = false)
```

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

### Explicit Default Configuration for Fixed Effect Model

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/configuration.md

This snippet demonstrates the equivalent of the default workflow by explicitly listing all configuration arguments for the `reg` function. This is useful for understanding the default behavior.

```julia
m = reg(df, @formula(y ~ x1 + x2 + fe(state) + fe(year)), 
    Vcov.robust();
    contrasts = Dict(),
    weights = nothing,
    save = :none,
    method = :cpu,
    double_precision = true,
    tol = 1e-6,
    maxiter = 10000,
    drop_singletons = true,
    progress_bar = true,
    subset = nothing,
    first_stage = true)
```

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

### Get Predictions from a Model

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Generate predictions using the `predict` function. For models with fixed effects, ensure `save = :fe` was used during model fitting.

```julia
# Requires save = :fe if model has fixed effects
m = reg(df, @formula(Sales ~ NDI + fe(State)), save = :fe)

# Predict on original data
yhat = predict(m, df)

# Predict on subset
yhat_subset = predict(m, df[1:10, :])
```

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

### Get Fixed Effects Degrees of Freedom

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/fixedeffectmodel.md

Retrieves the total degrees of freedom absorbed by all fixed effects in the model. This is useful for understanding the model's complexity.

```julia
m = reg(df, @formula(y ~ x + fe(state) + fe(year)))
dof_fes(m)      # Returns number of state levels + year levels (adjusted for nesting)
dof(m)          # Returns number of non-FE regressors
```

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

### Displaying Model Summary and Coefficient Table

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

View the model summary directly in the REPL or generate a formatted coefficient table for detailed analysis.

```julia
# Display model summary
m  # Just type m at REPL
```

```julia
# Pretty coefficient table
coeftable(m)
```

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

### Checking Convergence and Diagnostics

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/ADVANCED-USAGE.md

Shows how to enable a progress bar during model fitting and access convergence information such as the `converged` status, number of `iterations`, and `r2_within` from the fitted model object.

```julia
m = reg(df, @formula(y ~ x + fe(id) + fe(year)));
        progress_bar = true)

# Check convergence status
m.converged     # true/false
m.iterations    # Number of iterations
m.r2_within     # Within-FE R²

```

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

### GPU Acceleration for Residualization

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/partial_out.md

Utilizes NVIDIA GPU acceleration for demeaning calculations, which can significantly speed up residualization on large datasets. Requires CUDA.jl to be installed and an NVIDIA GPU to be available.

```julia
using CUDA

# Use NVIDIA GPU for faster demeaning on large problems
result_df, _, _, _, _ = partial_out(
    df,
    @formula(y ~ x + fe(id_1) + fe(id_2)),
    method = :CUDA
)
```

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

### Programmatic Formula Construction

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/reg.md

Demonstrates how to dynamically construct formulas using the `term` function from StatsModels.jl.

```APIDOC
## Programmatic Formula Construction

For dynamic formula building:

```julia
using StatsModels: term
reg(df, term(:y) ~ term(:x1) + term(:x2) + fe(:state) + fe(:state) & term(:year))
```
```

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

### Fit Model with DataFrames.jl

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Example of fitting a fixed effects model using data from a DataFrame.jl DataFrame, including adding a new column with random data and an ID for fixed effects.

```julia
using DataFrames
df |> @transform(@:y = @:x .+ rand(nrow(@)), @:id = 1:(nrow(@)))
m = reg(df, @formula(y ~ x + fe(id)))
```

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

### Debugging Convergence Issues

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/configuration.md

This snippet shows how to debug convergence issues by tightening the tolerance (`tol`), increasing the maximum number of iterations (`maxiter`), enabling the progress bar, and disabling singleton dropping.

```julia
# Tighten tolerance and show progress
m = reg(df, formula; tol = 1e-8, maxiter = 50000, 
        progress_bar = true, drop_singletons = false)
```

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

### Selecting Computation Method

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/reg.md

Choose the computation method for estimation, including CPU (multi-threaded) or GPU acceleration (CUDA or Metal). GPU acceleration requires the respective packages to be loaded.

```julia
using FixedEffectModels

df = DataFrame(y = rand(100), x = rand(100), id = rand(1:10, 100))

# Use CPU (default)
reg(df, @formula(y ~ x + fe(id)), method=:cpu)

# Use CUDA (requires using CUDA first)
# reg(df, @formula(y ~ x + fe(id)), method=:CUDA)

# Use Metal (requires using Metal first)
# reg(df, @formula(y ~ x + fe(id)), method=:Metal)

```

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

### Get Fixed Effect Estimates

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Extract fixed effect estimates using the `fe` function. Ensure `save = :fe` was used during model fitting. Use `keepkeys = true` to include original grouping variables.

```julia
# Must save FE
m = reg(df, @formula(Sales ~ NDI + fe(State)), save = :fe)

# Get FE estimates
fe_est = fe(m)

# Include original grouping variables
fe_est_full = fe(m, keepkeys = true)
```

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

### Programmatic Formula Construction

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/README.md

Shows how to construct a regression formula programmatically using the `term` function for variable names and `fe` for fixed effects.

```julia
reg(df, term(:Sales) ~ term(:NDI) + fe(:State) + fe(:Year))
```

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

### Get Fixed Effect Estimates from FixedEffectModel

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/fixedeffectmodel.md

Returns estimated fixed effect coefficients as a DataFrame. This is only available if the model was fitted with fixed effects saved. The `keepkeys` argument controls whether original grouping variables are included in the output.

```julia
fe_estimates = fe(m; keepkeys=false)
```

```julia
m = reg(df, @formula(y ~ x + fe(state) + fe(year)), save = :fe)
fe_df = fe(m)  # DataFrame with state and year effects
fe_df_full = fe(m, keepkeys = true)  # Include original state and year columns
```

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

### Compare Models with RegressionTables.jl

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Demonstrates how to fit two different fixed effects models and then compare them side-by-side using RegressionTables.jl.

```julia
using RegressionTables
m1 = reg(df, @formula(y ~ x))
m2 = reg(df, @formula(y ~ x + fe(id)))
regtable(m1, m2)
```

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

### Load Data with CSV.jl and Fit Model

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Shows how to read data from a CSV file into a DataFrame using CSV.jl and then fit a fixed effects model.

```julia
using CSV
df = CSV.read("data.csv", DataFrame)
m = reg(df, @formula(Sales ~ NDI + fe(State)))
```

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

### Source Code Organization

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/API-OVERVIEW.md

Illustrates the directory structure and main files within the FixedEffectModels.jl package.

```julia
src/
├── FixedEffectModels.jl     # Main module, exports
├── FixedEffectModel.jl      # Model type, methods, display
├── fit.jl                   # reg() function (main estimation)
├── partial_out.jl           # partial_out() residualization
└── utils/
    ├── formula.jl           # fe(), formula parsing
    ├── fixedeffects.jl      # FE utilities (drop_singletons, etc.)
    ├── basecol.jl           # Matrix inversion, collinearity handling
    └── tss.jl               # Statistics (F-stat, TSS)

```

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

### Creating Formulas Programmatically and Handling Contrasts

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Construct formulas dynamically using the `term` and `fe` functions. Define custom contrasts for categorical variables using a dictionary.

```julia
# Create formula programmatically
term(:y) ~ term(:x1) + fe(:state)
```

```julia
# Check if a formula contains IV
has_iv(formula)
```

```julia
# Check if a formula contains FE
has_fe(formula)
```

```julia
# Create contrasts
using StatsModels: DummyCoding, EffectsCoding
contrasts = Dict(:var => DummyCoding(base=1))
```

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

### Select Computational Method for Demeaning

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/configuration.md

Choose the device or algorithm for demeaning operations. Defaults to multi-threaded CPU.

```julia
reg(df, formula; method = :cpu)        # Multi-threaded CPU (default)
reg(df, formula; method = :CUDA)       # NVIDIA GPU
reg(df, formula; method = :Metal)      # Apple GPU
```

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

### Inspecting Model Results

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/README.md

Lists common methods from `StatsAPI` that can be used to inspect the results of an estimated `FixedEffectModel`, such as retrieving coefficients, variance-covariance matrix, and confidence intervals.

```julia
# Coefficients
coef(m::FixedEffectModel)
```

```julia
vcov(m::FixedEffectModel)
```

```julia
confint(m::FixedEffectModel)
```

```julia
coefnames(m::FixedEffectModel)
```

```julia
responsename(m::FixedEffectModel)
```

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

### Handling Singleton Fixed Effects Groups

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/ADVANCED-USAGE.md

Illustrates the default behavior of dropping fixed effect groups with only one observation ('singletons'). Shows how to explicitly control this behavior with the `drop_singletons` option.

```julia
df = DataFrame(
    y = [1.0, 2.0, 3.0],
    x = [1.0, 2.0, 3.0],
    id = [1, 2, 3]  # All singletons
)

# Singletons dropped by default
m = reg(df, @formula(y ~ x + fe(id)); drop_singletons = true)
nobs(m)  # 0 observations (all were singletons)

# Keep singletons (not recommended)
m_unsafe = reg(df, @formula(y ~ x + fe(id)); drop_singletons = false)
```

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

### Configure Convergence and Tolerance

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/configuration.md

Adjust parameters for the iterative demeaning algorithm, including convergence tolerance, maximum iterations, and handling of singleton observations.

```julia
reg(df, formula; 
    tol = 1e-6,           # Convergence tolerance (default)
    maxiter = 10000,      # Maximum iterations (default)
    drop_singletons = true # Drop singleton observations (default)
)
```

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

### Common One-Liners for Fixed Effects Tasks

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Execute common fixed effects modeling tasks efficiently with these one-liner commands, including generating summaries, checking robustness, comparing models, and extracting results.

```julia
# Quick summary
display(coeftable(reg(df, @formula(y ~ x))))
```

```julia
# Check robustness to clustering
reg(df, @formula(y ~ x + fe(id)), Vcov.cluster(:id)).coef
```

```julia
# Compare models
[r2(reg(df, @formula(y ~ x))), r2(reg(df, @formula(y ~ x + fe(id))))]
```

```julia
# Extract residuals quickly
residuals(reg(df, @formula(y ~ x), save = :residuals))
```

```julia
# Plot FE estimates
fe(reg(df, @formula(y ~ x + fe(id)), save = :fe))
```

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

### Configure partial_out() Function

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/configuration.md

Set various options for the partial_out function, which shares many parameters with reg() but has specific defaults and features.

```julia
partial_out(df, formula;
    weights = nothing,
    add_mean = false,
    maxiter = 10000,
    contrasts = Dict(),
    method = :cpu,
    double_precision = true,
    tol = 1e-6,
    align = true,
    drop_singletons = false)

```

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

### Predict with DataFrame Input

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/errors.md

Shows how to correctly use the `predict` function by providing a DataFrame or a similar table-like structure for new data. Using non-table inputs will result in an error.

```julia
yhat = predict(m, df)
yhat = predict(m, DataFrame(x = [1, 2, 3]))
```

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

### Partial Out Function Usage

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/partial_out.md

Demonstrates the basic usage of the partial_out function with different fixed effect specifications, including additive, interactive, and group-specific slopes.

```APIDOC
## `partial_out(df, @formula(y ~ fe(state) + fe(year)))`

### Description
Calculates residuals for a linear model controlling for additive fixed effects for 'state' and 'year'.

### Method
`partial_out`

### Endpoint
N/A (Julia function)

### Parameters
- **df**: DataFrame containing the data.
- **formula**: Formula specifying the dependent variable, independent variables, and fixed effects using `@formula` and `fe()` syntax.

## `partial_out(df, @formula(y ~ fe(state) & fe(year)))`

### Description
Calculates residuals for a linear model controlling for interactive fixed effects between 'state' and 'year'.

### Method
`partial_out`

### Endpoint
N/A (Julia function)

### Parameters
- **df**: DataFrame containing the data.
- **formula**: Formula specifying the dependent variable, independent variables, and interactive fixed effects using `@formula` and `fe()` syntax.

## `partial_out(df, @formula(y ~ fe(state) & time))`

### Description
Calculates residuals for a linear model controlling for group-specific slopes where 'state' interacts with 'time'.

### Method
`partial_out`

### Endpoint
N/A (Julia function)

### Parameters
- **df**: DataFrame containing the data.
- **formula**: Formula specifying the dependent variable, independent variables, and group-specific slopes using `@formula` and `fe()` syntax.

## `partial_out(df, @formula(y ~ fe(state) * time))`

### Description
Calculates residuals for a linear model controlling for both additive fixed effects ('state') and interactive effects ('state' * 'time').

### Method
`partial_out`

### Endpoint
N/A (Julia function)

### Parameters
- **df**: DataFrame containing the data.
- **formula**: Formula specifying the dependent variable, independent variables, and combined fixed effects using `@formula` and `fe()` syntax.
```

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

### Generating Publication-Quality Regression Tables

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/ADVANCED-USAGE.md

Shows how to use the RegressionTables.jl package to create a publication-quality HTML table summarizing multiple regression models, including those with fixed effects and robust standard errors.

```julia
using RegressionTables

m1 = reg(df, @formula(Sales ~ NDI), Vcov.robust())
m2 = reg(df, @formula(Sales ~ NDI + fe(State)), Vcov.robust())
m3 = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), Vcov.robust())

# Create publication-quality table
regtable(m1, m2, m3;
         renderSettings = htmlOutput("output.html"))
```

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

### Adjust Model Estimation for Very Small Samples

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/ADVANCED-USAGE.md

For very small samples relative to the number of regressors, convergence issues may arise. Consider relaxing the tolerance and keeping all data by setting `drop_singletons = false`.

```julia
# Few observations, many fixed effects
if nrow(df) < 100
    # May encounter convergence issues
    m = reg(df, @formula(y ~ x + fe(id));
            tol = 1e-5,      # Relax tolerance
            drop_singletons = false)  # Keep data
end
```

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

### Monitor Sample Attrition with drop_singletons

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/ADVANCED-USAGE.md

Monitor sample attrition by comparing the original sample size with the number of observations used in the model. `drop_singletons = true` removes observations that are singletons within a fixed effect group.

```julia
m = reg(df, @formula(Sales ~ NDI + fe(State)));

# Check which observations were used
esample = m.esample
nrow(df)        # Original sample size
nobs(m)         # Used sample size
nrow(df) - nobs(m)  # Observations dropped
```

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

### GPU Acceleration for Large Fixed Effects Models

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/ADVANCED-USAGE.md

Compares the performance of CPU and GPU estimation for large datasets using FixedEffectModels.jl. Ensures results are numerically consistent between methods.

```julia
using CUDA, FixedEffectModels

# Large dataset
N = 10_000_000
K = 100
df_large = DataFrame(
    y = randn(N),
    x = randn(N, K)...,  # K regressors
    id1 = rand(1:100_000, N),
    id2 = rand(1:10_000, N)
)

# CPU estimation (slow for very large problems)
@time m_cpu = reg(df_large, @formula(y ~ x1 + x2 + ... + fe(id1) + fe(id2)), 
                  method = :cpu, double_precision = false)

# GPU estimation (fast)
@time m_gpu = reg(df_large, @formula(y ~ x1 + x2 + ... + fe(id1) + fe(id2)), 
                  method = :CUDA, double_precision = false)

# Verify same results (up to numerical precision)
@assert maximum(abs.(coef(m_cpu) .- coef(m_gpu))) < 1e-4
```

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

### Accessing Model Coefficients and Statistics

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Retrieve various model statistics such as coefficients, standard errors, confidence intervals, R-squared, and log-likelihood using dedicated functions.

```julia
coef(m)
```

```julia
stderror(m)
```

```julia
confint(m, level=0.95)
```

```julia
coefnames(m)
```

```julia
responsename(m)
```

```julia
vcov(m)
```

```julia
r2(m)
```

```julia
adjr2(m)
```

```julia
nobs(m)
```

```julia
dof(m)
```

```julia
dof_residual(m)
```

```julia
dof_fes(m)
```

```julia
rss(m)
```

```julia
nulldeviance(m)
```

```julia
mss(m)
```

```julia
loglikelihood(m)
```

```julia
formula(m)
```

```julia
coeftable(m)
```

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

### Access Model Properties

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Demonstrates how to access key properties of a fitted fixed effects model, such as the number of observations, degrees of freedom, and R-squared.

```julia
@show nobs(m), dof(m), dof_residual(m), r2(m)
```

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

### Use GPU Acceleration for Regressions

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Leverage GPU acceleration for faster computations by specifying the `method` argument. Supports CUDA for NVIDIA GPUs and Metal for Apple GPUs.

```julia
using CUDA

# NVIDIA GPU
m = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), 
        method = :CUDA)

# Apple GPU
using Metal
m = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), 
        method = :Metal)
```

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

### Estimate Model with Metal GPU Acceleration

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/README.md

Use this snippet to estimate a fixed effect model using Metal-enabled GPUs for potential speedups. Ensure Metal is functional and the `method` is set to `:Metal`. Supports only single precision.

```julia
using Metal, FixedEffectModels, RDatasets
@assert Metal.functional()
df = dataset("plm", "Cigar")
reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), method = :Metal)
```

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

### Formula Syntax for Fixed Effects

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/README.md

Demonstrates constructing formulas with dependent variables, exogenous variables, and high-dimensional fixed effects indicated by the `fe` function. Multiple fixed effects can be added with `+`, interacted with `&`, or combined with variables using `*`.

```julia
reg(df, @formula(Sales ~ Price + fe(State) + fe(Year)))
```

```julia
reg(df, @formula(Sales ~ NDI + fe(State) + fe(State)&Year))
```

```julia
reg(df, @formula(Sales ~ NDI + fe(State)&fe(Year)))
```

```julia
reg(df, @formula(Sales ~ (Price ~ Pimin)))
```

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

### Access Fixed Effects for Models with FE

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/errors.md

Demonstrates the correct way to access fixed effect estimates from a model that was fitted with the `save = :fe` option. This is only applicable to models that actually include fixed effects.

```julia
m = reg(df, @formula(y ~ x + fe(id)), save = :fe)
fe_estimates = fe(m)
```

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

### Visualization Workflow with Residualization

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/partial_out.md

Demonstrates a workflow for visualizing conditional relationships by first partialing out confounding variables (fixed effects) from both the dependent and independent variables, and then plotting the residualized relationship.

```julia
using Gadfly

# Partial out both axes with respect to confounders
result, _, _, _, _ = partial_out(
    df,
    @formula(wage + experience ~ fe(education) + fe(year)),
    add_mean = true
)

# Now plot residualized relationship
plot(result, x = :experience, y = :wage, Geom.point, Geom.smooth)
```

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

### Comparing Standard Error Methods

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/ADVANCED-USAGE.md

Estimate a model and then extract coefficients and standard errors using different methods (simple and robust) without re-estimating the model.

```julia
m = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), Vcov.simple())

# Extract results with different SE estimators
coef_simple = coef(m)

# Re-estimate with robust SEs
m_robust = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), Vcov.robust())
coef_robust = coef(m_robust)  # Same coefficients, different SEs
stderror(m_robust)
```

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

### Common Workflow: Visualization

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/partial_out.md

Illustrates a common workflow for visualizing the conditional relationship between two variables after controlling for confounders using `partial_out` and plotting.

```APIDOC
## `partial_out(df, @formula(Sales + NDI ~ fe(Year)), add_mean = true)`

### Description
Residualizes 'Sales' and 'NDI' by controlling for 'Year' fixed effects, and adds back the original means of these variables.

### Method
`partial_out`

### Endpoint
N/A (Julia function)

### Parameters
- **df**: DataFrame containing the data.
- **formula**: Formula specifying multiple dependent variables ('Sales', 'NDI') and fixed effects ('Year').
- **add_mean**: Boolean, defaults to `true`. If `true`, the original means of the dependent variables are added back to the residuals.

## Plotting Conditional Relationship

### Description
Plots the conditional relationship between 'NDI' (x-axis) and 'Sales' (y-axis) using a linear model fit, after controlling for 'Year' fixed effects.

### Method
`plot` (from Gadfly.jl)

### Endpoint
N/A (Julia function)

### Parameters
- **residualized**: DataFrame returned by `partial_out`.
- **x**: Symbol representing the variable for the x-axis (:NDI).
- **y**: Symbol representing the variable for the y-axis (:Sales).
- **Geom.point**: Specifies to draw points for each observation.
- **Geom.smooth(method=:lm)**: Specifies to draw a smoothed line representing the linear model fit.
```

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

### coeftable

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/fixedeffectmodel.md

Creates a formatted coefficient table for a FixedEffectModel, including standard errors, t-statistics, p-values, and confidence intervals. The output is a `CoefTable` object that provides a nicely formatted display when printed.

```APIDOC
## coeftable(m::FixedEffectModel; level=0.95)

### Description
Creates a formatted coefficient table with standard errors, t-statistics, p-values, and confidence intervals. Returns a `CoefTable` object that displays nicely when printed.

### Parameters
- **m** (`FixedEffectModel`) - The fitted FixedEffectModel object.
- **level** (`Float64`) - Optional. The confidence level for the intervals, defaults to 0.95.

### Returns
`CoefTable` with columns: Estimate, Std. Error, t-stat, Pr(>|t|), Lower XX%, Upper XX%

### Example
```julia
m = reg(df, @formula(y ~ x + fe(id)))
coeftable(m)  # Display nicely
ct = coeftable(m, level = 0.90)  # 90% confidence intervals
```
```

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

### Accessing Model Convergence Status and Iterations

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Check if the demeaning algorithm converged and retrieve the number of iterations performed during model estimation.

```julia
# Did demeaning converge?
m.converged
```

```julia
# Iterations performed
m.iterations
```

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

### Efficiently Handling Categorical Variables with Many Levels

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/ADVANCED-USAGE.md

Recommends using `fe()` for categorical variables with many levels instead of explicit dummy variables (`C()`) for significantly improved performance.

```julia
# Instead of using explicit dummy variables,
# use fe() which is much faster:

# Slow (creates explicit dummies)
# m_slow = reg(df, @formula(y ~ x + C(category)))

# Fast (absorbs as fixed effect)
m_fast = reg(df, @formula(y ~ x + fe(category)))
```

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

### Custom Contrasts

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/reg.md

Shows how to specify custom contrasts for categorical variables.

```APIDOC
## Custom Contrasts

```julia
using StatsModels: DummyCoding
df.YearCat = categorical(df.Year)
m = reg(df, @formula(Sales ~ YearCat), contrasts = Dict(:YearCat => DummyCoding(base = 1980)))
```
```

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

### IV Formula with No Endogenous Variables

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/errors.md

Shows how to correctly specify both endogenous and instrumental variables in an IV formula. Using IV syntax without specifying endogenous variables leads to an error.

```julia
# Incorrect
@formula(y ~ x + (~ instrument))

# Correct
@formula(y ~ x + (endogenous ~ instrument))
```

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

### Handle Empty Sample in partial_out

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/errors.md

Addresses the 'sample is empty' error in `partial_out()`, which occurs when no observations remain after filtering. This is analogous to the same error in `reg()`.

```julia
# This error occurs if all observations are filtered out, e.g., due to missing values.
```

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

### Correct Formula Structure

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/errors.md

Demonstrates the correct way to structure formulas, especially when dealing with instrumental variables. Avoid nesting formulas on the left-hand side.

```julia
# Incorrect
@formula((y ~ x) ~ z)

# Correct
@formula(y ~ (x ~ z))
```

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

### Performance Tuning: Very Large Datasets

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/configuration.md

This snippet is for optimizing performance on very large datasets by trading precision for speed. It uses the GPU (`method = :CUDA`), disables the progress bar, and reduces precision.

```julia
# Trade precision for speed
m = reg(df, formula; method = :CUDA, double_precision = false, 
        progress_bar = false)
```

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

### GPU Acceleration for Fixed Effects Models

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/API-OVERVIEW.md

Utilize GPU acceleration for fixed effects models using CUDA for NVIDIA GPUs or Metal for Apple GPUs, potentially offering significant speedups.

```julia
using CUDA

# NVIDIA GPU (2-5× speedup for large problems)
m = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), 
        method = :CUDA)

# Apple GPU
using Metal
m = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), 
        method = :Metal)
```

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

### Specifying Computation Method

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/README.md

Indicates the hardware acceleration method for computation, supporting CPU, CUDA (NVIDIA GPUs), or Metal (Apple GPUs).

```julia
method = :cpu
```

```julia
method = :CUDA
```

```julia
method = :Metal
```

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

### Handling Nested Fixed Effects

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/fixed-effect-terms.md

Demonstrates how to specify potentially nested fixed effects (e.g., worker within firm) in a formula. The package automatically handles nesting and adjusts degrees of freedom.

```julia
# Worker nested within Firm (each worker belongs to one firm)
m = reg(df, @formula(y ~ x + fe(worker_id) + fe(firm_id)))

```

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

### Fixed Effect Interactions

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/fixed-effect-terms.md

Demonstrates various ways to combine fixed effects within regression formulas using operators like `+`, `&`, and `*`.

```APIDOC
## Fixed Effect Interactions

### Additive Fixed Effects

Multiple fixed effects with `+` operator:

```julia
# Separate State and Year effects
m = reg(df, @formula(y ~ x + fe(State) + fe(Year)))
```

The model absorbs State effects and Year effects independently. Coefficients estimate effects that vary across states and years but in an additive way.

### Interaction of Fixed Effects (Full Cells)

Use `&` operator to interact two fixed effects:

```julia
# Full State-Year interaction (separate effect per state-year combination)
m = reg(df, @formula(y ~ x + fe(State) & fe(Year)))
```

Creates a separate fixed effect for each combination of State and Year that appears in the data.

### Fixed Effect with Continuous Interaction (Group-Specific Slopes)

Interact a fixed effect with a continuous variable using `&`:

```julia
# State-specific slopes on Year (continuous)
m = reg(df, @formula(y ~ x + fe(State) & Year))
```

Estimates a separate slope coefficient on the continuous variable `Year` for each `State`, while constraining the overall intercept effect of `Year` to zero (absorbed by fixed effect).

### Combined Additive and Interaction Terms

Use `*` operator to get both additive and interaction effects:

```julia
# Equivalent to: fe(State) + fe(State)&Year
m = reg(df, @formula(y ~ x + fe(State) * Year))
```

Estimates State effects plus State-specific slopes on Year.

### Multiple Fixed Effect Interactions

Interact multiple fixed effects:

```julia
# Three-way State-Region-Year interaction
m = reg(df, @formula(y ~ x + fe(State) & fe(Region) & fe(Year)))
```

Each combination of State, Region, and Year gets its own fixed effect coefficient.
```

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

### Define Verbose Formula

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/QUICK-REFERENCE.md

Shows how to define a verbose formula for a fixed effects model.

```julia
formula(m)
```

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

### Use reg() for IV Models

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/errors.md

Demonstrates that `partial_out()` does not support instrumental variables. For models involving instrumental variables, the `reg()` function should be used instead.

```julia
# partial_out(df, @formula(y ~ (x ~ z))) # This would cause an error
m = reg(df, @formula(y ~ (x ~ z)))
```

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

### Main Functions

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/API-OVERVIEW.md

The primary functions available in the FixedEffectModels.jl package for estimation and data manipulation.

```APIDOC
## `reg()`

### Description
Estimate regression models with fixed effects and/or instrumental variables.

### Purpose
This is the core function for estimating various types of regression models, including OLS, 2SLS, and models with fixed effects.

### Returns
- `FixedEffectModel`: An object representing the fitted regression model.

## `partial_out()`

### Description
Residualize variables by removing the effect of specified fixed effects.

### Purpose
Useful for data preprocessing, visualization, or preparing variables for further analysis.

### Returns
- `DataFrame`: A DataFrame containing the residualized variables.

## `fe()`

### Description
Marks variables as fixed effects within regression formulas.

### Purpose
Used in conjunction with the `@formula` macro to specify which categorical variables should be treated as fixed effects.

### Returns
- `FixedEffectTerm`: An object representing a fixed effect term in a formula.
```

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

### Apply Observation Weights in reg()

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/configuration.md

Use the `weights` keyword argument to specify a column containing observation weights. Weights must be positive; missing or non-positive weights are automatically dropped.

```julia
reg(df, formula; weights = :population)
```

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

### Default Workflow for Fixed Effect Model

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/configuration.md

This snippet shows a typical call to the `reg` function with default settings for estimating a fixed effect model. It includes the formula and a robust covariance estimator.

```julia
m = reg(df, @formula(y ~ x1 + x2 + fe(state) + fe(year)), 
         Vcov.robust())
```

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

### GPU Acceleration for Regression

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/reg.md

Estimates a regression model using GPU acceleration by specifying `method = :CUDA`. Requires the CUDA package.

```Julia
using CUDA
m = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), method = :CUDA)
```

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

### Manual Prediction on Subset and Out-of-Sample Data

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/ADVANCED-USAGE.md

Explains how to perform predictions using a fitted FixedEffectModels.jl model on a subset of the original data or on new, out-of-sample data. Notes limitations for out-of-sample fixed effect prediction.

```julia
m = reg(df, @formula(Sales ~ NDI + fe(State)), save = :fe)

# Subset prediction
subset = df.Year .== 2000
pred_subset = predict(m, df[subset, :])

# Out-of-sample prediction (only works for observed state values)
new_data = DataFrame(State = ["CA", "NY"], NDI = [50.0, 45.0])
try
    pred_oos = predict(m, new_data)
catch e
    println("Out-of-sample FE prediction not supported")
end
```

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

### Basic Formula Structure with Fixed Effects

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/api-reference/fixed-effect-terms.md

Illustrates the basic structure of a formula including dependent variable, regressors, and fixed effects defined with fe().

```julia
@formula(y ~ x1 + x2 + fe(group) + fe(time))

```

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

### Julia FixedEffectModels Benchmark

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/benchmark/benchmark.md

Generates a large dataset and benchmarks different regression models with fixed effects and clustered standard errors using FixedEffectModels.jl.

```julia
using DataFrames, CategoricalArrays, FixedEffectModels
N = 10_000_000
K = 100
id1 = rand(1:(N/K), N)
id2 = rand(1:K, N)
x1 =  randn(N)
x2 =  randn(N)
y= 3 .* x1 .+ 2 .* x2 .+ sin.(id1) .+ cos.(id2).^2 .+ randn(N)
df = DataFrame(id1 = categorical(id1), id2 = categorical(id2), x1 = x1, x2 = x2, y = y)
@time reg(df, @formula(y ~ x1 + x2))
# 0.338749 seconds (450 allocations: 691.441 MiB, 2.30% gc time)
@time reg(df, @formula(y ~ x1 + x2 + fe(id1)))
#  0.463058 seconds (1.00 k allocations: 929.129 MiB, 13.31% gc time)
@time reg(df, @formula(y ~ x1 + x2 + fe(id1) + fe(id2)))
#  1.006031 seconds (3.22 k allocations: 1.057 GiB, 1.68% gc time)
@time reg(df, @formula(y ~ x1 + x2), Vcov.cluster(:id1))
# 0.380562 seconds (580 allocations: 771.606 MiB, 3.07% gc time)
@time reg(df, @formula(y ~ x1 + x2), Vcov.cluster(:id1, :id2))
#0.765847 seconds (719 allocations: 1.128 GiB, 2.01% gc time)
```

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

### Re-estimate Model to Save Residuals

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/errors.md

Demonstrates how to re-estimate a fixed effect model to save residuals. This is necessary when residuals are needed for post-estimation analysis but were not initially saved.

```julia
m = reg(df, @formula(y ~ x + fe(id)))
# residuals(m) # This would cause an error
m = reg(df, @formula(y ~ x + fe(id)), save = :residuals)
e = residuals(m)
```

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

### Visualization with Partialing Out

Source: https://github.com/fixedeffects/fixedeffectmodels.jl/blob/main/_autodocs/API-OVERVIEW.md

Obtain residuals from a regression model after partialing out fixed effects, suitable for plotting the conditional relationship between variables. Requires importing Plots.

```julia
# Residualize both outcome and regressor
residuals, _, _, _, _ = partial_out(
    df,
    @formula(Sales + NDI ~ fe(State) + fe(Year)),
    add_mean = true  # Add back means for interpretability
)

# Plot conditional relationship
using Plots
scatter(residuals.NDI, residuals.Sales)
```