{"id":239,"date":"2021-03-14T11:20:54","date_gmt":"2021-03-14T10:20:54","guid":{"rendered":"https:\/\/lorentzen.ch\/?p=239"},"modified":"2021-03-14T11:20:56","modified_gmt":"2021-03-14T10:20:56","slug":"a-beautiful-regression-formula","status":"publish","type":"post","link":"https:\/\/lorentzen.ch\/index.php\/2021\/03\/14\/a-beautiful-regression-formula\/","title":{"rendered":"A Beautiful Regression Formula"},"content":{"rendered":"\n

Lost in Translation between R and Python 4<\/h1>\n\n\n\n

Hello statistics aficionados<\/p>\n\n\n\n

This is the next article in our series “Lost in Translation between R and Python”<\/strong>. The aim of this series is to provide high-quality R and<\/strong> Python 3 code to achieve some non-trivial tasks. If you are to learn R, check out the R tab below. Similarly, if you are to learn Python, the Python tab will be your friend.<\/p>\n\n\n\n

The last one was a deep dive into historic mortality<\/a> rates.<\/p>\n\n\n\n

No Covid-19, no public data for a change: This post focusses on a real beauty, namely a decomposition of the R-squared in a linear regression model<\/p>\n\n\n\n

E(y) = \\alpha + \\sum_{j = 1}^p x_j \\beta_j<\/pre><\/div>\n\n\n\n
fitted by least-squares. If the response y and all p covariables are standardized to variance 1<\/strong> beforehand, then the R-squared<\/a> can be obtained as the cross-product of the fitted coefficients and the usual correlations between each covariable and the response:<\/p>\n\n\n\n
R^2 = \\sum_{j = 1}^p \\text{cor}(y, x_j)\\hat\\beta_j.<\/pre><\/div>\n\n\n\n
Two elegant derivations can be found in this answer<\/a> to the same question, written by the number 1 contributor to crossvalidated<\/a>: whuber. Look up a couple of his posts – and statistics will suddenly feel super easy and clear.<\/p>\n\n\n\n
Direct consequences of the formula are:<\/p>\n\n\n\n
  1. If a covariable is uncorrelated with the response, it cannot contribute to the R-squared, i.e. neither improve nor worsen. This is not obvious.<\/li>
  2. A correlated covariable only improves R-squared if its coefficient is non-zero. Put differently: if the effect of a covariable is already fully covered by the other covariables, it does not improve the R-squared. This is<\/em> somewhat obvious.<\/li><\/ol>\n\n\n\n
    Note that all formulas refer to in-sample calculations.<\/p>\n\n\n\n
    Since we do not want to bore you with math, we simply demonstrate the result with short R and Python codes based on the famous iris dataset. <\/p>\n\n\n
    \n\t\t\t
    \n\t\t\t\t
    \n\t\t\t\t
    R<\/div>\n\t\t\t<\/div>
    \n\t\t\t\t
    Python<\/div>\n\t\t\t<\/div><\/div>\n\t\t\t<\/div>\n\t\t\t
    \n\n
    y <- \"Sepal.Width\"\nx <- c(\"Sepal.Length\", \"Petal.Length\", \"Petal.Width\")\n\n# Scaled version of iris\niris2 <- data.frame(scale(iris[c(y, x)]))\n\n# Fit model \nfit <- lm(reformulate(x, y), data = iris2)\nsummary(fit) # multiple R-squared: 0.524\n(betas <- coef(fit)[x])\n# Sepal.Length Petal.Length  Petal.Width \n#    1.1533143   -2.3734841    0.9758767 \n\n# Correlations (scaling does not matter here)\n(cors <- cor(iris[, y], iris[x]))\n# Sepal.Length Petal.Length Petal.Width\n#   -0.1175698   -0.4284401  -0.3661259\n \n# The R-squared?\nsum(betas * cors) # 0.524<\/pre><\/div>\n\n<\/div>
    \n\n
    # Import packages\nimport numpy as np\nimport pandas as pd\nfrom sklearn import datasets\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\n\n# Load data\niris = datasets.load_iris(as_frame=True).data\nprint(\"The data:\", iris.head(3), sep = \"\\n\")\n\n# Specify response\nyvar = \"sepal width (cm)\"\n\n# Correlations of everyone with response\ncors = iris.corrwith(iris[yvar]).drop(yvar)\nprint(\"\\nCorrelations:\", cors, sep = \"\\n\")\n\n# Prepare scaled response and covariables\nX = StandardScaler().fit_transform(iris.drop(yvar, axis=1))\ny = StandardScaler().fit_transform(iris[[yvar]])\n\n# Fit linear regression\nOLS = LinearRegression().fit(X, y)\nbetas = OLS.coef_[0]\nprint(\"\\nScaled coefs:\", betas, sep = \"\\n\")\n\n# R-squared via scikit-learn: 0.524\nprint(f\"\\nUsual R-squared:\\t {OLS.score(X, y): .3f}\")\n\n# R-squared via decomposition: 0.524\nrsquared = betas @ cors.values\nprint(f\"Applying the formula:\\t {rsquared: .3f}\")<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n
    Indeed: the cross-product of coefficients and correlations equals the R-squared of 52%.<\/strong><\/p>\n\n\n\n
    The Python notebook and R code can be found at:<\/p>\n\n\n\n