Combining prediction intervals in regression

Having performed a linear regression, I can find the confidence interval for the response conditioned on a particular x value. However, I am interested in a C.I for the mean response for a set of N new observations. That is, I need to combine the N prediction intervals.

The closest post I could find was Calculating the mean using regression data, but it only handles the univariate case.

I tried deriving the standard error of the mean response below, but I'm not sure if this correct.

\begin{align} var(\hat{\bar{y}}) &= var \left( \frac{1}{n} \sum_i \hat{y}_i|x_1 \ldots x_n \right) \\ &= var \left( \frac{1}{n} \sum_i \hat{y}_i|x_i \right), \quad \text{where the } \hat{y_i}|x_i \text{ are independent} \\ &= \frac{1}{n^2} \sum_i var(\hat{y}_i|x_i) \\ \end{align}

where $var(\hat{y}_i|x_i) = \sqrt{\sigma^2 x_i^T (X^TX)^{-1}x_i}$ for $x_i$ in the training data and $var(\hat{y}_i|x^*_i) = \sqrt{\sigma^2 (1+ x_i^{*T} (X^TX)^{-1}x^*_i)}$ for $x^*_i$ in the test data.

Am I on the right track here? Also, is there an R implementation somewhere, or should I do it from scratch?

Edit: I am also reading up on Bayesian regression methods which specify the predictive distribution $P(y_i|x_i^*)$, and a credible interval for the response. We face a similar problem here, namely, how to compute the predictive distribution for the mean response $P(\overline{y}|x_1^* \ldots x_n^*)$?