2
$\begingroup$

I am using the statsmodels python package to perform multivariate linear regression. I want to produce 80% prediction interval bands as part of my result.

The statsmodels package can produce prediction intervals for a given alpha and new predictor(s). Fortunately my residuals are normally distributed so the conventional prediction interval for normally distributed residuals is valid. My understanding is that statsmodels uses the conventional prediction interval calculation:https://online.stat.psu.edu/stat501/lesson/3/3.3

I have attached a plot of a sample of my actual target variable / predicted equivalent. The shaded blue region represents my 80% prediction interval (I have bounded it to zero, because negative values are not possible). The model was trained on data of exactly the same shape i.e., Null, rising to some peak, and then returning to 0 again.

My question is why doesn't the prediction interval vary as much as I'd expect it to along the x-axis? I'd expect it to drop to zero in the first and final third of the x-axis, given that this happens for all of the sample data I used to train the model.

It is worth mentioning that the prediction interval I have is not constant - it is dependent on the predictor values.

I know I can use bootstrapping to generate a more accurate prediction interval, but I'm curious why this conventional statistical PI isn't working for my use case.

enter image description here

EDIT: I cannot post my training data, but here is an example code excerpt to show the method I was using (ordinary least-squares):

    from statsmodels.regression.linear_model import OLS, OLSResults

    """
    X_train is a 5-feature (5 column) training set
    y_train is a single-column of target variables

    X_test is a 5-feature testing set
    y_test is a single-column of un-seen target variables

    prediction_summary_frame is the dataframe of predictions produced by statsmodels. It
    has a row of predictions per sample in X_new, with columns: 
    
    mean
    mean_se
    obs_ci_lower
    obs_ci_upper
    mean_ci_lower
    mean_ci_upper
    
    Where mean_ci_upper and mean_ci_lower values represent  the confidence interval, and
    obs_ci_lower and obs_ci_upper represent the prediction intervals. 
    """

    model = OLS(X_train, y_train)
    ols_results = OLS.fit()

    prediction_summary_frame = ols_results.get_prediction(X_test).summary_frame(
        alpha=0.1
    )

     ax.plot(
            X_test.index,
            prediction_summary_frame["mean"],
            label="predicted",
        )
     ax.fill_between(
            X_test.index,
            prediction_summary_frame["obs_ci_upper"],
            prediction_summary_frame["obs_ci_lower"],
            alpha=0.4,
        )
     ax.plot(
            X_test.index,
            y_test,
            label="actual",
        )
     plt.show()

```
$\endgroup$
10
  • $\begingroup$ Can you edit your post to include a Minimal Working Example? $\endgroup$ Commented Aug 18, 2022 at 12:12
  • $\begingroup$ do you mean the code I am using to fit the model and produce predictions? $\endgroup$
    – Archie
    Commented Aug 18, 2022 at 12:26
  • 1
    $\begingroup$ The code and the dataset. Something that will produce the plot you are showing if we run it in a clean interpreter ("Working"), but please reduced as much as possible to the essentials ("Minimal"). $\endgroup$ Commented Aug 18, 2022 at 12:36
  • 1
    $\begingroup$ In a Poisson model the variance and prediction interval depends on the poisson rate. In that case the width of the interval would go to zero as the rate goes to zero. $\endgroup$
    – Josef
    Commented Aug 18, 2022 at 13:41
  • 1
    $\begingroup$ yes it does, but time is not one of the features, because temporal effects are captured in other features I am using $\endgroup$
    – Archie
    Commented Aug 19, 2022 at 6:50

1 Answer 1

2
$\begingroup$

Prediction interval for OLS contains two components, uncertainty about the predicted mean plus uncertainty of a new residual.

In OLS, the assumption is that the residual variance is constant, so the width of the prediction interval coming from the second component will be the same for all values of x.

The model that underlies OLS is linear and does not impose any nonnegativity constraints, so the prediction interval in the first and last part will be centered around mean zero, but with approximately the same width as in other parts because of the homoscedasticity.

A more appropriate model would be Poisson likelihood model or quasi-likelihood model. Poisson is inherently heteroscedastic, the variance is equal to the mean. With a log-link, exponential mean function, the prediction interval would collapse to {0} or {0, 1} set as the Poisson rate or mean goes to zero.

(However, statsmodels currently provides prediction intervals only for OLS, but not for nonlinear models like Poisson and non-gaussian GLM. Prediction intervals that ignore parameter uncertainty can be computed from the properties of the Poisson distribution, but combining those prediction intervals with parameter uncertainty requires more work, e.g. computing tolerance intervals as in https://www.statsmodels.org/dev/generated/statsmodels.stats.rates.tolerance_int_poisson.html )

$\endgroup$
1
  • $\begingroup$ Thank you for taking the time, this an absolutely fantastic answer! $\endgroup$
    – Archie
    Commented Aug 18, 2022 at 22:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.