Why can't we add all the individual Pearson's $r$'s in a multiple regression and calculate $R^2$ based on this sum? Is there an easy mathematical explanation to this as $r^2$ is squared and don't add up some way?
It is hard to know how to answer this. I guess I would just say that you can't add up the univariate Pearson's $r$'s and get $R$ or $R^2$. Bear in mind that $r$ is bound by $[-1,\ 1]$, and that $R^2$ is bound by $[0,\ 1]$. So it may be easy to see that simple addition would yield sums that violate those rules. Here is an example you can play with yourself:
## here are 3 simple variables x1, x2, & y: x1 = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20) x2 = c(2,1,4,3,6,5,8,7,10,9,12,11,14,13,16,15,18,17,20,19) y = c(3,1,2,6,4,5,9,7,8,12,10,11,15,13,14,18,16,17,20,19) ## here are their Pearson's correlations: cor(x1,y) # r = 0.9714286 cor(x2,y) # r = 0.9593985 cor(x1,x2) # r = 0.9849624 ## here I fit a multiple regression model & get R, & R^2: model = lm(y~x1+x2) cor(y, fitted(model)) # R = 0.9715432 summary(model)$r.squared # R^2 = 0.9438961 ## we can square R to get R^2: 0.9715432^2 # = 0.9438962 ## adding r(x1,y) + r(x2,y) does not equal either R or R^2, though: 0.9714286 + 0.9593985 # = 1.930827