Exercise8.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/choiitea/Choi_DSPN_S24/blob/main/Exercise8.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "szkkhiCZDF52"
      },
      "source": [
        "# Exercise 8:  Linear models, continued\n",
        "\n",
        "This homework assignment is designed to give you a deeper understanding of linear models. First, we'll dive into the math behind the closed-form solution of maximum likelihood estimation. **In the first section below, write your answers using Latex equation formatting.**\n",
        "\n",
        "*Note: Check out [this page](https://gtribello.github.io/mathNET/assets/notebook-writing.html) and [this page](https://towardsdatascience.com/write-markdown-latex-in-the-jupyter-notebook-10985edb91fd) for resources on how to do Latex formatting. You can also double click on the question cells in this notebook to see how math is formatted in the questions.*\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QJscNReoylRt"
      },
      "source": [
        "---\n",
        "## 1. Deriving the Maximum Likelihood Estimate for Simple Linear Regression (6 points)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nH82gwuymPi0"
      },
      "source": [
        "Using the mean squared error (MSE) as your objective function (the thing you're trying to minimize when you fit your model) allows for a closed form solution to finding the maximum likelihood estimate (MLE) of your model parameters in linear regression. Let’s consider the simple, single predictor variable model, i.e. simple linear regression :  $Y= \\beta_0 + \\beta_1 X $.\n",
        "\n",
        "a) Use algebra to show how you can expand out $MSE(\\beta_0, \\beta_1)$ to get from i to ii below.\n",
        "\n",
        "> _i)_ $E[ (Y-(\\beta_0 + \\beta_1 X))^2]$\n",
        "\n",
        "> _ii)_ $E[Y^2] -2 \\beta_0E[Y]-2 \\beta_1 Cov[X,Y]-2 \\beta_1 E[X]E[Y]+ \\beta_0^2 +2 \\beta_0 \\beta_1 E[X]+\\beta_1^2 Var[X]+ \\beta_1^2 (E[X])^2$\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dn2hveNho-Of"
      },
      "source": [
        "**Answer:**\n",
        "\n",
        "> $E[ (Y-(\\beta_0 + \\beta_1 X))^2]$\\\n",
        "> $=E[(Y^2-2Y(\\beta_0+\\beta_1X)+(\\beta_0+\\beta_1X)^2)]$\\\n",
        "> $=E[Y^2]+E[-2\\beta_0Y-2\\beta_1XY]+E[(\\beta_0+\\beta_1X)^2]$\\\n",
        "> $=E[Y^2]-2E[\\beta_0Y+\\beta_1XY]+E[(\\beta_0+\\beta_1X)^2]$\\\n",
        "> $=E[Y^2]-2(E[\\beta_0Y]+E[\\beta_1XY])+E[(\\beta_0+\\beta_1X)^2]$\\\n",
        "> $=E[Y^2]-2\\beta_0E[Y]-2\\beta_1E[XY]+E[(\\beta_0)^2+2\\beta_0\\beta_1X+(\\beta_1)^2X^2]$\\\n",
        "> $=E[Y^2]-2\\beta_0E[Y]-2\\beta_1E[XY]+E[(\\beta_0)^2]+E[2\\beta_0\\beta_1X]+E[(\\beta_1)^2X^2]$\\\n",
        "> $=E[Y^2]-2\\beta_0E[Y]-2\\beta_1E[XY]+(\\beta_0)^2+2\\beta_0\\beta_1E[X]+(\\beta_1)^2E[X^2]$\n",
        "> $=E[Y^2]-2\\beta_0E[Y]-2\\beta_1(Cov(x,Y)+E[X]E[Y])+(\\beta_0)^2+2\\beta_0\\beta_1E[X]+(\\beta_1)^2E[X^2]$\n",
        "> $=E[Y^2]-2\\beta_0E[Y]-2\\beta_1(Cov(x,Y)+E[X]E[Y])+(\\beta_0)^2+2\\beta_0\\beta_1E[X]+(\\beta_1)^2(Var(X)+E[X]^2)$\\\n",
        "> $=E[Y^2]-2\\beta_0E[Y]-2\\beta_1Cov[X,Y]-2 \\beta_1 E[X]E[Y]+\\beta_0^2+2 \\beta_0\\beta_1E[X]+\\beta_1^2Var[X]+\\beta_1^2(E[X])^2$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GCr46r9xwRXP"
      },
      "source": [
        "b) Prove that the MLE of $\\beta_0$ is $E[Y]- \\beta_1 E[X]$ by taking the derivative of _ii_ above, with respect to $\\beta_0$, setting the derivative to zero, and solving for $\\beta_0$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ul-PZyLbwTCQ"
      },
      "source": [
        "**Answer:**\n",
        "\n",
        "> To get the MLE of $\\beta_0$, we take the derivative of MSE with respect to $\\beta_0$ and find the solution for $\\beta_0$ where drivative of MSE with respect to $\\beta_0$ equals 0.\\\n",
        "> MSE is $E[Y^2]-2\\beta_0E[Y]-2\\beta_1Cov[X,Y]-2 \\beta_1 E[X]E[Y]+\\beta_0^2+2 \\beta_0\\beta_1E[X]+\\beta_1^2Var[X]+\\beta_1^2(E[X])^2$\\\n",
        "> So,\\\n",
        "> $\\frac{d(MSE)}{d\\beta_0}=\\frac{d}{d\\beta_0}\\left(E[Y^2] -2\\beta_0E[Y]-2 \\beta_1 Cov[X,Y]-2 \\beta_1 E[X]E[Y]+ \\beta_0^2 +2 \\beta_0 \\beta_1 E[X]+\\beta_1^2 Var[X]+ \\beta_1^2 (E[X])^2 \\right)$\\\n",
        "> $=\\frac{d}{d\\beta_0}\\left(-2\\beta_0E[Y]+\\beta_0^2 +2\\beta_0\\beta_1E[X]\\right)$\\\n",
        "> $=-2E[Y]+2\\beta_0+2\\beta_1E[X]$\\\n",
        "> Now, setting this derivative equal to 0 to find the value of $\\beta_0$ that minimizes the MSE,\\\n",
        "> $-2E[Y]+2\\beta_0+2\\beta_1E[X]=0$\\\n",
        "> $2\\beta_0=2E[Y]-2\\beta_1E[X]$\\\n",
        "> $\\beta_0=E[Y]-\\beta_1E[X]$\\\n",
        "> This shows that MLE of $\\beta_0$ is $E[Y]-\\beta_1E[X]$, derived by minimizing the MSE\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-uv4Z7afw4gB"
      },
      "source": [
        "c) Prove that the MLE for $\\beta_1$ is $Cov[X,Y]/Var[X]$ by taking the derivative of equation _ii_ above, with respect to $\\beta_1$, setting the derivative to zero, and solving for $\\beta_1$. *Hint: after you've simplified / expanded a bit, plug in the solution for $\\beta_0$ from part b.*"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sWTFZ6ZSw6sh"
      },
      "source": [
        "**Answer:**\n",
        "\n",
        "> $\\frac{d(MSE)}{d\\beta_1}=\\frac{d}{d\\beta_1}\\left(E[Y^2] -2\\beta_0E[Y]-2 \\beta_1 Cov[X,Y]-2 \\beta_1 E[X]E[Y]+ \\beta_0^2 +2 \\beta_0 \\beta_1 E[X]+\\beta_1^2 Var[X]+ \\beta_1^2 (E[X])^2 \\right)$\\\n",
        "> $=\\frac{d}{d\\beta_1}\\left(-2\\beta_1Cov[X,Y]-2 \\beta_1 E[X]E[Y]+2 \\beta_0 \\beta_1 E[X]+\\beta_1^2 Var[X]+ \\beta_1^2 (E[X])^2 \\right)$\\\n",
        "> $=-2Cov[X,Y]-2E[X]E[Y]+2\\beta_0E[X]+2\\beta_1Var[X]+2\\beta_1(E[X])^2$\\\n",
        "> Previously, we showed that the MLE of $\\beta_0$ is $E[Y]-\\beta_1E[X]$, so we can substitute $\\beta_0$ into the derivative.\\\n",
        "> $=-2Cov[X,Y]-2E[X]E[Y]+2(E[Y]-\\beta_1E[X])E[X]+2\\beta_1Var[X]+2\\beta_1(E[X])^2$\\\n",
        "> $=-2Cov[X,Y]-2E[X]E[Y]+2E[X]E[Y]-2\\beta_1(E[X])^2+2\\beta_1Var[X]+2\\beta_1(E[X]^2)$\\\n",
        "> $=-2Cov[X,Y]+2\\beta_1Var[X]$\\\n",
        "> $=-2Cov[X,Y]+2\\beta_1Var[X] = 0$\\\n",
        "> $2\\beta_1Var[X]=2Cov[X,Y]$\\\n",
        "> $\\beta_1Var[X]=Cov[X,Y]$\\\n",
        "> $\\beta_1=\\frac{Cov[X,Y]}{Var[X]}$\\\n",
        "> This shows that the MLE for $\\beta_1$ is $\\frac{Cov[X,Y]}{Var[X]}$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "66X264ZpDF58"
      },
      "source": [
        "---\n",
        "## 2. Connecting to data (4 points)\n",
        "\n",
        "Now let's connect this to some real data. Once again we'll be using the  **unrestricted_trimmed_1_7_2020_10_50_44.csv** file from the *Homework/hcp_data* folder in the class GitHub repository.\n",
        "\n",
        "​\n",
        "This data is a portion of the [Human Connectome Project database](http://www.humanconnectomeproject.org/). It provides measures of cognitive tasks and brain morphology measurements from 1206 participants. The full description of each variable is provided in the **HCP_S1200_DataDictionary_April_20_2018.csv** file in the *Homework/hcp_data* folder in the class GitHub repository.\n",
        "\n",
        "a) Use the `setwd` and `read.csv` functions to load data from the **unrestricted_trimmed_1_7_2020_10_50_44.csv** file. Then use the `tidyverse` tools make a new dataframe `d1` that only inclues the subject ID (`Subject`), Flanker Task performance (`Flanker_Unadj`), and total grey matter volume (`FS_Total_GM_Vol`) variables and remove all _NA_ values.\n",
        "\n",
        "Use the `head` function to look at the first few rows of each data frame."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "PZ0lngBjDF58",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 465
        },
        "outputId": "3906f095-fc53-4039-d2a9-1f673db3fe6b"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "── \u001b[1mAttaching core tidyverse packages\u001b[22m ──────────────────────── tidyverse 2.0.0 ──\n",
            "\u001b[32m✔\u001b[39m \u001b[34mdplyr    \u001b[39m 1.1.4     \u001b[32m✔\u001b[39m \u001b[34mreadr    \u001b[39m 2.1.5\n",
            "\u001b[32m✔\u001b[39m \u001b[34mforcats  \u001b[39m 1.0.0     \u001b[32m✔\u001b[39m \u001b[34mstringr  \u001b[39m 1.5.1\n",
            "\u001b[32m✔\u001b[39m \u001b[34mggplot2  \u001b[39m 3.4.4     \u001b[32m✔\u001b[39m \u001b[34mtibble   \u001b[39m 3.2.1\n",
            "\u001b[32m✔\u001b[39m \u001b[34mlubridate\u001b[39m 1.9.3     \u001b[32m✔\u001b[39m \u001b[34mtidyr    \u001b[39m 1.3.1\n",
            "\u001b[32m✔\u001b[39m \u001b[34mpurrr    \u001b[39m 1.0.2     \n",
            "── \u001b[1mConflicts\u001b[22m ────────────────────────────────────────── tidyverse_conflicts() ──\n",
            "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mfilter()\u001b[39m masks \u001b[34mstats\u001b[39m::filter()\n",
            "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mlag()\u001b[39m    masks \u001b[34mstats\u001b[39m::lag()\n",
            "\u001b[36mℹ\u001b[39m Use the conflicted package (\u001b[3m\u001b[34m<http://conflicted.r-lib.org/>\u001b[39m\u001b[23m) to force all conflicts to become errors\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<table class=\"dataframe\">\n",
              "<caption>A data.frame: 6 × 3</caption>\n",
              "<thead>\n",
              "\t<tr><th></th><th scope=col>Subject</th><th scope=col>Flanker_Unadj</th><th scope=col>FS_Total_GM_Vol</th></tr>\n",
              "\t<tr><th></th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
              "</thead>\n",
              "<tbody>\n",
              "\t<tr><th scope=row>1</th><td>100206</td><td>130.42</td><td>807245</td></tr>\n",
              "\t<tr><th scope=row>2</th><td>100307</td><td>112.56</td><td>664124</td></tr>\n",
              "\t<tr><th scope=row>3</th><td>100408</td><td>121.18</td><td>726206</td></tr>\n",
              "\t<tr><th scope=row>4</th><td>100610</td><td>126.53</td><td>762308</td></tr>\n",
              "\t<tr><th scope=row>5</th><td>101006</td><td>101.85</td><td>579632</td></tr>\n",
              "\t<tr><th scope=row>6</th><td>101107</td><td>107.04</td><td>665024</td></tr>\n",
              "</tbody>\n",
              "</table>\n"
            ],
            "text/markdown": "\nA data.frame: 6 × 3\n\n| <!--/--> | Subject &lt;int&gt; | Flanker_Unadj &lt;dbl&gt; | FS_Total_GM_Vol &lt;int&gt; |\n|---|---|---|---|\n| 1 | 100206 | 130.42 | 807245 |\n| 2 | 100307 | 112.56 | 664124 |\n| 3 | 100408 | 121.18 | 726206 |\n| 4 | 100610 | 126.53 | 762308 |\n| 5 | 101006 | 101.85 | 579632 |\n| 6 | 101107 | 107.04 | 665024 |\n\n",
            "text/latex": "A data.frame: 6 × 3\n\\begin{tabular}{r|lll}\n  & Subject & Flanker\\_Unadj & FS\\_Total\\_GM\\_Vol\\\\\n  & <int> & <dbl> & <int>\\\\\n\\hline\n\t1 & 100206 & 130.42 & 807245\\\\\n\t2 & 100307 & 112.56 & 664124\\\\\n\t3 & 100408 & 121.18 & 726206\\\\\n\t4 & 100610 & 126.53 & 762308\\\\\n\t5 & 101006 & 101.85 & 579632\\\\\n\t6 & 101107 & 107.04 & 665024\\\\\n\\end{tabular}\n",
            "text/plain": [
              "  Subject Flanker_Unadj FS_Total_GM_Vol\n",
              "1 100206  130.42        807245         \n",
              "2 100307  112.56        664124         \n",
              "3 100408  121.18        726206         \n",
              "4 100610  126.53        762308         \n",
              "5 101006  101.85        579632         \n",
              "6 101107  107.04        665024         "
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "# INSERT CODE HERE\n",
        "library(tidyverse)\n",
        "\n",
        "system(\"gdown --id 1hywRmGdvhbDYTrQRyl1_bLJsq-T3GJq2\")\n",
        "\n",
        "d <- read.csv(\"unrestricted_trimmed_1_7_2020_10_50_44.csv\")\n",
        "\n",
        "d1 <- d %>%\n",
        "  select(Subject, Flanker_Unadj, FS_Total_GM_Vol) %>%\n",
        "  drop_na()\n",
        "\n",
        "head(d1)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J3owDQ0U2Ewn"
      },
      "source": [
        "b) Now we're going to see if the solutions we proved above actually line up with the model fit that R gives us (it should...). Calculate what the $\\beta_0$ and $\\beta_1$ coefficients should be for a simple linear regression model using `Flanker_Unadj` as $Y$ and `FS_Total_GM_Vol` as $X$. Use the formulas we derived above ($\\beta_1 = Cov[XY]/Var[X]$ , $\\beta_0 = E[Y] - \\beta_1E[X]$). Then use `lm()` to compare the coefficients you calculated with the ones R gives you."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "id": "mWvD8shRDF5_",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 374
        },
        "outputId": "602464da-0485-443f-fe8f-bff269c3ca17"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[1] 90.25646\n",
            "[1] 3.109965e-05\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "\n",
              "Call:\n",
              "lm(formula = Flanker_Unadj ~ FS_Total_GM_Vol, data = d1)\n",
              "\n",
              "Residuals:\n",
              "     Min       1Q   Median       3Q      Max \n",
              "-28.2280  -6.5969  -0.0119   6.5497  31.0108 \n",
              "\n",
              "Coefficients:\n",
              "                 Estimate Std. Error t value Pr(>|t|)    \n",
              "(Intercept)     9.026e+01  3.041e+00  29.677  < 2e-16 ***\n",
              "FS_Total_GM_Vol 3.110e-05  4.412e-06   7.048 3.18e-12 ***\n",
              "---\n",
              "Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
              "\n",
              "Residual standard error: 9.864 on 1111 degrees of freedom\n",
              "Multiple R-squared:  0.0428,\tAdjusted R-squared:  0.04194 \n",
              "F-statistic: 49.68 on 1 and 1111 DF,  p-value: 3.177e-12\n"
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "beta_1 <- cov(d1$FS_Total_GM_Vol, d1$Flanker_Unadj) / var(d1$FS_Total_GM_Vol) # calculate \\beta_1\n",
        "beta_0 <- mean(d1$Flanker_Unadj) - beta_1 * mean(d1$FS_Total_GM_Vol) # calculate \\beta_0\n",
        "print(beta_0)\n",
        "print(beta_1)\n",
        "# now compare these to lm summary output\n",
        "summary(lm(Flanker_Unadj ~ FS_Total_GM_Vol, data = d1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xcnXbsZvDF6B"
      },
      "source": [
        "**DUE:** 5pm EST, Feb 28, 2024"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aG5swCweDF6B"
      },
      "source": [
        "**IMPORTANT** Did you collaborate with anyone on this assignment? If so, list their names here.\n",
        "> *Someone's Name*"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "R",
      "language": "R",
      "name": "ir"
    },
    "language_info": {
      "codemirror_mode": "r",
      "file_extension": ".r",
      "mimetype": "text/x-r-source",
      "name": "R",
      "pygments_lexer": "r",
      "version": "4.2.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}