### Setup for GLM.jl Examples

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

This code block sets up the necessary packages and environment for running examples in the GLM.jl documentation. It includes using CategoricalArrays, DataFrames, Distributions, GLM, RDatasets, and Optim.

```julia
using CategoricalArrays, DataFrames, Distributions, GLM, RDatasets, Optim
```

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

### Generalized Linear Model Setup in R

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Example of setting up data and fitting a Poisson generalized linear model in R using glm().

```r
glm> ## Dobson (1990) Page 93: Randomized Controlled Trial : (slightly modified)
glm> counts <- c(18,17,15,20,10,21,25,13,13)

glm> outcome <- gl(3,1,9)

glm> treatment <- gl(3,3)

glm> print(d.AD <- data.frame(treatment, outcome, counts))
  treatment outcome counts
1         1       1     18
2         1       2     17
3         1       3     15
4         2       1     20
5         2       2     10
6         2       3     21
7         3       1     25
8         3       2     13
9         3       3     13

glm> glm.D93 <- glm(counts ~ outcome + treatment, family=poisson())

glm> anova(glm.D93)
Analysis of Deviance Table

Model: poisson, link: log


```

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

### Install GLM.jl Package

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Use Pkg.add to install the GLM package and its dependencies.

```julia
Pkg.add("GLM")
```

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

### GLM with Custom Starting Values and Convergence Control

Source: https://context7.com/juliastats/glm.jl/llms.txt

Fit a GLM model with custom starting coefficients and control convergence parameters. This allows for fine-tuning the optimization process.

```julia
glm(@formula(Y ~ X1 + X2), binary_data, Binomial();
    start=[0.0, 0.0, 0.0],  # starting coefficients
    maxiter=50,             # max iterations
    atol=1e-8,              # absolute tolerance
    rtol=1e-8,              # relative tolerance
    minstepfac=0.001)       # minimum step factor for line search
```

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

### Refitting Models with `fit!`

Source: https://context7.com/juliastats/glm.jl/llms.txt

Demonstrates how to refit an existing model using the `fit!` function. For `LinearModel`, this recomputes coefficients. For `GeneralizedLinearModel`, it allows specifying convergence parameters and starting values.

```julia
model = lm(@formula(Y ~ X), data)
fit!(model)

glm_model = GeneralizedLinearModel(...)
fit!(glm_model;
     maxiter=50,
     minstepfac=0.0001,
     atol=1e-8,
     rtol=1e-8,
     start=[0.0, 0.0])
```

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

### Simple Linear Model in R

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Example of fitting a simple linear model using lm() in R with the Formaldehyde dataset.

```r
> coef(summary(lm(optden ~ carb, Formaldehyde)))
               Estimate  Std. Error    t value     Pr(>|t|)
(Intercept) 0.005085714 0.007833679  0.6492115 5.515953e-01
carb        0.876285714 0.013534536 64.7444207 3.409192e-07
```

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

### Example GLM Debug Output

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

This is an example of the debug output generated when fitting a generalized linear model with debugging enabled. It shows iteration progress, deviance, and change in deviance.

```julia
┌ Debug: Iteration: 1, deviance: 5.129147109764238, diff.dev.:0.05057195315968688
└ @ GLM ~/.julia/dev/GLM/src/glmfit.jl:418
```

```julia
┌ Debug: Iteration: 2, deviance: 5.129141077001254, diff.dev.:6.032762984276019e-6
└ @ GLM ~/.julia/dev/GLM/src/glmfit.jl:418
```

```julia
┌ Debug: Iteration: 3, deviance: 5.129141077001143, diff.dev.:1.1102230246251565e-13
└ @ GLM ~/.julia/dev/GLM/src/glmfit.jl:418
```

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

### Probit Regression

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Example of performing probit regression on a small dataset. This is suitable for binary outcomes.

```julia
julia> data = DataFrame(X=[1,2,2], Y=[1,0,1])
3×2 DataFrame
 Row │ X      Y     
     │ Int64  Int64 
─────┼──────────────
   1 │     1      1
   2 │     2      0
   3 │     2      1

julia> probit = glm(@formula(Y ~ X), data, Binomial(), ProbitLink())
GeneralizedLinearModel

Y ~ 1 + X

Coefficients:
────────────────────────────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)  Lower 95%  Upper 95%
────────────────────────────────────────────────────────────────────────
(Intercept)   9.63839     293.909   0.03    0.9738   -566.414    585.69
X            -4.81919     146.957  -0.03    0.9738   -292.849    283.211
────────────────────────────────────────────────────────────────────────
```

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

### Simple Linear Model in Julia

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Corresponds to the R example, fitting a linear model using GLM.jl with the Formaldehyde dataset. Requires importing GLM and RDatasets.

```julia
julia> using GLM, RDatasets

julia> form = dataset("datasets", "Formaldehyde")
6×2 DataFrame
 Row │ Carb     OptDen  
     │ Float64  Float64 
─────┼──────────────────
   1 │     0.1    0.086
   2 │     0.3    0.269
   3 │     0.5    0.446
   4 │     0.6    0.538
   5 │     0.7    0.626
   6 │     0.9    0.782

julia> lm1 = fit(LinearModel, @formula(OptDen ~ Carb), form)
LinearModel

OptDen ~ 1 + Carb

Coefficients:
───────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      t  Pr(>|t|)   Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────────
(Intercept)  0.00508571  0.00783368   0.65    0.5516  -0.0166641  0.0268355
Carb         0.876286    0.0135345   64.74    <1e-06   0.838708   0.913864
───────────────────────────────────────────────────────────────────────────
```

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

### Complex Linear Model in R

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Example of fitting a multiple linear regression model in R using lm() with the LifeCycleSavings dataset.

```r
> coef(summary(lm(sr ~ pop15 + pop75 + dpi + ddpi, LifeCycleSavings)))
                 Estimate   Std. Error    t value     Pr(>|t|)
(Intercept) 28.5660865407 7.3545161062  3.8841558 0.0003338249
pop15       -0.4611931471 0.1446422248 -3.1885098 0.0026030189
pop75       -1.6914976767 1.0835989307 -1.5609998 0.1255297940
dpi         -0.0003369019 0.0009311072 -0.3618293 0.7191731554
ddpi         0.4096949279 0.1961971276  2.0881801 0.0424711387
```

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

### Complex Linear Model in Julia

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Corresponds to the R example, fitting a multiple linear regression model using GLM.jl with the LifeCycleSavings dataset. Requires importing RDatasets.

```julia
julia> LifeCycleSavings = dataset("datasets", "LifeCycleSavings")
50×6 DataFrame
 Row │ Country         SR       Pop15    Pop75    DPI      DDPI
     │ String15        Float64  Float64  Float64  Float64  Float64
─────┼─────────────────────────────────────────────────────────────
   1 │ Australia         11.43    29.35     2.87  2329.68     2.87
   2 │ Austria           12.07    23.32     4.41  1507.99     3.93
   3 │ Belgium           13.17    23.8      4.43  2108.47     3.82
   4 │ Bolivia            5.75    41.89     1.67   189.13     0.22
   5 │ Brazil            12.88    42.19     0.83   728.47     4.56
   6 │ Canada             8.79    31.72     2.85  2982.88     2.43
   7 │ Chile              0.6     39.74     1.34   662.86     2.67
   8 │ China             11.9     44.75     0.67   289.52     6.51
  ⋮  │       ⋮            ⋮        ⋮        ⋮        ⋮        ⋮
  44 │ United States      7.56    29.81     3.43  4001.89     2.45
  45 │ Venezuela          9.22    46.4      0.9    813.39     0.53
  46 │ Zambia            18.56    45.25     0.56   138.33     5.14
  47 │ Jamaica            7.72    41.12     1.73   380.47    10.23
  48 │ Uruguay            9.24    28.13     2.72   766.54     1.88
  49 │ Libya              8.89    43.69     2.07   123.58    16.71
  50 │ Malaysia           4.71    47.2      0.66   242.69     5.08
                                                    35 rows omitted

julia> fm2 = lm(@formula(SR ~ Pop15 + Pop75 + DPI + DDPI), LifeCycleSavings)
LinearModel

SR ~ 1 + Pop15 + Pop75 + DPI + DDPI

Coefficients:
─────────────────────────────────────────────────────────────────────────────────
                    Coef.   Std. Error      t  Pr(>|t|)    Lower 95%    Upper 95%
─────────────────────────────────────────────────────────────────────────────────
(Intercept)  28.5661       7.35452       3.88    0.0003  13.7533      43.3788
Pop15        -0.461193     0.144642     -3.19    0.0026  -0.752518    -0.169869
Pop75        -1.6915       1.0836       -1.56    0.1255  -3.87398      0.490983
DPI          -0.000336902  0.000931107  -0.36    0.7192  -0.00221225   0.00153844
DDPI          0.409695     0.196197      2.09    0.0425   0.0145336    0.804856
─────────────────────────────────────────────────────────────────────────────────
```

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

### Model Constructors

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/api.md

Demonstrates how to fit a LinearModel using the `fit` function and the shorthand `lm` function. Both methods produce identical results for linear models.

```APIDOC
## Model Constructors

The most general approach to fitting a model is with the `fit` function, as in

```jldoctest constructors
julia> using RDatasets

julia> df = RDatasets.dataset("mlmRev", "Oxboys");

julia> fit(LinearModel, hcat(ones(nrow(df)), df.Age), df.Height)
LinearModel

Coefficients:
─────────────────────────────────────────────────────────────────
        Coef.  Std. Error       t  Pr(>|t|)  Lower 95%  Upper 95%
─────────────────────────────────────────────────────────────────
x1  149.372      0.528565  282.60    <1e-99  148.33     150.413
x2    6.52102    0.816987    7.98    <1e-13    4.91136    8.13068
─────────────────────────────────────────────────────────────────
```

This model can also be fit as

```jldoctest constructors
julia> lm(hcat(ones(nrow(df)), df.Age), df.Height)
LinearModel

Coefficients:
─────────────────────────────────────────────────────────────────
        Coef.  Std. Error       t  Pr(>|t|)  Lower 95%  Upper 95%
─────────────────────────────────────────────────────────────────
x1  149.372      0.528565  282.60    <1e-99  148.33     150.413
x2    6.52102    0.816987    7.98    <1e-13    4.91136    8.13068
─────────────────────────────────────────────────────────────────
```

### Methods

- `fit(model_type, X, y)`: Fits a specified model type to the data.
- `lm(X, y)`: A shorthand for fitting a linear model.
```

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

### Fit Negative Binomial Model with Initial Shape Parameter

Source: https://context7.com/juliastats/glm.jl/llms.txt

Demonstrates how to fit a Negative Binomial model, optionally providing an initial estimate for the shape parameter θ.

```julia
theta = nb_model.rr.d.r

nb_fixed = glm(@formula(Days ~ Eth + Sex + Age + Lrn), quine,
               NegativeBinomial(2.0), LogLink())

nb_init = negbin(@formula(Days ~ Eth + Sex + Age + Lrn), quine, LogLink();
                 initialθ=1.5)

nb_precise = negbin(@formula(Days ~ Eth + Sex + Age + Lrn), quine, LogLink();
                    maxiter=100,
                    atol=1e-8,
                    rtol=1e-8)

nb_canonical = negbin(@formula(Days ~ Eth + Sex + Age + Lrn), quine)
```

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

### Model Selection using AIC

Source: https://context7.com/juliastats/glm.jl/llms.txt

Demonstrates how to select the best model from a list of candidate models by comparing their AIC values. The model with the minimum AIC is chosen as the preferred model.

```julia
models = [model1, model2, model3]
aics = aic.(models)
best_model = models[argmin(aics)]
```

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

### Fitting Linear Models with Weights

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Demonstrates how to fit a linear model using weighted observations.

```APIDOC
## POST /api/models/linear/weighted

### Description
Fits a linear model to the provided data, incorporating observation weights.

### Method
POST

### Endpoint
/api/models/linear/weighted

### Parameters
#### Request Body
- **formula** (string) - Required - The model formula (e.g., "y ~ x").
- **data** (object) - Required - The dataset containing the variables.
- **weights** (string) - Required - The name of the column in the data to be used as weights.

### Request Example
```json
{
  "formula": "y ~ x",
  "data": {
    "y": [1, 2, 3, 4, 5],
    "x": [10, 20, 30, 40, 50]
  },
  "weights": "weights"
}
```

### Response
#### Success Response (200)
- **model_type** (string) - The type of model fitted (e.g., "LinearModel").
- **coefficients** (object) - Coefficients of the fitted model.
- **fit_statistics** (object) - Statistics related to the model fit.

#### Response Example
```json
{
  "model_type": "LinearModel",
  "coefficients": {
    "(Intercept)": 0.51673,
    "x": -0.0478667
  },
  "fit_statistics": {
    "r2": 0.9399,
    "adj_r2": 0.935,
    "aic": 100.5,
    "bic": 102.1
  }
}
```
```

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

### Get Canonical Link for a Distribution

Source: https://context7.com/juliastats/glm.jl/llms.txt

Determine the canonical link function for a given distribution family. The canonical link is automatically used by default in GLM.jl.

```julia
canonicallink(Binomial())        # LogitLink()
canonicallink(Poisson())         # LogLink()
canonicallink(Normal())          # IdentityLink()
canonicallink(Gamma())           # InverseLink()
canonicallink(NegativeBinomial(1.0))  # NegativeBinomialLink(1.0)
```

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

### Fitting Linear Model with `fit`

Source: https://context7.com/juliastats/glm.jl/llms.txt

Demonstrates fitting a linear model using the generic `fit` function, which is equivalent to using the `lm` function. This method allows for explicit specification of the model type.

```julia
using GLM, DataFrames
data = DataFrame(X=[1,2,3,4,5], Y=[2,4,5,4,6])
model1 = fit(LinearModel, @formula(Y ~ X), data)
model2 = lm(@formula(Y ~ X), data)
```

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

### Fitting Models with Matrix Interface

Source: https://context7.com/juliastats/glm.jl/llms.txt

Illustrates fitting linear and generalized linear models using a matrix-based interface, where the design matrix (X) and response vector (y) are provided directly. This is an alternative to the formula interface.

```julia
X = hcat(ones(5), [1,2,3,4,5])
y = [2,4,5,4,6]
model_mat = fit(LinearModel, X, y)
glM_mat = fit(GeneralizedLinearModel, X, y, Normal(), IdentityLink())
```

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

### Linear Regression with QR Method

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Demonstrates linear regression using the QR method, which is generally more robust to multicollinearity than the Cholesky method.

```julia
julia> mdl = lm(@formula(rdintens ~ sales + sales^2), rdchem; method=:qr)
LinearModel

rdintens ~ 1 + sales + :(sales ^ 2)

Coefficients:
─────────────────────────────────────────────────────────────────────────────────
                    Coef.   Std. Error      t  Pr(>|t|)    Lower 95%    Upper 95%
─────────────────────────────────────────────────────────────────────────────────
(Intercept)   2.61251      0.429442      6.08    <1e-05   1.73421     3.49082
sales         0.000300571  0.000139295   2.16    0.0394   1.56805e-5  0.000585462
sales ^ 2    -6.94594e-9   3.72614e-9   -1.86    0.0725  -1.45667e-8  6.7487e-10
─────────────────────────────────────────────────────────────────────────────────
```

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

### Fitting Generalized Linear Model with `fit`

Source: https://context7.com/juliastats/glm.jl/llms.txt

Shows how to fit a generalized linear model using the generic `fit` function, which is equivalent to using the `glm` function. This requires specifying the model type, formula, data, distribution, and optionally a link function.

```julia
glm_model1 = fit(GeneralizedLinearModel, @formula(Y ~ X), data, Normal())
glm_model2 = glm(@formula(Y ~ X), data, Normal())
```

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

### Predict Response Values

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Generate predicted response values from a fitted model using the `predict` function. Provide `test_data` for new covariate values, or omit it to get fitted values.

```julia
test_data = DataFrame(X=[4])

round.(predict(mdl, test_data); digits=8)
```

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

### Compare Cholesky and QR Factorization for Ill-Conditioned Matrix

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Compares the results of fitting a linear model using Cholesky and QR factorization methods on an ill-conditioned design matrix. The QR method is generally preferred for better numerical stability and accuracy in such cases, as indicated by a higher log-likelihood.

```julia
X = [-0.4011512997627107 0.6368622664511552;
            -0.0808472925693535 0.12835204623364604;
            -0.16931095045225217 0.2687956795496601;
            -0.4110745650568839 0.6526163576003452;
            -0.4035951747670475 0.6407421349445884;
            -0.4649907741370211 0.7382129928076485;
            -0.15772708898883683 0.25040532268222715;
            -0.38144358562952446 0.6055745630707645;
            -0.1012787681395544 0.16078875117643368;
            -0.2741403589052255 0.4352214984054432]

y = [4.362866166172215,
            0.8792840060172619,
            1.8414020451091684,
            4.470790758717395,
            4.3894454833815395,
            5.0571760643993455,
            1.7154177874916376,
            4.148527704012107,
            1.1014936742570425,
            2.9815131910316097]

modelqr = lm(X, y; method=:qr)

modelchol = lm(X, y; method=:cholesky)

loglikelihood(modelqr) > loglikelihood(modelchol)
```

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

### Fit Linear Model with AnalyticWeights

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Use `aweights` to fit a linear model with analytic weights, suitable for aggregate values with differing variances.

```jldoctest
julia> using StableRNGs, DataFrames, StatsBase, GLM

julia> data = DataFrame(y = rand(StableRNG(1), 100), x = randn(StableRNG(2), 100), weights = repeat([1, 2, 3, 4], 25))

julia> m_aweights = lm(@formula(y ~ x), data, weights=aweights(data.weights))
LinearModel

y ~ 1 + x

Coefficients:
──────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      t  Pr(>|t|)  Lower 95%  Upper 95%
──────────────────────────────────────────────────────────────────────────
(Intercept)   0.51673     0.0270707  19.09    <1e-34   0.463009  0.570451
x            -0.0478667   0.0308395  -1.55    0.1239  -0.109067  0.0133333
──────────────────────────────────────────────────────────────────────────
```

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

### Fit OLS Model and Extract Statistics

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Demonstrates fitting an Ordinary Least Squares (OLS) model using `lm` and extracting various statistics like coefficients, standard errors, R-squared, AIC, BIC, and log-likelihood. Requires DataFrames and GLM packages.

```julia
using DataFrames, GLM, StatsBase

data = DataFrame(X=[1,2,3], Y=[2,4,7])
olds = lm(@formula(Y ~ X), data)

round.(stderror(ols), digits=5)

round.(predict(ols), digits=5)

round.(confint(ols); digits=5)

round(r2(ols); digits=5)

round(adjr2(ols); digits=5)

deviance(ols)

dof(ols)

dof_residual(ols)

round(aic(ols); digits=5)

round(aicc(ols); digits=5)

round(bic(ols); digits=5)

dispersion(ols)

loglikelihood(ols)

nullloglikelihood(ols)

round.(vcov(ols); digits=5)
```

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

### Handle Categorical Variables and Contrasts in GLM.jl

Source: https://context7.com/juliastats/glm.jl/llms.txt

Shows how GLM.jl automatically handles categorical variables using dummy coding and how to specify custom contrast coding systems like EffectsCoding and HelmertCoding.

```julia
using GLM, DataFrames, CategoricalArrays, StatsModels

data = DataFrame(
    Y = rand(100),
    X = randn(100),
    Group = categorical(repeat(["A", "B", "C", "D"], 25))
)

model = lm(@formula(Y ~ X + Group), data)

model_dummy = lm(@formula(Y ~ X + Group), data;
                 contrasts=Dict(:Group => DummyCoding()))

model_effects = lm(@formula(Y ~ X + Group), data;
                   contrasts=Dict(:Group => EffectsCoding()))

model_helmert = lm(@formula(Y ~ X + Group), data;
                   contrasts=Dict(:Group => HelmertCoding()))

model_ref = lm(@formula(Y ~ X + Group), data;
               contrasts=Dict(:Group => DummyCoding(; base="C")))

data2 = DataFrame(
    Y = rand(100),
    Factor1 = categorical(repeat(["Low", "High"], 50)),
    Factor2 = categorical(repeat(["Control", "Treatment"], inner=50))
)

model_multi = lm(@formula(Y ~ Factor1 * Factor2), data2;
                 contrasts=(
                     :Factor1 => DummyCoding(),
                     :Factor2 => EffectsCoding()
                 ))
```

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

### Fit Linear Model with Weights

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Fit a linear model using `lm` with specified weights. Ensure weights are applied using `aweights` before fitting.

```julia
data.weights = aweights(data.weights)
lm(@formula(y ~ x), data, weights=:weights)
```

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

### Fit Linear Model using `fit`

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/api.md

Use the `fit` function to construct and fit a LinearModel. Requires specifying the model type, design matrix, and response vector.

```julia
using RDatasets

df = RDatasets.dataset("mlmRev", "Oxboys")

fit(LinearModel, hcat(ones(nrow(df)), df.Age), df.Height)
```

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

### Fit Linear Model using `lm`

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/api.md

The `lm` function provides a convenient shorthand for fitting a LinearModel. It takes the design matrix and response vector as arguments.

```julia
lm(hcat(ones(nrow(df)), df.Age), df.Height)
```

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

### Optimize PowerLink Lambda for Linear Regression

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Finds the optimal lambda for a PowerLink transformation in a linear model by minimizing the BIC. Requires GLM, RDatasets, StatsBase, DataFrames, and Optim packages. The optimal lambda is then used to fit the best model.

```julia
using GLM, RDatasets, StatsBase, DataFrames, Optim

trees = DataFrame(dataset("datasets", "trees"));

bic_glm(λ) = bic(glm(@formula(Volume ~ Height + Girth), trees, Normal(), PowerLink(λ)));

optimal_bic = optimize(bic_glm, -1.0, 1.0);

round(optimal_bic.minimizer, digits = 5) # Optimal λ

glm(@formula(Volume ~ Height + Girth), trees, Normal(), PowerLink(optimal_bic.minimizer)) # Best model

round(optimal_bic.minimum, digits=5)
```

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

### Perform Weighted Regression in GLM.jl

Source: https://context7.com/juliastats/glm.jl/llms.txt

Illustrates how to perform linear regression with different types of weights (Frequency, Analytic, Probability, Unit) and how weights affect variance estimation.

```julia
using GLM, DataFrames, StatsBase

data = DataFrame(
    Y = rand(100),
    X = randn(100),
    freq = rand(1:5, 100),      # frequency weights
    prec = rand(100),           # precision/analytic weights
    prob = rand(100)            # sampling probability weights
)

model_unweighted = lm(@formula(Y ~ X), data)

model_freq = lm(@formula(Y ~ X), data; weights=fweights(data.freq))

model_analytic = lm(@formula(Y ~ X), data; weights=aweights(data.prec))

model_prob = lm(@formula(Y ~ X), data; weights=pweights(data.prob))

data.weights = aweights(data.prec)
model_by_name = lm(@formula(Y ~ X), data; weights=:weights)

glM(@formula(Y ~ X), data, Normal(); weights=fweights(data.freq))

coef(model_freq)
coef(model_analytic)
coef(model_prob)

stderror(model_freq)
stderror(model_analytic)
stderror(model_prob)

loglikelihood(model_freq)
loglikelihood(model_analytic)
```

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

### Generate Coefficient Summary Tables in Julia

Source: https://context7.com/juliastats/glm.jl/llms.txt

The `coeftable` function returns a formatted table of coefficient estimates, including standard errors, test statistics, p-values, and confidence intervals. Custom confidence levels can be specified.

```julia
using GLM, DataFrames

data = DataFrame(X1=randn(100), X2=randn(100), Y=randn(100))

# Linear model coefficient table
lm_model = lm(@formula(Y ~ X1 + X2), data)
ct = coeftable(lm_model)
# CoefTable with columns: Coef., Std. Error, t, Pr(>|t|), Lower 95%, Upper 95%

# Access table components
ct.cols      # column data
ct.colnms    # column names
ct.rownms    # row names (coefficient names)

# Custom confidence level
coeftable(lm_model; level=0.99)  # 99% confidence intervals
```

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

### Fit Linear Model with ProbabilityWeights

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Use `pweights` to fit a linear model with probability weights, correcting for under- or over-sampling.

```jldoctest
julia> m_pweights = lm(@formula(y ~ x), data, weights=pweights(data.weights))
LinearModel

y ~ 1 + x

Coefficients:
───────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      t  Pr(>|t|)  Lower 95%   Upper 95%
───────────────────────────────────────────────────────────────────────────
(Intercept)   0.51673     0.0287193  17.99    <1e-32   0.459737  0.573722
x            -0.0478667   0.0265532  -1.80    0.0745  -0.100561  0.00482739
───────────────────────────────────────────────────────────────────────────
```

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

### Fit OLS Model using QR Factorization

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Fits an Ordinary Least Squares (OLS) model using the `lm` function with the QR factorization method. This method is more numerically stable than the default Cholesky factorization for ill-conditioned matrices.

```julia
data = DataFrame(X=[1,2,3], Y=[2,4,7])
olds = lm(@formula(Y ~ X), data; method=:qr)
```

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

### Calculate Log-Likelihood with AnalyticWeights

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Calculate the log-likelihood for a model fitted with analytic weights. Note the filtered output for precision.

```jldoctest
julia> loglikelihood(m_aweights)
-16.296307561384253
```

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

### Fit Negative Binomial Model using glm

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Use the `glm` function with `NegativeBinomial` distribution and `LogLink` to fit a negative binomial regression model. This method requires specifying the dispersion parameter `r` during model creation.

```jldoctest
julia> using GLM, RDatasets

julia> quine = dataset("MASS", "quine")
146×5 DataFrame
 Row │ Eth   Sex   Age   Lrn   Days
     │ Cat…  Cat…  Cat…  Cat…  Int32
─────┼───────────────────────────────
   1 │ A     M     F0    SL        2
   2 │ A     M     F0    SL       11
   3 │ A     M     F0    SL       14
   4 │ A     M     F0    AL        5
   5 │ A     M     F0    AL        5
   6 │ A     M     F0    AL       13
   7 │ A     M     F0    AL       20
   8 │ A     M     F0    AL       22
  ⋮  │  ⋮     ⋮     ⋮     ⋮      ⋮
 140 │ N     F     F3    AL        3
 141 │ N     F     F3    AL        3
 142 │ N     F     F3    AL        5
 143 │ N     F     F3    AL       15
 144 │ N     F     F3    AL       18
 145 │ N     F     F3    AL       22
 146 │ N     F     F3    AL       37
                     131 rows omitted

julia> nbrmodel = glm(@formula(Days ~ Eth+Sex+Age+Lrn), quine, NegativeBinomial(2.0), LogLink())
GeneralizedLinearModel

Days ~ 1 + Eth + Sex + Age + Lrn

Coefficients:
────────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      z  Pr(>|z|)   Lower 95%   Upper 95%
────────────────────────────────────────────────────────────────────────────
(Intercept)   2.88645      0.227144  12.71    <1e-36   2.44125     3.33164
Eth: N       -0.567515     0.152449  -3.72    0.0002  -0.86631    -0.26872
Sex: M        0.0870771    0.159025   0.55    0.5840  -0.224606    0.398761
Age: F1      -0.445076     0.239087  -1.86    0.0627  -0.913678    0.0235251
Age: F2       0.0927999    0.234502   0.40    0.6923  -0.366816    0.552416
Age: F3       0.359485     0.246586   1.46    0.1449  -0.123814    0.842784
Lrn: SL       0.296768     0.185934   1.60    0.1105  -0.0676559   0.661191
────────────────────────────────────────────────────────────────────────────
```

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

### Model Methods

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/api.md

Provides a list of methods available for working with fitted GLM.jl models, including diagnostics and prediction.

```APIDOC
## Model methods

- `cooksdistance(model)`: Computes Cook's distance for a model.
- `StatsBase.deviance(model)`: Calculates the deviance of a model.
- `GLM.dispersion(model)`: Returns the dispersion parameter of a model.
- `GLM.ftest(model)`: Performs an F-test on the model.
- `StatsBase.nobs(model)`: Returns the number of observations in the model.
- `StatsBase.nulldeviance(model)`: Calculates the null deviance of a model.
- `StatsBase.predict(model, X)`: Predicts values using the fitted model.
```

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

### Fit Linear Model with Categorical Variable (DummyCoding)

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Fit a linear model using `lm` with a categorical variable specified using explicit DummyCoding. This method is useful when not using DataFrames or when a specific contrast coding is desired.

```julia
using StableRNGs

data = DataFrame(y = rand(StableRNG(1), 100), x = repeat([1, 2, 3, 4], 25))

lm(@formula(y ~ x), data, contrasts = Dict(:x => DummyCoding()))

```

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

### Model Methods API

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Provides access to various methods applicable to fitted statistical models.

```APIDOC
## GET /api/models/{model_id}/methods

### Description
Retrieves results from common statistical methods applied to a fitted model.

### Method
GET

### Endpoint
/api/models/{model_id}/methods

### Parameters
#### Path Parameters
- **model_id** (string) - Required - The unique identifier for the fitted model.

#### Query Parameters
- **method** (string) - Optional - The specific method to call (e.g., `coef`, `r2`, `aic`, `predict`, `cooksdistance`). If omitted, returns a list of available methods.
- **newX** (object) - Optional - For the `predict` method, provides new covariate values for prediction.

### Response
#### Success Response (200)
- **method_name** (string) - The name of the method called.
- **result** (any) - The result of the method call. The type depends on the method.

#### Response Example (for `coef` method)
```json
{
  "method_name": "coef",
  "result": [
    -0.66666667,
    2.5
  ]
}
```

#### Response Example (for `r2` method)
```json
{
  "method_name": "r2",
  "result": 0.98684211
}
```

#### Response Example (for `aic` method)
```json
{
  "method_name": "aic",
  "result": 5.84251593
}
```

#### Response Example (for `predict` method with `newX`)
```json
{
  "method_name": "predict",
  "result": [
    9.33333333
  ]
}
```

#### Response Example (for `cooksdistance` method)
```json
{
  "method_name": "cooksdistance",
  "result": [
    2.5,
    0.25,
    2.5
  ]
}
```
```

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

### Calculate Log-Likelihood with FrequencyWeights

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Calculate the log-likelihood for a model fitted with frequency weights. Note the filtered output for precision.

```jldoctest
julia> loglikelihood(m_fweights)
-25.518609617564483
```

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

### Compare Three Models with F-test

Source: https://context7.com/juliastats/glm.jl/llms.txt

Use the `ftest` function to compare three models. This can be done in either forward or backward direction.

```julia
ftest(null_model, treatment_model, full_model)
```

```julia
ftest(full_model, treatment_model, null_model)
```

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

### Enable GLM Debugging Output

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Set the JULIA_DEBUG environment variable to GLM to enable detailed output from the fitting steps, useful for diagnosing failed model fits.

```julia
ENV["JULIA_DEBUG"] = "GLM"
```

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

### Linear Regression with Cholesky Method

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/examples.md

Demonstrates linear regression using the Cholesky method, which can detect multicollinearity. Note that multicollinearity may lead to NaN values in coefficients.

```julia
julia> y = [9.42190647125244, 2.084805727005, 3.9376676082611, 2.61976027488708, 4.04761934280395, 2.15384602546691,
            2.66240668296813, 4.39475727081298, 5.74520826339721, 3.59616208076477, 1.54265284538269, 2.59368276596069,
            1.80476510524749, 1.69270837306976, 3.04201245307922, 2.18389105796813, 2.73844122886657, 2.88134002685546,
            2.46666669845581, 3.80616021156311, 5.12149810791015, 6.80378007888793, 3.73669862747192, 1.21332454681396,
            2.54629635810852, 5.1612901687622, 1.86798071861267, 1.21465551853179, 6.31019830703735, 1.02669405937194, 
            2.50623273849487, 1.5936255455017];

julia> x = [4570.2001953125, 2830, 596.799987792968, 133.600006103515, 42, 390, 93.9000015258789, 907.900024414062,
            19773, 39709, 2936.5, 2513.80004882812, 1124.80004882812, 921.599975585937, 2432.60009765625, 6754,
            1066.30004882812, 3199.89990234375, 150, 509.700012207031, 1452.69995117187, 8995, 1212.30004882812,
            906.599975585937, 2592, 201.5, 2617.80004882812, 502.200012207031, 2824, 292.200012207031, 7621, 1631.5];

julia> rdchem = DataFrame(rdintens=y, sales=x);

julia> mdl = lm(@formula(rdintens ~ sales + sales^2), rdchem; method=:cholesky)
LinearModel

rdintens ~ 1 + sales + :(sales ^ 2)

Coefficients:
───────────────────────────────────────────────────────────────────────────────────────
                    Coef.     Std. Error       t  Pr(>|t|)      Lower 95%     Upper 95%
───────────────────────────────────────────────────────────────────────────────────────
(Intercept)   0.0          NaN            NaN       NaN     NaN            NaN
sales         0.000852509    0.000156784    5.44    <1e-05    0.000532313    0.00117271
sales ^ 2    -1.97385e-8     4.56287e-9    -4.33    0.0002   -2.90571e-8    -1.04199e-8
───────────────────────────────────────────────────────────────────────────────────────
```

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

### Fit Linear Model with FrequencyWeights

Source: https://github.com/juliastats/glm.jl/blob/master/docs/src/index.md

Use `fweights` to fit a linear model with frequency weights, representing the number of times each observation was seen.

```jldoctest
julia> m_fweights = lm(@formula(y ~ x), data, weights=fweights(data.weights))
LinearModel

y ~ 1 + x

Coefficients:
─────────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      t  Pr(>|t|)   Lower 95%    Upper 95%
─────────────────────────────────────────────────────────────────────────────
(Intercept)   0.51673     0.0170172  30.37    <1e-84   0.483213    0.550246
x            -0.0478667   0.0193863  -2.47    0.0142  -0.0860494  -0.00968394
─────────────────────────────────────────────────────────────────────────────
```

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

### Displaying GLM Coefficient Table

Source: https://context7.com/juliastats/glm.jl/llms.txt

Prints a formatted coefficient table for a linear model, showing coefficients, standard errors, t-statistics, p-values, and confidence intervals. This is useful for a quick overview of model results.

```julia
println(coeftable(lm_model))
```