Having issues interpreting my prediction results for my negative binomial regression model

A little background, here is my model:

mod1 = glm.nb(Goals~Defense, data = Messi.Liga)
summary(mod1)


And here is the output:

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.05943 1.11817 2.736 0.00622 **
Defense -0.03954 0.01498 -2.639 0.00831 **


Just so you get an idea of what the data frames look like, goals are the number of goals scored in a single game and Defense is the Defensive Rating of the opposition.

Prediction model:

BPLpredictions = data.frame(predict(mod1, bpl.df, type = "response", se.fit = TRUE))


Everything seems to be going smoothly but i've seen blog posts where they will take the exponential of the "fit" results and use that for interpretation but when I do that, I get completely unreasonable results. What's odd is my results prior to exponentiating "fit" are logical.

Also, to create my confidence interval, I've been using the following code:

BPLpredictions["upr"] = BPLpredictions$fit + (z * BPLpredictions$se.fit)
BPLpredictions["lwr"] = BPLpredictions$fit - (z * BPLpredictions$se.fit)


If you think I should be doing it differently, please feel free to give recommendations.

Thanks.

The question whether or not to exponentiate, is a question of what scale is employed for the normal approximation for the confidence intervals.

In negative binomial regressions you usually use a log-link, i.e., $\log(\mu) = \eta = x^\top \beta$, where $\mu$ is the expectation of the response and $\eta = x^\top \beta$ is the corresponding linear predictor. Given fitted regression coefficients $\hat \beta$ one can easily obtain the fitted linear predictor $\hat \eta = x^\top \hat \beta$ (called type = "link" for glm.nb()) and the corresponding expected response $\hat \mu = \exp(\hat \eta)$ (called type = "response" for glm.nb()).

However, for obtaining the corresponding confidence intervals it makes a difference whether the normal approximation is applied for $\hat \eta$ or for $\hat \mu$. Asymptotically, both versions agree but in finite samples they are different. The predict() method for glm.nb() can return standard errors on either scale but some might argue that it is more natural to apply the normal approximation on the linear predictor scale. The reason is simply that the linear predictor lives on the real line (like the normal distribution) and subsequently applying the inverse link $\exp(\cdot)$ assures positivity of the interval. So I would guess that the blog post you are referring to took this approach.

In practice, however, the difference may be rather small. As a simple reproducible example (that is a really non-sensical but requires no special packages) consider:

library("MASS")
m <- glm.nb(dist ~ speed, data = cars)


So we model how stopping distance increases with speed. (It's not really count data...but it's easy to visualize and does not require any further R packages here.) Then we can add the fitted values and standard errors on the response ($\mu$) scale and on the link ($\eta$) scale:

cars <- cbind(cars, "resp" = predict(m, type = "response", se.fit = TRUE)[1:2])
cars <- cbind(cars, "link" = predict(m, type = "link", se.fit = TRUE)[1:2])


Plotting the data, the fitted expectation, and the confidence intervals from the response scale (in red) yields the figure on the left below.

This is obtained with the usual approximation adding 2 times the standard errors:

plot(dist ~ speed, data = cars)
lines(resp.fit ~ speed, data = cars, col = 2)
lines(resp.fit + 2 * resp.se.fit ~ speed, data = cars, col = 2)
lines(resp.fit - 2 * resp.se.fit ~ speed, data = cars, col = 2)


Adding also the exponentiated fitted values and confidence interval from the link scale (in blue) then yields the plot on the right.

lines(exp(link.fit) ~ speed, data = cars, col = 4)
lines(exp(link.fit + 2 * link.se.fit) ~ speed, data = cars, col = 4)
lines(exp(link.fit - 2 * link.se.fit) ~ speed, data = cars, col = 4)


You can see that the fitted responses exactly coincide but the confidence intervals are slightly different, especially for higher values. Of course, the differences may or may not be larger on other data - but generally they will decrease with the sample size.