Fit a binomial or multinomial logistic regression model using iteratively re-weighted least squares.

The type double function is imsls_d_logistic_regression.

Pointer to an array of length n_coefficients × n_classes containing the estimated coefficients. The function fits a full model, where n_coefficients = 1 + n_independent. The optional arguments IMSLS_NO_INTERCEPT, IMSLS_X_INDICES, and IMSLS_X_INTERACTIONS may be used to specify different models. Note that the last column (column n_classes) represents the reference class and is set to all zeros.

where

and

are unknown coefficients that are to be estimated.

Solving for the linear component η1 results in the log-odds or logit transformation of π1 (x):

Given a set of N observations (yi, xi), where yi follows a binomial (n, π) distribution with parameters n = 1 and π = π1 (xi), the likelihood and log-likelihood are, respectively,

The log-likelihood in terms of the parameters, {β01, β1}, is therefore

where

With a binary outcome, only one probability needs to be modeled. The second probability can be obtained from the constraint, π1 (x) + π2(x) = 1. If each yi is the number of successes in ni independent trials, the log-likelihood becomes

or

See optional argument IMSLS_FREQUENCIES to set frequencies ni > 1.

To test the significance of the model, the log-likelihood of the fitted model is compared to that of an intercept‑only model. In particular, G = ‑2 (l(β01) - l(β01, β1)) is a likelihood-ratio test statistic and under the null hypothesis, H0 : β11 = β12 = … = β1p = 0, G is distributed as chi-squared with p‑1 degrees of freedom. A significant result suggests that at least one parameter in the model is non‑zero. See Hosmer and Lemeshow (2000) for further discussion.

In the multinomial case, the response vector is yi = (yi1, yi2, …, yiK)', where yik = 1 when the i-th observation belongs to class k and yik = 0, otherwise. Furthermore, because the outcomes are mutually exclusive,

and π1 (x) + π2 (x) +--- + πK (x) = 1. The last class K serves as the baseline or reference class in the sense that it is not modeled directly but found from

If there are multiple trials, ni > 1, then the constraint on the responses is

The log‑likelihood in the multinomial case becomes

or

The constraint

is handled by setting ηK = 0 for the K‑th class, and then the log‑likelihood is

Note that for the multinomial case, the log‑odds (or logit) is

Note that each of the logits involve the odds ratio of being in class l versus class K, the reference class. Maximum likelihood estimates can be obtained by solving the score equation for each parameter:

To solve the score equations, the function employs a method known as iteratively re-weighted least squares or IRLS. In this case the IRLS is equivalent to the Newton-Raphson algorithm (Hastie, et. al., 2009, Thisted, 1988).

Consider the full vector of parameters

the Newton-Raphson iteration is

where H denotes the Hessian matrix, i.e., the matrix of second partial derivatives defined by

and

and G denotes the gradient vector, the vector of first partial derivatives,

Both the gradient and the Hessian are evaluated at the most recent estimate of the parameters, βn. The iteration continues until convergence or until maximum iterations are reached. Following the theory of maximum likelihood estimation (Kendall and Stuart, 1979), standard errors are obtained from Fisher’s information matrix (-H)-1 evaluated at the final estimates.

When the IMSLS_NEXT_RESULTS option is specified, the function combines estimates of the same model from separate fits using the method presented in Xi, Lin, and Chen (2008). To illustrate, let β1 and β2 be the MLE’s from separate fits to two different sets of data, and let H1 and H2 be the associated Hessian matrices. Then the combined estimate,

approximates the MLE of the combined data set. The model structure, Imsls_f_model **next_model contains the combined estimates as well as other elements. See Table 1: Imsls_f_model Data Structure below.

Iteration stops when the estimates converge within tolerance, when maximum iterations are reached, or when the gradient becomes within tolerance of 0, whichever event occurs first. When the gradient converges before the coefficient estimate converges, a condition in the data known as complete or quasi-complete separation may be present. Separation in the data means that one or more independent variable perfectly predicts the response. When detected, the function stops the iteration, issues a warning, and returns the current values of the model estimates. Some of the coefficient estimates and standard errors may not be reliable. Furthermore, overflow issues may occur before the gradient converges. In such cases the program issues a fatal error.

The first example is from Prentice (1976) and involves the mortality of beetles after five hours exposure to eight different concentrations of carbon disulphide. The table below lists the number of beetles exposed (N) to each concentration level of carbon disulphide (x, given as log dosage) and the number of deaths which result (y):

The number of deaths at each concentration level is the binomial response (n_classes = 2) and the log-dosage is the single independent variable. Note that this example illustrates the GROUP_COUNTS format for y and the optional argument IMSLS_FREQUENCIES.

In this example the response is a multinomial random variable with 4 outcome classes. The 3 independent variables represent 2 categorical variables and 1 continuous variable. A subset of 2 independent variables along with the intercept defines the logistic regression model. A test of significance is performed.

Example 3 uses the same data as in Example 2 and an additional set of 50 observations using the same data generating process. The model structure includes all 3 independent variables and an intercept, and a single model fit is approximated from two separate model fits. Example 3 also includes a fit on the full data set for comparison purposes.