Tengo la siguiente regresión lineal:

import statsmodels.formula.api as sm

model = sm.ols(formula = 'a ~ b + c', data = data).fit()

Quiero agregar un término cuadrático para b en este modelo.

¿Hay una manera simple de hacer esto con statsmodels.ols? ¿Hay un paquete mejor que debería usar para lograr esto?