Say we want to do OLS on $m$ samples $(x_i,y_i)$, where $(x_1,y_1), \ldots, (x_{m-1},y_{m-1})$ have the relationship $y_i = \beta \cdot x_i$, but $(x_m,y_m)$ does not. Specifically, $(x_m,y_m)$ is an outlier where $||y_m - \beta \cdot x_i|| = c$.
Can we quantify the effect that this outlier will have on our estimation of $\hat \beta$ on this data? What about when there are multiple such outliers?
Take a look at influence functions. Influence functions measure the effect of a point vector on an estimate of a parameter; in your case, the pair (x,y). It was Originally developed by Hampel in his thesis on robustness.
Mallows and Gnanadesikan showed how to use it in outlier detection. In the x-y plane, you can construct contours of constant influence and thereby determine the direction to look for outliers in a multivariate sense (bivariate in your case). Contours of constant influence are hyperbolae for the bivariate correlation. Because the estimate of the slope parameter in simple linear regression is directly proportional to the bivariate correlation the two influence functions are closely related.
I was motivated by the results in Gnandesikan's multivariate book to apply the influence function for correlation to detect outliers in data validation studies. You can see my paper in the American Journal of Management and Mathematical Science (1983). It provides a "user-friendly" account of this approach. Another reference is Gnanadesikan's book published by Wiley. This is the second edition. The first edition was published in 1977.
See also my paper from Taylor and Francis (current publisher of the journal).
My ORNL technical report from 1979 which the paper is based on can be found for free as a pdf file online.
dfbeta is automatically in R as part of the stats package. This deletion diagnostic gives the relative change in magnitude of the regression coefficient from deleting the $i$-th observation for all $n$ observations in an lm or glm model. It's also closely connected with the jackknife robust error estimator.
$\endgroup$
Using the model form $\boldsymbol{Y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}$, the ordinary least square (OLS) estimator can be written as:
$$\begin{equation} \begin{aligned} \hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\text{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\text{T} \boldsymbol{Y} \\[6pt] &= (\boldsymbol{X}^\text{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\text{T} (\boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}) \\[6pt] &= \boldsymbol{\beta} + (\boldsymbol{X}^\text{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\text{T} \boldsymbol{\varepsilon}. \\[6pt] &= \boldsymbol{\beta} + \sum_{i=1}^n \Big[ (\boldsymbol{X}^\text{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\text{T} \Big] \boldsymbol{\varepsilon}. \\[6pt] \end{aligned} \end{equation}$$
If we define the matrix $\boldsymbol{M} \equiv (\boldsymbol{X}^\text{T} \boldsymbol{X})^{-1} \boldsymbol{X}^\text{T}$ (called the pseudo-inverse of $\boldsymbol{X}$) then we can write the OLS estimator as:
$$\hat{\boldsymbol{\beta}} = \boldsymbol{\beta} + \begin{bmatrix} m_{1,1} & m_{1,2} & \cdots & m_{1,n} \\ m_{2,1} & m_{2,2} & \cdots & m_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ m_{p,1} & m_{p,2} & \cdots & m_{p,n} \\ \end{bmatrix} \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix} = \begin{bmatrix} \beta_1 + \sum m_{1,i} \varepsilon_i \\ \beta_2 + \sum m_{2,i} \varepsilon_i \\ \vdots \\ \beta_p + \sum m_{p,i} \varepsilon_i \end{bmatrix}.$$
In this form you can see that each elements of the OLS estimator includes a linear combination of the error terms in the model, so an outlier for a single value will affect all of the estimated coefficients (unless the relevant weighting in the above expression is zero).
Illustrating a model with a single outlier: To illustrate the situation described in your question, suppose we have standard error terms $\epsilon_1,...,\epsilon_n \sim \text{IID N}(0, \sigma^2)$ and an additional large deviation $\xi \sim \text{N}(0, \kappa^2)$ for the last observation. We define the error terms for the model as:
$$\varepsilon_i = \begin{cases} \epsilon_i & & \text{for } i = 1,...,n-1, \\ \epsilon_i + \xi & & \text{for } i = n. \end{cases}$$
In this case we can write the OLS estimator as:
$$\hat{\boldsymbol{\beta}} = \begin{bmatrix} \beta_1 + \sum m_{1,i} \epsilon_i + m_{1,n} \cdot \xi \\ \beta_2 + \sum m_{2,i} \epsilon_i + m_{2,n} \cdot \xi \\ \vdots \\ \beta_p + \sum m_{p,i} \epsilon_i + m_{p,n} \cdot \xi \end{bmatrix}.$$
Is there a version of the correlation coefficient that is less-sensitive to outliers? discusses estimating the effect of an outlier on the correlation coefficient . This leads naturally to possibly converting the robust correlation coefficient to a robust regression coefficient which is what I believe you are asking for.
