irls Subroutine

public subroutine irls(method, x, y, beta, weights, options, info, status, resid, stats, alpha, cov)

Fit a GLM from initial coefficients. The working response and weights are z = eta + (y - mu) * g'(mu) and W = 1 / (g'(mu)**2 * V(mu)). The implementation solves the equivalent increment problem sqrt(W) X delta = sqrt(W) (z - eta), using column-scaled, pivoted QR rather than normal equations. Weights and right-hand sides are normalized in log space to avoid overflow.

No intercept, offset, binomial trial counts, or prior observation weights are added. A common dispersion factor cancels from updates. Robust weights act on raw residuals without scale estimation; set options%use_robust_weighting false for ordinary likelihood fitting.

Initial predictors must be in the family's domain. Proposed updates are halved up to 64 times to preserve that domain and avoid increasing deviance (with current robust weights held fixed). Only an unhalved step can declare convergence. There is no guarantee of convergence, especially for separated binomial data. The iteration limit returns the last valid coefficients; supply info to distinguish convergence from an exhausted iteration budget.

Size mismatches stop with FS_MATRIX_SIZE_ERROR. Empty or nonfinite inputs, invalid family/options/responses/callbacks, unrepresentable arithmetic, and invalid mean domains stop with FS_INVALID_INPUT_ERROR. Too few observations or numerical rank deficiency of the weighted, scaled design stop with FS_UNDERDEFINED_PROBLEM_ERROR. Missing data must be removed or imputed by the caller.

Arguments

Type IntentOptional Attributes Name
integer(kind=int32), intent(in) :: method

FS_GLM_BINOMIAL (logit), FS_GLM_POISSON (log), or FS_GLM_GAMMA (reciprocal).

real(kind=real64), intent(in), dimension(:,:) :: x

Finite N-by-P design, N >= P > 0. Weighted columns must be numerically independent at epsilon * max(N,P) relative scale.

real(kind=real64), intent(in), dimension(:) :: y

N finite responses: [0,1] for binomial, >= 0 for Poisson,

0 for Gamma. Fractional responses are allowed; binomial proportions are treated as unit-trial quasi-likelihood data.

real(kind=real64), intent(inout), dimension(:) :: beta

P finite starting coefficients; receives the last valid iterate. Zero is valid for binomial/Poisson, but Gamma requires X beta > 0.

procedure(robust_weight_function), intent(in), optional, pointer :: weights

Robust callback used only when enabled. Absent/disassociated pointers select Tukey's biweight on raw response residuals.

type(glm_options), intent(in), optional :: options

Solver controls; omission uses the documented component defaults.

type(glm_convergence_info), intent(out), optional :: info

Convergence flag, completed update count, and last maximum absolute coefficient change (not the response residual).

procedure(glm_iteration_update), intent(in), optional, pointer :: status

An optional pointer to a routine that can be used to extract iteration information.

real(kind=real64), intent(out), optional, dimension(:) :: resid

An optional P-element array where the residual will be written.

type(glm_coefficient_statistics), intent(out), optional, dimension(:) :: stats

P coefficient statistics, computed only after convergence; otherwise all fields are NaN. Ordinary binomial/Poisson fits use unit-dispersion information covariance and normal Wald inference. Ordinary Gamma fits use Pearson dispersion and approximate t inference with N-P degrees of freedom. Robust fits use an HC1 empirical sandwich and normal Wald inference, treating rows as independent. Robust callbacks must be deterministic, smooth near the solution, and observation-wise for this inference to apply. Gamma and robust inference require N > P. Fractional responses do not automatically enable overdispersion estimation.

real(kind=real64), intent(in), optional :: alpha

Finite significance level in (0,1), default 0.05. Intervals have confidence level 1-alpha. Validated even when stats is omitted.

real(kind=real64), intent(out), optional, dimension(:,:) :: cov

P-by-P coefficient covariance in original beta order and units. Uses the same final-model information/Pearson or HC1 sandwich calculation as stats. May be requested independently of stats. All entries are NaN if fitting does not converge. Gamma and robust covariance require N > P; unrepresentable entries stop with FS_INVALID_INPUT_ERROR. Not a response covariance matrix.