How can I estimate Poisson quasi-MLE models in Python?

Question

How do I estimate Poisson quasi maximum likelihood (MLE) models in python? Standard Poisson MLE assumes that the conditional mean of the variable equals its variance. Quasi-MLE permits relaxes this assumption. Following this example,

#load packages

import numpy as np
import statsmodels.api as sm
from scipy import stats
from matplotlib import pyplot as plt

 #load sample data

 data = sm.datasets.star98.load()
 data.exog = sm.add_constant(data.exog, prepend=False)

 #Estimate Poisson MLE
 glm_poisson = sm.GLM(data.endog.NABOVE, data.exog, family=sm.families.Poisson())
 res = glm_poisson.fit(scale="X2")
 print(res.summary())

I'd like to do something like the below, which works for the binomial distribution, except for a more general Poisson. That is, define a different variance (e.g. conditional variance=a*(conditional mean), where a is a parameter to be estimated) and estimate a model.

     class vf(sm.families.varfuncs.VarianceFunction):
         def __call__(self, mu):
         return mu ** 2 * (1 - mu) ** 2

        def deriv(self, mu):
        return 2 * mu - 6 * mu ** 2 + 4 * mu ** 3

    bin = sm.families.Binomial()
    bin.variance = vf()
    model2 = sm.GLM.from_formula("blotch ~ 0 + C(variety) + C(site)", family=bin, data=df)
    result2 = model2.fit(scale="X2")
    print(result2.summary())

I knew this wasn't going to work, but to give you the general idea, I'd like to make the Poisson() variance more general

 class vf(sm.families.varfuncs.VarianceFunction):
     def __call__(self, mu):
     return mu *z

     def deriv(self, mu):
     return z
 pois = sm.families.Poisson()
 pois.variance=vf()
 glm_poisson1 = sm.GLM(data.endog.NABOVE, data.exog, family=sm.families.pois())
 res = glm_poisson1.fit(scale="X2")
 print(res.summary())