Regression slope question Suppose there are $n$ feature vectors $X_1, ...X_n$ which predict $Y$. Suppose we linear regress $X_1$ on $Y$, $X_2$ on $Y$ and so on and obtain all positive slopes $\beta_1>0, \beta_2>0,..\beta_n >0$. Now we perform multiple linear regression $Y = \sum_i\beta_iX_i$, and estimate slopes $\beta_i$, is it necessary that all $\beta_i>0$?
I am looking for a rigorous proof of why it is not necessary $\beta_i>0$ for all $i$.
 A: The answer is no: the univariate regressions can all have positive coefficients, yet the multiple regression can involve negative coefficients.
You ask for intuition.  One route to that is through the geometric interpretation of the Ordinary Least Squares solution, as expressed by the Normal equations: they state that the coefficients in the regression of a column of $n$ data $y$ on explanatory columns $x_i$ are the coefficients of the orthogonal projection of the vector $y$ onto the span of the vectors $x_i$ in $n$-dimensional Euclidean space $R^n.$
Although $n$ dimensions sounds like a lot, you don't need more than two dimensions to visualize useful examples: just limit yourself to $2$ explanatory variables $x_1$ and $x_2.$ They determine a plane (at least when they're not collinear) and we may suppose $y$ is located in that very plane.
Here, to illustrate, is a situation where $y= x_1 - x_2.$  As indicated at the bottom right, although these vectors "live" in $R^n,$ we can draw them faithfully within the plane of the paper.

The gray area at the left is the locus of points whose projection onto some of the $x_i$ has a negative coefficient.  It's the union of all the half-planes behind the $x_i.$  Provided $y$ lies outside that gray area, its regression (i.e., projection) onto any of the variables $x_i$ individually is guaranteed to have a positive value.  Nevertheless, whenever $y$ is outside the infinite wedge lying between the $x_1$ and $x_2$ direction (their positive cone $\{\alpha x_1 + \beta x_2\mid \alpha\ge 0, \beta\ge 0 \}$), it (obviously) cannot have all its coefficients positive. That leaves great scope for counterexamples.
A: No. Let's do a simulation to develop a counterexample.
set.seed(2021)
library(MASS)
N <- 100
X <- MASS::mvrnorm(N, c(0, 0), matrix(c(1, 0.95, 0.95, 1), 2, 2))
x1 <- X[, 1]
x2 <- X[, 2]
y <- x1 - x2 + rnorm(N)
L1 <- lm(y ~ x1)
L2 <- lm(y ~ x2)
L <- lm(y ~ x1 + x2)

This happens if you tack on + 0 in the lm calls to indicate that the intercept should be excluded, too.
set.seed(2021)
library(MASS)
N <- 100
X <- MASS::mvrnorm(N, c(0, 0), matrix(c(1, 0.95, 0.95, 1), 2, 2))
x1 <- X[, 1]
x2 <- X[, 2]
y <- x1 - x2 + rnorm(N)
L1 <- lm(y ~ x1 + 0)
L2 <- lm(y ~ x2 + 0)
L <- lm(y ~ x1 + x2 + 0)

However, the usual advice is to include an intercept unless you have a good reason not to.
