GLM.jl

The lm function uses Cholesky factorization by default, which is fast but can be numerically unstable for ill-conditioned design matrices or aggressive in detecting multicollinearity.

To improve numerical stability or handle multicollinearity, use the method keyword argument:

• method=:cholesky (Default): Fast, but sensitive to ill-conditioned matrices and multicollinearity.
• method=:qr: More stable for ill-conditioned design matrices and less aggressive in dropping collinear columns.

If you encounter models where coefficients appear as NaN or parameters are not estimated due to multicollinearity detection, try switching to method=:qr.

# Using QR decomposition for better stability
ols = lm(@formula(Y ~ X), data; method=:qr)

# Using Cholesky decomposition (default)
ols = lm(@formula(Y ~ X), data; method=:cholesky)

To install GLM.jl and its dependencies (including the Distributions package), use the Julia package manager:

Pkg.add("GLM")

If a generalized linear model fit fails, you can enable detailed iteration output by setting the JULIA_DEBUG environment variable to GLM. This will output the deviance and the change in deviance for each iteration of the fitting process.

# Run this in your terminal before starting Julia, or set it within your environment
ENV["JULIA_DEBUG"] = "GLM"

Negative binomial regression can be performed in two ways using the GLM.jl package:

1. Using glm: Pass NegativeBinomial(theta) as the family argument to glm, where theta is the dispersion parameter.
2. Using negbin: A specialized convenience function negbin that estimates the dispersion parameter automatically.

Both methods require a link function, such as LogLink().

using GLM, RDatasets

quine = dataset("MASS", "quine")

# Method 1: Using glm with a specified theta
nbrmodel = glm(@formula(Days ~ Eth+Sex+Age+Lrn), quine, NegativeBinomial(2.0), LogLink())

# Method 2: Using the negbin convenience function (estimates theta)
nbrmodel = negbin(@formula(Days ~ Eth+Sex+Age+Lrn), quine, LogLink())

# Accessing the estimated theta from the model
println("Estimated theta = ", round(nbrmodel.rr.d.r, digits=5))

To perform a Probit regression, use the glm function with a Binomial() family and a ProbitLink() link function. This is typically used when the response variable is binary (0 or 1).

using GLM, DataFrames

data = DataFrame(X=[1,2,2], Y=[1,0,1])
probit = glm(@formula(Y ~ X), data, Binomial(), ProbitLink())

For convenience, you can store pre-constructed weight objects directly in a DataFrame column. You can then pass the column name as a symbol to the weights keyword argument in lm or glm.

# Pre-construct the weights in the DataFrame
data.weights = aweights(data.weights);

# Pass the column name as a symbol
m_aweights = lm(@formula(y ~ x), data, weights=:weights)

You can use PowerLink(λ) as a link function in a glm call to model data where the relationship between the mean and the linear predictor follows a power law. To find the optimal power parameter λ, you can wrap the bic_glm calculation (using the bic function) in an optimization routine like Optim.optimize to minimize the Bayesian Information Criterion (BIC).

using GLM, RDatasets, StatsBase, DataFrames, Optim

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

# Define a function to calculate BIC for a given λ
bic_glm(λ) = bic(glm(@formula(Volume ~ Height + Girth), trees, Normal(), PowerLink(λ)));

# Optimize λ to minimize BIC
optimal_bic = optimize(bic_glm, -1.0, 1.0);

# Use the optimal λ to fit the best model
optimal_λ = round(optimal_bic.minimizer, digits = 5)
model = glm(@formula(Volume ~ Height + Girth), trees, Normal(), PowerLink(optimal_λ))

Once a model is fitted, you can use several methods to inspect the results, perform predictions, or check model fit:

• StatsBase.predict(model, X): Predict response values for new data X.
• StatsBase.nobs(model): Returns the number of observations.
• StatsBase.deviance(model): Returns the deviance.
• StatsBase.nulldeviance(model): Returns the null deviance.
• GLM.ftest(model): Performs an F-test.
• GLM.dispersion(model): Returns the dispersion parameter.
• cooksdistance(model): Calculates Cook's distance for each observation.