Puedes descargar su ejemplo desde https://web.stanford.edu/~hastie/CASI_files/DATA/kidney.txt y fácilmente replicar sus resultados.
> kidney <- read.table("kidney.txt", header=TRUE)
> str(kidney)
'data.frame': 157 obs. of 2 variables:
$ age: int 18 19 19 20 21 21 21 22 22 22 ...
$ tot: num 2.44 3.86 -1.22 2.3 0.98 -0.5 2.74 -0.12 -1.21 0.99 ...
> fit <- lm(tot ~ age, data=kidney)
> fit$coefficients
(Intercept) age
2.86002680 -0.07858842
En cuanto a los errores estándar, los errores estándar de los valores ajustados, $\text{se}(\hat{y})$, son diferentes de los errores estándar de los coeficientes, $\text{se}(\hat\beta)$.
La matriz modelo $X$ es:
> X <- model.matrix(fit)
> head(X)
(Intercept) age
1 1 18
2 1 19
3 1 19
4 1 20
5 1 21
6 1 21
Colocando $S=(X^TX)^{-1}$, $\text{cov}(\hat\beta)=\sigma^2_yS$ (ver mi respuesta a esta pregunta). Dado un único valor ajustado, $\hat{y}_h$ y la fila correspondiente $h$ de $X$, por ejemplo $$y_1=2.44,\qquad x_1=\begin{bmatrix}1 \\ 18\end{bmatrix}$$ la varianza de $\hat{y}_h$ es: $$\text{var}(\hat{y}_h)=\text{var}(x_h^T\hat\beta)=x_h^T\text{cov}(\hat\beta)x_h=x_h^T(S\sigma^2_y)x_h =\sigma^2_y(x_h^TSx_h)$$ Estimas $\sigma^2_y$ por la media cuadrática residual, RMS, el error estándar de $\hat{y}_h$ es: $$\text{se}(\hat{y}_h)=\sqrt{RMS(x_h^TSx_h)}$$ y depende de $x_h$.
Cuando solo hay una variable independiente, $$S=(X^TX)^{-1}=\frac{1}{n\sum(x_i-\bar{x})^2} \begin{bmatrix}\sum x_i^2 & -\sum x_i \\ -\sum x_i & n \end{bmatrix}$$ y \begin{align*} x_h^T(X^TX)^{-1}x_h &=\frac{\sum x_i^2-2x_hn\bar{x}+nx_h^2}{n\sum(x_i-\bar{x})^2}=\frac{\sum x_i^2 -n\bar{x}^2+n(x_h-\bar{x})^2}{n\sum(x_i-\bar{x})^2}\\ &=\frac1n+\frac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2} \end{align*} (Recuerda que $\sum(x_i-\bar{x})^2=\sum x_i^2-n\bar{x}^2$).
La "versión extendida de la fórmula (1.2)" (que es simplemente el error estándar de una media) es: $$\text{se}(\hat{y}_h)=\left[RMS\left(\frac1n+\frac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)\right]^{\frac12}$$ Por cierto, así es como se calculan los intervalos de confianza.
Ver Kutner, Nachtsheim, Neter & Li, Applied Linear Statistical Models, McGraw-Hill, 2005, §2.4, o Seber & Lee, Linear Regression Analysis, John Wiley & Sons, 2003, §6.1.3.
En R:
> S <- solve(t(X) %*% X)
> RMS <- summary(fit)$sigma^2
> x_h <- matrix(c(1, 20), ncol=1) # primer error estándar en la Tabla 1.1
> y_h_se <- sqrt(RMS * (t(x_h) %*% S %*% x_h)); y_h_se
[,1]
[1,] 0.2066481
> x_h <- matrix(c(1, 80), ncol=1) # último error estándar en la Tabla 1.1
> y_h_se <- sqrt(RMS * (t(x_h) %*% S %*% x_h)); y_h_se
[,1]
[1,] 0.420226
EDITAR
Si estás interesado en el error estándar de $\hat{y}_{h(new)}=\hat\alpha+\hat\beta x_{h(new)}$, cuando $x_{h(new)}$ es una nueva observación, no sabes cuál sería $\hat{y}_h$ en una regresión en $n+1$ puntos. Diferentes muestras darían diferentes predicciones, por lo que debes tener en cuenta la desviación de $\hat{y}_{h(new)}$ alrededor de $\hat{y}_h=\hat\alpha+\hat\beta x_h$: $$\text{var}[y_{h(new)}-\hat{y}_h]=\text{var}(y_{h(new)})+\text{var}(\hat{y}_h)$$ Entonces, la varianza de tu predicción tiene dos componentes: la varianza de $y$, que estimas por RMS, y la varianza de la distribución muestral de $\hat{y}_h$, $RMS(x_h^TSx_h)$:
$$RMS + RMS\left(\frac1n+\frac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)$$ La "versión extendida de la fórmula (1.2)" se convierte en: $$\text{se}(\hat{y}_{h(new)})=\left[RMS\left(1+\frac1n+\frac{(x_{h(new)}-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)\right]^{\frac12}$$ Ver Kutner, Nachtsheim, Neter & Li, Applied Linear Statistical Models, McGraw-Hill, 2005, §2.5, o https://online.stat.psu.edu/stat501/lesson/3/3.3.
