How do I fix this convergence error? Python 3 statsmodels

user_8465696 0
Viewed 2067

This question pertains to Python 3 statsmodels and its general linear model class.

Whenever I have an array of values for my endogenous variable such that the values are more than an order of magnitude apart, the GLM does not converge and it throws an exception. Here is a coded example of what I mean.

import pandas as pd
import pyarrow.parquet as pq
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import math 

col = ["a", \
       "b", \
       "c", \
       "d", \
       "e", \
       "f", \
       "g", \
       "h"]

df = pd.DataFrame(np.random.randint(low=1, high=100, size=(20, 8)), columns=col)
df["a"] = [0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01]
df2 = pd.DataFrame(np.random.randint(low=1, high=100, size=(20, 8)), columns=col)
df2["a"] = np.random.randint(low=10000, high=99999, size=(20, 1))
df3 = pd.DataFrame(np.random.randint(low=1, high=100, size=(20, 8)), columns=col)
df3["a"] = [0.01, \
            0.01, \
            0.01, \
            0.01, \
            0.01, \
            0.01, \
            np.random.randint(low=10000, high=99999), \
            0.01, \
            0.01, \
            0.01, \
            0.01, \
            0.01, \
            np.random.randint(low=10000, high=99999), \
            0.01, \
            0.01, \
            0.01, \
            0.01, \
            0.01, \
            0.01, \
            0.01]
try:
    actual = df[["a"]]

    fml1 = "a ~ log(b) + c + d + e + f + g"

    data1 = df[["b", "c", "d", "e", "f", "g"]]

    model = sm.GLM(actual, data1, formula=fml1, family=sm.families.Tweedie(link_power=1.1)).fit()
    model_pred = model.predict()
    print("SUCCESS")
except:
    print("FAILURE")
try:
    actual = df2[["a"]]

    fml1 = "a ~ log(b) + c + d + e + f + g"

    data1 = df2[["b", "c", "d", "e", "f", "g"]]

    model = sm.GLM(actual, data1, formula=fml1, family=sm.families.Tweedie(link_power=1.1)).fit()
    model_pred = model.predict()
    print("SUCCESS")
except:
    print("FAILURE")
try:
    actual = df3[["a"]]

    fml1 = "a ~ log(b) + c + d + e + f + g"

    data1 = df3[["b", "c", "d", "e", "f", "g"]]

    model = sm.GLM(actual, data1, formula=fml1, family=sm.families.Tweedie(link_power=1.1)).fit()
    model_pred = model.predict()
    print("SUCCESS")
except:
    print("FAILURE")

If you run this code, you should get an exception only on the last set of data. Why is this? How do I get the GLM to converge? Are there alternatives?

statsmodels pandas matplotlib numpy python
1 Answers

It seems that fitting parameters of Tweedie distributions is not easy. Indeed, a set of parameter w is valid only if for all observations x_i the positivity of the dot products is verified, i.e. enter image description here, otherwise the value that is used in prediction is undefined as a negative real number cannot be elevated to a non-integer power value.

Hence in most optimizers, this relationship should be maintained across all iterations and it can be tough to maintain especially if the data contains values with different orders of magnitude.

Then, I see two main solutions to deal with this problem

  • The easiest one: Force your coefficients w to be positive. As in your case all observations x_i are positive you will remain in the feasible set. This can be done using Newton solver and a callback for example:

        model = sm.GLM(actual, data1, formula=fml1,
                   family=sm.families.Tweedie(link_power=1.1))

    def callback(x):
        x[x < 0] = 0

    result = model.fit(method='newton', disp=True, start_params=np.ones(6),
                       callback=callback)

    This will converge every time but will reach a solution in which all coefficients are positive, i.e. there is no inhibition effect.

  • Another solution might be looking at the conjugate solvers. For some reasons, they perform better with these constraints. This can be done using Conjugate Gradient "cg" and Newton Conjugate Gradient "ncg" methods. They might not converge every time but chances are they will. You can play with the start_params vector it doesn't but it is not an exact science. In your case, you can try the following setting:

    model = sm.GLM(actual, data1, formula=fml1,
               family=sm.families.Tweedie(link_power=1.1))
result = model.fit(method='cg', disp=True, start_params=0.1 * np.ones(6))

PS: I'm not an expert on Tweedie distributions but I have been working on other Poisson like distributions that are facing the same problems.

user_20157620
Related