This answer responds only to a small part of the OP's post, where he posits that outliers can be detected when the absolute studentised residuals are greater than three.
Diagnosing 'outliers' using studentised residuals: I note that you are diagnosing 'outliers' as points that have an absolute studentised residual greater than three. The hyperlink you provide to James, Witten, Hastie and Tibshirani (2017) notes these points as 'possible outliers' (p. 97). They refer to removing data points only in the case of actual measurement error, and also note that an outlier many indicate a deficiency in the model. This 'rule-of-thumb' is also commonly taught in introductory courses, and is in a lot of textbooks. It is extremely unfortunate that this rule is still being taught; it is an idea that should have died more than fifty years ago.
This rule-of-thumb for outliers is completely wrong - it has no basis in statistical theory.
Distribution of maximum absolute studentised residual: When you look at extreme values based on absolute studentised residuals, you must bear in mind that you are cherry-picking a maximum value from a set of underlying random variables; it is not appropriate to compare this value to the tails of the distribution for a single underlying random variable. You need to consider the distribution of the statistic you are actually looking at.
Under the standard Gaussian regression assumption that $\varepsilon_1, ..., \varepsilon_n \sim \text{IID N}(0, \sigma^2)$ your (externalised) studentised residuals are only weakly dependent and have distribution:
$$S_1, ..., S_n \sim \text{Student T}(df_{Res}).$$
For large $n$, the distribution of the maximum absolute studentised residual is well-approximated by an extreme value distribution which depends on the value $n$:
$$M_n \equiv \max_{i = 1,...,n} |S_i| \sim \text{EV}(n),$$
As $n$ gets bigger, this distribution shifts to the right and its expected value becomes bigger, without any upper bound. It has approximate mean $\mathbb{E}(M_n) \approx \Phi^{-1}(1 - 1/n)$, which tends to infinity as $n \rightarrow \infty$. This is hardly surprising: as you take more and more approximately IID random variables from an unbounded distribution, the (absolute) biggest one tends to get bigger and bigger, without any upper bound. So, the idea that you can diagnose an 'outlier' by observing that $M>3$, without any consideration of the sample size, is just absurd.$\dagger$
There are various formal tests for outliers that take account of sample size and thereby correctly use the fact that the maximum studentised residual has a distribution that increases with $n$. The most well-known test is Grubbs' test, so that is a good place to start. Proper outlier tests operate by comparing the maximum value to an appropriate extreme value distribution that approximates the distribution of that maximum value, under the model assumptions. If you find an outlier that is diagnosed with Grubbs' test then this tells you that this extreme value is so extreme that it is unlikely to have come from an underlying Gaussian regression model, even after taking account of the sample size.
What does it actually mean to find one or more 'outliers'? Assuming you have tested correctly, using a test that accounts for sample size, you may find that you are able to identify one or more points whose values are so extreme that they are 'outliers' with respect to the model assumption of normally distributed error terms. So, what does this actually mean?
Unless you have reason to believe that there are actual measurement errors in your data, all this really means is that the underlying errors in your model have a distribution with fatter tails than the normal distribution. The normal distribution has thin tails, and often data does not fit this well. If your data fails Grubbs' test, or some other appropriate outlier test, you might consider changing your model to one with an error distribution that has fatter tails. A simple but effective change is to use a GLM with a generalised error distribution, which adds an additional parameter that allows for additional kurtosis.
It is a very bad idea to remove data points purely because they are 'outliers' in comparison to some assumed error distribution. If you do this, you are effectively requiring the data to meet your model assumptions rather than requiring your model assumptions to conform to the reality of the data. Filtering out data points that have high absolute studentised residuals means that you will systematically underestimate the variability in your data. Unless there is reason to believe that there is actual measurement error in the data point, you should keep it in your data. If you have outliers from this then you should adapt your model to allow for higher kurtosis in your error distribution.
$\dagger$ Again, to be clear, no derogatory comment about this rule should be interpreted as a derogatory comment about the OP. This is a rule that keeps getting taught in statistics courses and is also found in many (otherwise excellent) statistics textbooks. It is rarely corrected when it comes up, and a concerted effort by the statistics profession is needed to kill it off.
