How to calculate uncertainty in the sum of the predictions I am using dummy data to ask my question:
I have a linear model y~x.
library(tidyverse)
library(ggpubr)
theme_set(theme_pubr())
model <- lm(sales ~ youtube, data = marketing)

Let's assume I am predicting into a new dataset:
pred=predict(model,marketing[1:20,], interval="prediction")

I want to calculate the sum of all sales of predicted values (sum(pred$fit)) and also calculate the uncertainty in the estimate. I am a bit unsure how to calculate uncertainty
Am I correct to calculate this uncertainty by: lwr.sum=(sum(lwr^2))^0.5 and upr.sum=(sum(upr^2))^0.5
UPDATE: I tried to calculate using my method, but my mean is not inside the lwr.sum and upr.sum
 A: The predicted values will be strongly correlated because they are all based on the same two estimated coefficients (slope and intercept).
To reason through this, we need some notation.  Let $x_i$ be an array of $n$ new values of the explanatory variable youtube; let $\hat\beta_0$ and $\hat\beta_1$ be the coefficient estimates, so that the predicted value for $x_i$ is
$$\hat y_i  = \hat\beta_0+\hat\beta_1 x_i,$$
and let $\Sigma$ be their variance-covariance matrix (obtainable in R with the function vcov(model)).
The sum of the predictions therefore is
$$\hat p = \sum_{i=1}^n \hat y_i = \sum_{i=1}^n \hat\beta_0 + \hat\beta_1 x_i = n\hat\beta_0 + S_x\hat\beta_1$$
where I have written $S_x$ for the sum of the $x_i.$  The variance is
$$\operatorname{Var}(\hat p) = \operatorname{Var}(n\hat\beta_0) + \operatorname{Var}(S_x\hat\beta_1) + 2\operatorname{Cov}(n\hat\beta_0, S_x\hat\beta_1) = \pmatrix{n&S_x}\Sigma\pmatrix{n\\S_x}.$$
(The sum of the variances looks a little like the formulas proposed in the question; but usually the covariance term cannot be neglected: it will be zero only when the mean of the original data used for the estimates is zero.)
(I wrote this in matrix notation to indicate the generalization to multiple explanatory variables if you are doing multiple regression.)
The square root of $\operatorname{Var}(\hat p)$ is the standard error of this prediction.  Use it in the usual way to erect confidence or prediction limits around that value.
This R code shows the steps in detail: see the three lines at the end where p.hat and se are computed.  The is written to work without change for multiple regression.
#
# Create some data.
#
X <- data.frame(x = seq(1, 19, by = 2)) # Explanatory values
beta <- c(5, 1/2)                       # True coefficients
sigma <- 2                              # Error s.d.
X <- within(X, y <- beta[1] + beta[2] * x + rnorm(nrow(X), 0, sigma))
#
# Fit an OLS model.
#
mod <- lm(y ~ x, X)
#
# Specify the new explanatory variables.
#
X.New <- data.frame(x = seq(3, 8))
#
# Predict the results.
#
X.New$y <- predict(mod, X.New)
#
# Compute the sum of the predictions and their standard error.
#
p.hat <- sum(X.New$y)
v <- colSums(model.matrix(formula(mod), X.New))
se <- sqrt(v %*% vcov(mod) %*% v)
c(`Predicted total` = signif(p.hat, 3),
  `Standard error` = signif(se, 3))

