Support for multinomial logistic regression

[pzivich]

Is your feature request related to a problem? Please describe. Unrelated to problem.

Describe the solution you'd like Addition of estimating equations for multinomial logistic regression. This is currently a major missing component for the regression support of deli.

Ideally, I would integrate into ee_regression as an option. However, y is a bit of a unique matter since it takes a set of values and predicted probabilities are for each level. By processing y from the categories into a dummy variable behind the scenes, it becomes opaque to users. Instead, I am going to force users to dummy code y outside of the built-in estimating equation. This aligns the predictions a bit better. So, multinomial logit models are going to be separate ee_... objects

This addition will also require reconsidering how the prediction function works. A vector of predictions is no longer sufficient. We will need a matrix instead.

Describe alternatives you've considered Composing the logistic models together (doesn't work because it is a different model). We need new estimating functions for this.

Additional context Below is some rough code to implement. Note that y has a different data structure from the other regression models.

def ee_mlogit(theta, X, y):
    X = np.asarray(X)
    y = np.asarray(y)
    theta = np.asarray(theta)
    n_y_vals = y.shape[1]
    n_x_vals = X.shape[1]

    # Compute the denominator
    start_index = 0
    exp_pred_y = []
    denom = 1
    for i in range(1, n_y_vals):
        end_index = start_index + n_x_vals
        beta_i = theta[start_index: end_index]
        pred_y = np.exp(np.dot(X, beta_i))
        denom = denom + pred_y
        exp_pred_y.append(pred_y)
        start_index = end_index

    efuncs = []
    yhat_ref = y[:, 0] - 1/denom
    for i in range(1, n_y_vals):
        yhat_i = yhat_ref + (y[:, i] - exp_pred_y[i-1]/denom)
        residual = yhat_i[:, None] * X
        efuncs.append(residual.T)

    return np.vstack(efuncs)