Adding a variable with no relationship with $Y$ increases $R^2$ I'm confused by something I've found by adding a variable with no relationship with a DV, using multiple regression with four predictors and one DV (Y). If I regress $Y$ onto $X_1,~ X_2,$ and $X_3$, the multiple $R$ is less than if I add a 4th predictor with no relationship with $Y$. I didn't think this was possible. I've done this via a simulation and also more manually, with each shown below. What's even more confusing is that if I specify the 4th variable to have a correlation of $.2$ with the DV, the $R^2$ is less than if the 4th variable has a correlation of $0$ with the DV. How is this possible?
### via simulation ###

library(MASS)
library(psych)

rx12 = .2
rx13 = .25
rx14 = .3
rx23 = .35
rx24 = .3
rx34 = .4

rx1y = .15
rx2y = .25
rx3y = .2
rx4y = 0

corr_matrix  <- matrix(c(1, rx12, rx13, rx14, rx1y, 
                         rx12, 1, rx23, rx24, rx2y,
                         rx13, rx23, 1, rx34, rx3y,
                         rx14, rx24, rx34, 1, rx4y,
                         rx1y, rx2y, rx3y, rx4y, 1), nrow=5)
corr_matrix #this shows the correlation is zero#

set.seed(33)
data = as.data.frame (mvrnorm(n=1000, mu=c(.0, .0, .0, .0, 0), Sigma=corr_matrix, empirical=TRUE)) 
psych::corr.test(data)$r #this shows the correlation is zero#

summary(lm(V5 ~ V1 + V2 + V3, data=data)) #R^2 = .0833
summary(lm(V5 ~ V1 + V2 + V3 + V4, data=data)) #R^2 = .1044

### matrix multiplication with all 4 variables ###
corr_matrix_x  <- matrix(c(1, rx12, rx13, rx14,
                           rx12, 1, rx23, rx24, 
                           rx13, rx23, 1, rx34,
                           rx14, rx24, rx34, 1), nrow=4)
corr_matrix_y  <- matrix(c(rx1y, rx2y, rx3y, rx4y), nrow=4)
corr_matrix_y #this shows the correlation is zero#

x_inverse <- solve(corr_matrix_x)
betas <- as.matrix(x_inverse %*% corr_matrix_y)
t(betas) %*% corr_matrix_y   #R^2 = .1044

### 3 variables ###
corr_matrix_x  <- matrix(c(1, rx12, rx13, 
                           rx12, 1, rx23, 
                           rx13, rx23, 1), nrow=3)
corr_matrix_y  <- matrix(c(rx1y, rx2y, rx3y), nrow=3)

x_inverse <- solve(corr_matrix_x)
betas <- as.matrix(x_inverse %*% corr_matrix_y)
t(betas) %*% corr_matrix_y #R^2 = .0833

Follow-up:
My friend visually examined the issue of R^2 being larger when the X4 correlation with Y is zero or small (~ <.2), and then only increasing at correlations ~>.2. Image below:

 A: This will always be the case until you add enough variables so that the number of independent variables equals the number of observations and the $\mathrm{R}^2 = 1$.  It happens because there's a difference between the theoretical correlation between two random variables and the correlation between samples drawn from their distributions; the samples will, by the nature of randomness, (almost) certainly have nonzero correlations, and will therefore reduce the residual variance slightly, improving $\mathrm{R}^2$.
An example:
y <- rnorm(20)
x <- matrix(rnorm(400),20,20)

summary(lm(y~x-1))
Call:
lm(formula = y ~ x - 1)

Residuals:
ALL 20 residuals are 0: no residual degrees of freedom!

Coefficients:
    Estimate Std. Error t value Pr(>|t|)
x1  -10.7069        NaN     NaN      NaN
x2   11.3354        NaN     NaN      NaN
... and so on ...
x20 -11.5178        NaN     NaN      NaN

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:    NaN 
F-statistic:   NaN on 20 and 0 DF,  p-value: NA

Note the Multiple R-squared at the bottom equals 1!.
A bivariate example may help too:
x <- rnorm(100)
y <- rnorm(100)
cor(x,y)
[1] -0.09046082

A: *

*A variable that has zero correlation with $Y$ can improve the model.

*The correlation with $Y$ does not indicate by how much the model will improve.

*A less strong correlation can be better.


The following example illustrates what can be going on.
Let the true relationship be
$Y = a + b X_1 + \epsilon$
Let the variable that we use in the regression be instead of $X_1$ something slightly changed $Z_1 = X_1 + Z_2$ where $Z_2$ has no relationship with $Y$ but it correlates with $Z_1$.
Now a fit with $Z_1$ will result in a less good $R^2$ than a fit with $X_1$. Adding the variable $Z_2$ to the regression can correct for this by reducing some of the part of $Z_1$ that is not correctly modelling the $Y$ variable. The extra variable $Z_2$ does not need to have a direct relationship with $Y$, it can also work by having a relationship with $Z_1$ and $Y$ combined.

What's even more confusing is that if I specify the 4th variable to have a correlation of $.2$ with the DV, the $R^2$ is less than if the 4th variable has a correlation of $0$ with the DV.

We can visualize this when we consider two explanatory variables $x_1$ and $x_2$ such that we can imagine the vectors in 3D.

The fitted vector $\hat{y}$ will be a vector inside the plane spanned by $x_1$ and $x_2$. The square of the correlation between $\hat{y}$ and $y$ is the $R^2$ value. This value is highest ($R^2 = 1$) when $y$ is inside the plane.
Now consider a given correlation between $y$ and $x_1$ then the vector $y$ will lie on a circle around the vector $x_1$. For different points on that circle the vector $y$ will have different correlations with $x_2$. In the example the vector $y$ has zero correlation with $x_2$ when the vector is inside the plane. This is the worst case example where 0 correlation results in the highest $R^2$. Changing the correlation between $y$ and $x_2$ away from 0 (by choosing another position on the circle) will decrease the $R^2$ value.
A: You should consider the so-called Adjusted R-squared that imposes a penalty for considering additional independent variables to a model. The Un-adjusted  R-squared can never fall when a new independent variable is added to a regression equation: this is because Sum of Squared Residuals (SSR) never goes up (and usually falls) as more independent variables are added (assuming we use the same set of observations). Adjusted R-squared can be calculated as follows:
$$
\overline{R}^2=1-\frac{(1-R^2)(n-1)}{n-k-1}
$$
where $R^2$ is the Un-adjusted R-squared, $k$ is the number of independent variables and $n$ the sample size.
A: It would be so nice if the regression machinery could just pick up the lack of a relationship. Unfortunately, all it does is see potential explainers of the variability, and the algorithm figures out the coefficients that result in the smallest sum of squared errors. That can result in fitting features to the noise (error term). This is related to, if not a light version of, the “overfitting” discussed in machine learning circles.
