Inverted prediction with regressor?

[ghost]

Hi all,

First of all, thanks for your work on this awesome repository!

I'm having some issues with an 'exploding' forecast that seems to be kind of inverted. I'm creating a Prophet model for 'Object A' that doesn't have a lot of data. Therefore, there's no seasonality found in the data and thus only a trend available. To compensate for this, I'm extracting the seasonality component from a model that was trained on 'Group A' to which 'Object A' belongs (and which has more data available). I'm then adding the seasonality output from the 'Group A' model as a regressor to the model of 'Object A'.

So to summarise: The seasonality of 'Group A' (to which 'Object A' belongs) is added as a regressor to 'Object A'. You can find a small schematic overview of this below.

You can see the prediction of 'Object A', together with the components, in the screenshot below.

(Left top graph is the trend of 'Object A'; left bottom graph is the regressor and thus the seasonality of 'Group A'; right graph is the prediction of 'Object A')

The issue here is that the prediction is kind of 'exploding' in the wrong direction. You can see that the extra regressor is going up around 2021-03/2021-04, but the actual prediction is going down. Around 2021-06, the regressor is going down again, while the prediction is actually going up.

I don't really understand how this is happening as this way of working (extracting seasonality data from a model and using it in another model) is working fine for all my other models. I suppose it has something to do with the negative trend, but I still don't feel like it is expected behaviour?

The regressor is added to the model using the following code:

model.add_regressor('Regressor Name', mode='multiplicative')

The model itself is simply created with the following code:

model = Prophet()

Extra remarks

• The data itself is of weekly format.
• The data can not go into negative (should always be positive). I'm currently using the 'clipping' approach from Strategies for positive predictions #1668 where I simply clip predicted negative values to 0.
• I'm using Prophet 0.7.1 together with Python 3.7
• The seasonality which I'm using of 'Group A' is the following:
  

Thank you in advance for any input and/or solution!

[ghost]

I would like to add an additional example. The concept is exactly the same as explained in my main post, but it's for a different object, let's call it 'Object B', and this object belongs to a different group called 'Group B'.

You can see the prediction of 'Object B' and its component below:

(Left top graph is the trend of 'Object B'; left bottom graph is the regressor and thus the seasonality of 'Group B'; right graph is the prediction of 'Object B')

The prediction kind of 'exploded' here (predicting a very high number), but I found that I can fix this by using a lower number for prior_scale when adding the regressor.

The seasonality of 'Group B' that was used in the model of 'Object B' can be found below:

The strange thing here (and what makes it differ from the first example in my original post) is that Prophet seems to have deliberately 'inverted' the seasonality values of 'Group B' (the graph of the seasonality of 'Group B' is the inverse of the regressor that was used in 'Object B'). So it feels like Prophet noticed that the trend was (almost) completely negative and inverted the values used in the regressor such that - and - would become + making the actual predicted values correct? And this behaviour is not happening in the example that is in the original post..

[hansukyang]

I think it would help to apply a bit of common sense here and consider that there is really no magic behind these tools. In the first example, the model is trying to forecast one year ahead using only six months of data. In the second, there seems to be only two data points for one year of prediction - really not enough data points to generate any meaningful prediction.

As a comparison, one-day ahead hourly forecasts for electrical grid use models trained on three years worth of data - training period to forecast horizon ratio of 1000 (365 days x 3 years) to 1.

[ghost]

Thanks for your answer @hansukyang ! I know that the data is limited in some cases, and this is exactly the reason why I tried to introduce the regressor which is based on data of at least a few years. I'm really not expecting a perfect prediction with only a few datapoints, but I'm finding the behaviour of the regressor in case of a negative trend very strange, which is why I opened this issue.

[hansukyang]

I see, thanks for the clarification. I saw similar behaviour when I had small number of data points and I explained to myself that it's probably like fitting a parabola to two points. For my own use, I find the linear growth trend somewhat misleading so I was very happy when growth='flat' option was introduced recently.

[bletham]

@JeremyKeustersML6 I think the issue here is a little subtle, and it isn't one that has come up much before.
The multiplicative model for Prophet is basically

y(t) = trend(t) * (1 + seasonality(t) + beta * regressor(t)) + noise

so the magnitude of the regressor effect is trend(t) * beta * regressor(t) and the % effect as shown in the components plot is 100% * beta * regressor(t).

The key thing to note here is that the sign of the magnitude of the regressor effect will flip when the sign of the trend flips. Let's take beta>0 for now. Then, with a positive trend an increase in the regressor will produce an increase in y. With a negative trend, an increase in the regressor will produce a decrease in y. This is what you're seeing.

The core of the issue is really the fact that we are getting a negative trend when we know that is impossible. The clipping approach from #1668 ensures positive yhat, but doesn't ensure positive trend and so we have this issue. I'd be really interested to see how the ProphetPos class that I proposed there works in this time series; it would ensure the trend stays positive so the regressor sign would not flip (it would in all likelihood just stay really close to 0).

The other thing you could try would be an additive regressor, instead of multiplicative. Then the model is

y(t) = trend(t) + seasonality(t) + beta * regressor(t) + noise

and the effect of the regressor is no longer changed by the sign of the trend.

[ghost]

Thank you both @hansukyang and @bletham for your answers.

For my own use, I find the linear growth trend somewhat misleading so I was very happy when growth='flat' option was introduced recently.

This is an interesting option I didn't know about, thanks.

---

The core of the issue is really the fact that we are getting a negative trend when we know that is impossible. The clipping approach from #1668 ensures positive yhat, but doesn't ensure positive trend and so we have this issue. I'd be really interested to see how the ProphetPos class that I proposed there works in this time series; it would ensure the trend stays positive so the regressor sign would not flip (it would in all likelihood just stay really close to 0).

Thanks for the good explanation and the suggestion. I tried the ProphetPos on Object A and B:

ProphetPos on Object A:

ProphetPos on Object B:

There are two issues with these results:

• The class indeed seems to be fixing the negative trend issue, although there's still a small negative trend peak in Object A?
• The issue in my use-case is that the trend stays on 0, making the regressor redundant. I'm adding this regressor specifically for items that don't have a lot of datapoints (and no seasonality yet) so that they can still have some kind of prediction, so ideally I would have a solution where the trend remains slightly higher than 0. I tried changing this myself in the ProphetPos class, but as soon as I change the 0 values to a different number, it seems to break.

---

An another note, I also changed my code such that the seasonality/forecast of 'Group A' and 'Group B' is forecasted with seasonality_mode='multiplicative' instead of with the default additive (makes more sense). When doing this, the first issue of the negative trend on Object A is already smaller, but of course the second issue remains. You can see the results below.

ProphetPos on Object A (with multiplicative seasonality when forecasting Group A, of which its seasonality is used as regressor in the forecast of Object A):

(So in an ideal case scenario, there would be a small positive prediction around 05-2021)

ProphetPos on Object B (with multiplicative seasonality when forecasting Group B, of which its seasonality is used as regressor in the forecast of Object B):

---

The other thing you could try would be an additive regressor, instead of multiplicative. Then the model is

y(t) = trend(t) + seasonality(t) + beta * regressor(t) + noise

and the effect of the regressor is no longer changed by the sign of the trend.

I also tried this, but it performs in general worse on my data.

For Object A, it pushes the prediction into negative numbers without any perspective to go to positive again:

For Object B, it actually suddenly generates a positive trend:

[bletham]

Thanks for sharing those results. Yeah it's a bit tricky here since the trend does go to 0 early on, and then it ends up being a bit sticky at 0. I think @hansukyang's suggestion of a flat trend would be probably be the best chance then at improving the fit since that would force the model to treat the initial decrease as coming from the regressor, and the trend would stay positive.

[ghost]

Thanks again for your answer @bletham ! I'll do some experiments with the flat trend and see if it leads to something!

Regarding the ProphetPos, would it be technically possible to change this class such that the trend stays at a value slightly higher than 0?

[bletham]

I think in these lines:

        while min(trend[indx_future:]) < 0:
            indx_neg = indx_future + np.argmax(trend[indx_future:] < 0)

if you replace 0 with something greater than 0 that would do it.
However, the trend is fit in a scaled space, where everything has been divided by m.y_scale. So if you want the trend to saturate at 5, then you would need to in the trend code above set it to 5 / m.y_scale (except m isn't being passed into that function, so you'd need to hack it in somewhere). I think the combination of this class with an additive regressor may be worth giving a shot too.

[ghost]

I think in these lines:

        while min(trend[indx_future:]) < 0:
            indx_neg = indx_future + np.argmax(trend[indx_future:] < 0)

if you replace 0 with something greater than 0 that would do it.
However, the trend is fit in a scaled space, where everything has been divided by m.y_scale. So if you want the trend to saturate at 5, then you would need to in the trend code above set it to 5 / m.y_scale (except m isn't being passed into that function, so you'd need to hack it in somewhere). I think the combination of this class with an additive regressor may be worth giving a shot too.

Thanks again for your answer! I indeed tried something like this already in the past, but the program would just be stuck on some kind if infinite loop, so that probably has to do with the m.y_scale you're talking of. I'll try to hack m into that function, thanks a lot for this!

---

On the other hand, I was pleasantly surprised with the results that I got when using the growth='flat' solution from @hansukyang !

For Object A:

For Object B:

The results make sense, as I wanted more 'weight' on the regressor. It also actually didn't make any sense to get a real trend based on only a few datapoints, so this is kind of solved here. I will still experiment by adding eg. the trend of 'Group A' and 'Group B' respectively (instead of only adding their seasonalities) to still have some kind of 'real trend' in there!

---

Thank you again to @bletham and @hansukyang ! You've been both really helpful! I might still update this issue in the future with additional findings I encounter.

[ghost]

Just a small update, I succeeded incorporating the y_scale parameter in the following way.

I renamed the existing piecewise_linear method in the ProphetPos class to piecewise_linear_scale (since we need an extra variable in the signature, it's bad practice to overwrite the base method) and add an extra variable y_scale in the signature. I then also changed the two lines that have a 0 to 1/y_scale:

    @staticmethod
    def piecewise_linear_y_scale(t, deltas, k, m, changepoint_ts, y_scale):
        """Evaluate the piecewise linear function, keeping the trend
        positive.

        Parameters
        ----------
        t: np.array of times on which the function is evaluated.
        deltas: np.array of rate changes at each changepoint.
        k: Float initial rate.
        m: Float initial offset.
        changepoint_ts: np.array of changepoint times.

        Returns
        -------
        Vector trend(t).
        """
        # Intercept changes
        gammas = -changepoint_ts * deltas
        # Get cumulative slope and intercept at each t
        k_t = k * np.ones_like(t)
        m_t = m * np.ones_like(t)
        for s, t_s in enumerate(changepoint_ts):
            indx = t >= t_s
            k_t[indx] += deltas[s]
            m_t[indx] += gammas[s]
        trend = k_t * t + m_t
        if max(t) <= 1:
            return trend
        # Add additional deltas to force future trend to be positive
        indx_future = np.argmax(t >= 1)
        while min(trend[indx_future:]) < 1/y_scale:
            indx_neg = indx_future + np.argmax(trend[indx_future:] < 1/y_scale)
            k_t[indx_neg:] -= k_t[indx_neg]
            m_t[indx_neg:] -= m_t[indx_neg]
            trend = k_t * t + m_t
        return trend

I'm then overriding the predict_trend() method from the Prophet class and changing the piecewise_linear(...) function call to piecewise_linear_y_scale(..., y_scale):

    def predict_trend(self, df):
        """Predict trend using the prophet model.

        Parameters
        ----------
        df: Prediction dataframe.

        Returns
        -------
        Vector with trend on prediction dates.
        """
        k = np.nanmean(self.params['k'])
        m = np.nanmean(self.params['m'])
        deltas = np.nanmean(self.params['delta'], axis=0)

        t = np.array(df['t'])
        if self.growth == 'linear':
            trend = self.piecewise_linear_y_scale(t, deltas, k, m,
                                          self.changepoints_t, self.y_scale)
        elif self.growth == 'logistic':
            cap = df['cap_scaled']
            trend = self.piecewise_logistic(
                t, cap, deltas, k, m, self.changepoints_t)
        elif self.growth == 'flat':
            # constant trend
            trend = self.flat_trend(t, m)

        return trend * self.y_scale + df['floor']

The issue however is that the min(trend[indx_future:]) in the piecewise_linear_y_scale method keeps being 0, causing an infinite loop.

• The reason for it being 0 seems to be because the k_t and 'm_t' variables are 0 for all values after indx_neg. And of course, when doing k_t * t + m_t, you get 0 if k_t and m_t are 0.
• Then again, the reason for the k_t variable to be 0 for all values after indx_neg seems to be because all values of k_t are the same before it enters the loop (eg. -0.22501 in my example case). So this means that k_t[indx_neg:] -= k_t[indx_neg] will result in 0 for all values after indx_neg. The same applies to m_t.

Do you have any suggestions/pointers on what to fix/where to look, since I'm not familiar with all the different variables? A quick and dirty solution that actually works is to change at the end of the bottom loop (right after trend = k_t * t + m_t) all 0 values to 1/y_scale, but this feels very 'hacky' and probably doesn't work for all cases: trend[trend == 0] = 1/y_scale.

Thanks in advance! 😄

[bletham]

Ah, yes, I should have noticed that sooner. The

            k_t[indx_neg:] -= k_t[indx_neg]
            m_t[indx_neg:] -= m_t[indx_neg]

is resetting the trend to have 0 slope (k_t) and 0 offset (m_t) at indx_neg so that the trend is flat at 0. To have the trend be flat at 1/y_scale, What we really want is to have 0 slope and offset 1/y_scale. So

            k_t[indx_neg:] -= k_t[indx_neg]
            m_t[indx_neg:] = m_t[indx_neg:] - m_t[indx_neg] + 1/y_scale

should do the job. I would worry a little bit about relying on < for the float comparison with 1/y_scale, so I'd probably add a little tolerance in there, like add an extra 1e-6 or something to m_t.