Is it possible to create a 95% confidence interval around an abline using the basic plot function/ include a simple slope-intercept formula in ggplot?

Question

I am running a Phylogenetic Generalized Least Squares (PGLS) model, and want to show the correlation when I correct for relatedness as a line on my graph. Because of this, I cannot use the normal geom_smooth() function in ggplot. The normal abline() function works fine to make the line, as does the geom_abline()in ggplot, but neither of them allow me to show the 95% confidence interval around this line. I have tried geom_smooth(), but perhaps I am just missing something about it, but I keep getting errors. Is it possible to have the 95% confidence interval around a line created by abline() or geom_abline()? Or to create a formula in geom_smooth() which is just the intercept and slope? For interest, my intercept is 1.91253 and my slope is 0.0001783464.

I do not know how to share data in such a way that it looks decent on here, but the data is one column for species names, one column for the trait in question (in this case, colour dichromatism as scored by the Vorobyev-Osorio method), and another column for an environmental variable (in this case, average altitude).

Confidence interval is 1.96 times the Standard Error of the PGLS model, which in this case is 0.865730432.

I tried using geom_abline() and geom_ribbon() and got an error. I first tried to save the abline which I wanted to as its own vector doing

ablineCDG <- as.vector(abline(a = coef(CDGreaterOne)[1], b = coef(CDGreaterOne)[2]))

And then ran

ggplot(SigCD, aes(x=Altitude_Reported, y=Colour_discriminability_Absolute)) +
    geom_point()+
    geom_abline(slope = 0.0001783464, intercept = 1.91253)+
    geom_ribbon(aes(ymin = ablineCDG - 0.865730432, ymax = ablineCDG + 0.865730432), fill = "grey70")+
    labs(y= "Vorobyev-Osorio Colour Discrimination Score", x = "Altitude")

Which gave me the following error

Error in `geom_ribbon()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 3rd layer.
Caused by error in `check_aesthetics()`:
! Aesthetics must be either length 1 or the same as the data (131)
✖ Fix the following mappings: `ymin` and `ymax`
Run `rlang::last_trace()` to see where the error occurred.
> rlang::last_trace()
<error/rlang_error>
Error in `geom_ribbon()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 3rd layer.
Caused by error in `check_aesthetics()`:
! Aesthetics must be either length 1 or the same as the data (131)
✖ Fix the following mappings: `ymin` and `ymax`
---
Backtrace:
     ▆
  1. ├─base (local) `<fn>`(x)
  2. └─ggplot2:::print.ggplot(x)
  3.   ├─ggplot2::ggplot_build(x)
  4.   └─ggplot2:::ggplot_build.ggplot(x)
  5.     └─ggplot2:::by_layer(...)
  6.       ├─rlang::try_fetch(...)
  7.       │ ├─base::tryCatch(...)
  8.       │ │ └─base (local) tryCatchList(expr, classes, parentenv, handlers)
  9.       │ │   └─base (local) tryCatchOne(expr, names, parentenv, handlers[[1L]])
 10.       │ │     └─base (local) doTryCatch(return(expr), name, parentenv, handler)
 11.       │ └─base::withCallingHandlers(...)
 12.       └─ggplot2 (local) f(l = layers[[i]], d = data[[i]])
 13.         └─l$compute_aesthetics(d, plot)
 14.           └─ggplot2 (local) compute_aesthetics(..., self = self)
 15.             └─ggplot2:::check_aesthetics(evaled, n)

When I tried to run this using the geom_smooth() function, I tried to use the same formula as I had for the geom_ribbon() example above, namely setting the desired line as a vector and then plus and minus with the 95% confidence interval. It thus reads as:

ggplot(SigCD, aes(x=Altitude_Reported, y=Colour_discriminability_Absolute)) +
    geom_point()+
    geom_smooth(formula = ablineCDG, level = 0.95, color = "black")+ 
    labs(y= "Vorobyev-Osorio Colour Discrimination Score", x = "Altitude")

But I got the following message:

geom_smooth() using method = 'loess' Warning message: Computation failed in stat_smooth()

I hope this includes all of the information which is needed to solve this issue.

Edit 2

Running the vcov() function on the fitted model gives the following output

                   (Intercept) Altitude_Reported
(Intercept)        1.950982e-01     -2.144661e-05
Altitude_Reported -2.144661e-05      2.003417e-08

How would I go about storing the 95% confidence interval for every value x?