.ipynb_checkpoints/section5.1-checkpoint.ipynb

   1 {
   2  "cells": [
   3   {
   4    "cell_type": "markdown",
   5    "metadata": {},
   6    "source": [
   7     "# Section 5.1: Using the model"
   8    ]
   9   },
  10   {
  11    "cell_type": "markdown",
  12    "metadata": {},
  13    "source": [
  14     "## Imports and defintions"
  15    ]
  16   },
  17   {
  18    "cell_type": "code",
  19    "execution_count": 1,
  20    "metadata": {},
  21    "outputs": [
  22     {
  23      "name": "stderr",
  24      "output_type": "stream",
  25      "text": [
  26       "── Attaching packages ─────────────────────────────────────── tidyverse 1.2.1 ──\n",
  27       "✔ ggplot2 2.2.1     ✔ purrr   0.2.4\n",
  28       "✔ tibble  1.4.2     ✔ dplyr   0.7.4\n",
  29       "✔ tidyr   0.8.0     ✔ stringr 1.2.0\n",
  30       "✔ readr   1.1.1     ✔ forcats 0.2.0\n",
  31       "── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──\n",
  32       "✖ dplyr::filter() masks stats::filter()\n",
  33       "✖ dplyr::lag()    masks stats::lag()\n"
  34      ]
  35     }
  36    ],
  37    "source": [
  38     "library(tidyverse)\n",
  39     "# library(cowplot)"
  40    ]
  41   },
  42   {
  43    "cell_type": "code",
  44    "execution_count": 2,
  45    "metadata": {},
  46    "outputs": [],
  47    "source": [
  48     "# Multiple plot function\n",
  49     "#\n",
  50     "# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)\n",
  51     "# - cols:   Number of columns in layout\n",
  52     "# - layout: A matrix specifying the layout. If present, 'cols' is ignored.\n",
  53     "#\n",
  54     "# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),\n",
  55     "# then plot 1 will go in the upper left, 2 will go in the upper right, and\n",
  56     "# 3 will go all the way across the bottom.\n",
  57     "#\n",
  58     "multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {\n",
  59     "  library(grid)\n",
  60     "\n",
  61     "  # Make a list from the ... arguments and plotlist\n",
  62     "  plots <- c(list(...), plotlist)\n",
  63     "\n",
  64     "  numPlots = length(plots)\n",
  65     "\n",
  66     "  # If layout is NULL, then use 'cols' to determine layout\n",
  67     "  if (is.null(layout)) {\n",
  68     "    # Make the panel\n",
  69     "    # ncol: Number of columns of plots\n",
  70     "    # nrow: Number of rows needed, calculated from # of cols\n",
  71     "    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),\n",
  72     "                    ncol = cols, nrow = ceiling(numPlots/cols))\n",
  73     "  }\n",
  74     "\n",
  75     " if (numPlots==1) {\n",
  76     "    print(plots[[1]])\n",
  77     "\n",
  78     "  } else {\n",
  79     "    # Set up the page\n",
  80     "    grid.newpage()\n",
  81     "    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))\n",
  82     "\n",
  83     "    # Make each plot, in the correct location\n",
  84     "    for (i in 1:numPlots) {\n",
  85     "      # Get the i,j matrix positions of the regions that contain this subplot\n",
  86     "      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))\n",
  87     "\n",
  88     "      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,\n",
  89     "                                      layout.pos.col = matchidx$col))\n",
  90     "    }\n",
  91     "  }\n",
  92     "}"
  93    ]
  94   },
  95   {
  96    "cell_type": "markdown",
  97    "metadata": {},
  98    "source": [
  99     "## Modelling abrasion loss\n",
 100     "\n",
 101     "This example concerns the dataset `rubber`, which you first met in Exercise 3.1. The experiment introduced in that exercise concerned an investigation of the resistance of rubber to abrasion (the abrasion loss) and how that resistance depended on various attributes of the rubber. In Exercise 3.1, the only explanatory variable we analysed was hardness; but in Exercise 3.19, a second explanatory variable, tensile strength, was taken into account too. There are 30 datapoints."
 102    ]
 103   },
 104   {
 105    "cell_type": "code",
 106    "execution_count": 3,
 107    "metadata": {
 108     "scrolled": true
 109    },
 110    "outputs": [
 111     {
 112      "data": {
 113       "text/html": [
 114        "<table>\n",
 115        "<thead><tr><th scope=col>loss</th><th scope=col>hardness</th><th scope=col>strength</th></tr></thead>\n",
 116        "<tbody>\n",
 117        "\t<tr><td>372</td><td>45 </td><td>162</td></tr>\n",
 118        "\t<tr><td>206</td><td>55 </td><td>233</td></tr>\n",
 119        "\t<tr><td>175</td><td>61 </td><td>232</td></tr>\n",
 120        "\t<tr><td>154</td><td>66 </td><td>231</td></tr>\n",
 121        "\t<tr><td>136</td><td>71 </td><td>231</td></tr>\n",
 122        "\t<tr><td>112</td><td>71 </td><td>237</td></tr>\n",
 123        "\t<tr><td> 55</td><td>81 </td><td>224</td></tr>\n",
 124        "\t<tr><td> 45</td><td>86 </td><td>219</td></tr>\n",
 125        "\t<tr><td>221</td><td>53 </td><td>203</td></tr>\n",
 126        "\t<tr><td>166</td><td>60 </td><td>189</td></tr>\n",
 127        "\t<tr><td>164</td><td>64 </td><td>210</td></tr>\n",
 128        "\t<tr><td>113</td><td>68 </td><td>210</td></tr>\n",
 129        "\t<tr><td> 82</td><td>79 </td><td>196</td></tr>\n",
 130        "\t<tr><td> 32</td><td>81 </td><td>180</td></tr>\n",
 131        "\t<tr><td>228</td><td>56 </td><td>200</td></tr>\n",
 132        "\t<tr><td>196</td><td>68 </td><td>173</td></tr>\n",
 133        "\t<tr><td>128</td><td>75 </td><td>188</td></tr>\n",
 134        "\t<tr><td> 97</td><td>83 </td><td>161</td></tr>\n",
 135        "\t<tr><td> 64</td><td>88 </td><td>119</td></tr>\n",
 136        "\t<tr><td>249</td><td>59 </td><td>161</td></tr>\n",
 137        "\t<tr><td>219</td><td>71 </td><td>151</td></tr>\n",
 138        "\t<tr><td>186</td><td>80 </td><td>165</td></tr>\n",
 139        "\t<tr><td>155</td><td>82 </td><td>151</td></tr>\n",
 140        "\t<tr><td>114</td><td>89 </td><td>128</td></tr>\n",
 141        "\t<tr><td>341</td><td>51 </td><td>161</td></tr>\n",
 142        "\t<tr><td>340</td><td>59 </td><td>146</td></tr>\n",
 143        "\t<tr><td>283</td><td>65 </td><td>148</td></tr>\n",
 144        "\t<tr><td>267</td><td>74 </td><td>144</td></tr>\n",
 145        "\t<tr><td>215</td><td>81 </td><td>134</td></tr>\n",
 146        "\t<tr><td>148</td><td>86 </td><td>127</td></tr>\n",
 147        "</tbody>\n",
 148        "</table>\n"
 149       ],
 150       "text/latex": [
 151        "\\begin{tabular}{r|lll}\n",
 152        " loss & hardness & strength\\\\\n",
 153        "\\hline\n",
 154        "\t 372 & 45  & 162\\\\\n",
 155        "\t 206 & 55  & 233\\\\\n",
 156        "\t 175 & 61  & 232\\\\\n",
 157        "\t 154 & 66  & 231\\\\\n",
 158        "\t 136 & 71  & 231\\\\\n",
 159        "\t 112 & 71  & 237\\\\\n",
 160        "\t  55 & 81  & 224\\\\\n",
 161        "\t  45 & 86  & 219\\\\\n",
 162        "\t 221 & 53  & 203\\\\\n",
 163        "\t 166 & 60  & 189\\\\\n",
 164        "\t 164 & 64  & 210\\\\\n",
 165        "\t 113 & 68  & 210\\\\\n",
 166        "\t  82 & 79  & 196\\\\\n",
 167        "\t  32 & 81  & 180\\\\\n",
 168        "\t 228 & 56  & 200\\\\\n",
 169        "\t 196 & 68  & 173\\\\\n",
 170        "\t 128 & 75  & 188\\\\\n",
 171        "\t  97 & 83  & 161\\\\\n",
 172        "\t  64 & 88  & 119\\\\\n",
 173        "\t 249 & 59  & 161\\\\\n",
 174        "\t 219 & 71  & 151\\\\\n",
 175        "\t 186 & 80  & 165\\\\\n",
 176        "\t 155 & 82  & 151\\\\\n",
 177        "\t 114 & 89  & 128\\\\\n",
 178        "\t 341 & 51  & 161\\\\\n",
 179        "\t 340 & 59  & 146\\\\\n",
 180        "\t 283 & 65  & 148\\\\\n",
 181        "\t 267 & 74  & 144\\\\\n",
 182        "\t 215 & 81  & 134\\\\\n",
 183        "\t 148 & 86  & 127\\\\\n",
 184        "\\end{tabular}\n"
 185       ],
 186       "text/markdown": [
 187        "\n",
 188        "loss | hardness | strength | \n",
 189        "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
 190        "| 372 | 45  | 162 | \n",
 191        "| 206 | 55  | 233 | \n",
 192        "| 175 | 61  | 232 | \n",
 193        "| 154 | 66  | 231 | \n",
 194        "| 136 | 71  | 231 | \n",
 195        "| 112 | 71  | 237 | \n",
 196        "|  55 | 81  | 224 | \n",
 197        "|  45 | 86  | 219 | \n",
 198        "| 221 | 53  | 203 | \n",
 199        "| 166 | 60  | 189 | \n",
 200        "| 164 | 64  | 210 | \n",
 201        "| 113 | 68  | 210 | \n",
 202        "|  82 | 79  | 196 | \n",
 203        "|  32 | 81  | 180 | \n",
 204        "| 228 | 56  | 200 | \n",
 205        "| 196 | 68  | 173 | \n",
 206        "| 128 | 75  | 188 | \n",
 207        "|  97 | 83  | 161 | \n",
 208        "|  64 | 88  | 119 | \n",
 209        "| 249 | 59  | 161 | \n",
 210        "| 219 | 71  | 151 | \n",
 211        "| 186 | 80  | 165 | \n",
 212        "| 155 | 82  | 151 | \n",
 213        "| 114 | 89  | 128 | \n",
 214        "| 341 | 51  | 161 | \n",
 215        "| 340 | 59  | 146 | \n",
 216        "| 283 | 65  | 148 | \n",
 217        "| 267 | 74  | 144 | \n",
 218        "| 215 | 81  | 134 | \n",
 219        "| 148 | 86  | 127 | \n",
 220        "\n",
 221        "\n"
 222       ],
 223       "text/plain": [
 224        "   loss hardness strength\n",
 225        "1  372  45       162     \n",
 226        "2  206  55       233     \n",
 227        "3  175  61       232     \n",
 228        "4  154  66       231     \n",
 229        "5  136  71       231     \n",
 230        "6  112  71       237     \n",
 231        "7   55  81       224     \n",
 232        "8   45  86       219     \n",
 233        "9  221  53       203     \n",
 234        "10 166  60       189     \n",
 235        "11 164  64       210     \n",
 236        "12 113  68       210     \n",
 237        "13  82  79       196     \n",
 238        "14  32  81       180     \n",
 239        "15 228  56       200     \n",
 240        "16 196  68       173     \n",
 241        "17 128  75       188     \n",
 242        "18  97  83       161     \n",
 243        "19  64  88       119     \n",
 244        "20 249  59       161     \n",
 245        "21 219  71       151     \n",
 246        "22 186  80       165     \n",
 247        "23 155  82       151     \n",
 248        "24 114  89       128     \n",
 249        "25 341  51       161     \n",
 250        "26 340  59       146     \n",
 251        "27 283  65       148     \n",
 252        "28 267  74       144     \n",
 253        "29 215  81       134     \n",
 254        "30 148  86       127     "
 255       ]
 256      },
 257      "metadata": {},
 258      "output_type": "display_data"
 259     }
 260    ],
 261    "source": [
 262     "# Load the data from the file.\n",
 263     "rubber <- read.csv('rubber.csv')\n",
 264     "rubber"
 265    ]
 266   },
 267   {
 268    "cell_type": "markdown",
 269    "metadata": {},
 270    "source": [
 271     "The figures below show scatterplots of abrasion loss against hardness and of abrasion loss against strength. (A scatterplot of abrasion loss against hardness also appeared in Solution 3.1.) Figure (a) suggests a strong decreasing linear relationship of abrasion loss with hardness. Figure (b) is much less indicative of any strong dependence of abrasion loss on tensile strength."
 272    ]
 273   },
 274   {
 275    "cell_type": "code",
 276    "execution_count": 4,
 277    "metadata": {},
 278    "outputs": [],
 279    "source": [
 280     "# ggplot(rubber, aes(x=hardness, y=loss)) + geom_point()"
 281    ]
 282   },
 283   {
 284    "cell_type": "code",
 285    "execution_count": 5,
 286    "metadata": {
 287     "scrolled": false
 288    },
 289    "outputs": [
 290     {
 291      "data": {
 292       "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAADAFBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUW\nFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJyco\nKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6\nOjo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tM\nTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1e\nXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29w\ncHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGC\ngoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OU\nlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWm\npqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4\nuLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnK\nysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc\n3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u\n7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////i\nsF19AAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3deWAU5f3H8UHkVFGRKt6I1Wqx\ntYrHT6uotSq2JgICYkQDEiqtiFRKQargLVgFq4IKolWrVpSKgggVKiCHHEFAlPu+QhZyn5tN\nnt/MPjNPssnm2ZnneYadmf28/8huJsM3z2b3RfbKrkYQQtJpyV4AQkEIkBBSECAhpCBAQkhB\ngISQggAJIQUBEkIKAiSEFCQNqSwvUQWVCXdxUkl1idJ54QKV0/KrFZ/aUqXjKqrzVY7LD6uc\nlles+KytVHrW5sU7a9VBKg0lKi+ccBcnFZNipfOqDqucdoioPbUlJUrHVZJDKscdrlI5LVSk\n+KwN5ykdF++sBSQWIIkHSIDEAiTxAAmQWIAkHiABEguQxAMkQGIBkniABEgsQBIPkACJBUji\nARIgsQBJPEACJBYgiQdIgMQCJPEA6QhC2jfz34tVnhxAkgiQJEoupNlna5p22y51JweQJAIk\niZIKaVN7zegedScHkCQCJImSCunFqCPt6B3KTg4gSQRIEiUV0nAKSVuh7OQAkkSAJFFSIf2D\nOmqu7kYSIEkESBIlFdL2s6KQ7ld3cgBJIkCSKLn32i3opDu6e5+6kwNIEgGSREl+HCln8Zff\nqzw5gCQRIEmEZzZwAyTxAAmQWIAkHiABEguQxAMkQGIBkniABEgsQBIPkACJBUjiARIgsQBJ\nPEACJBYgiQdIgMQCJPEACZBYgCQeIAESC5DEAyRAYgGSeIAESCxAEg+QAIkFSOIBEiCxAEk8\nQAIkFiCJB0jqIFVGElaTeBcHVZNqpfPUri5CFJ9axSeWKB2X8mdtWB2k0kOJyqtKuIuTSkix\n0nlVeSqnHSZhleMOlST++TqpkhxWOU7xWVtMSpTOC+crHRfvrFUIKeFvRFy1kwhX7STy1VU7\nQIoJkCQCJG6AJBEgSQRIvOpBWvn44Ff2yswDJPEAKTCQXm2haVqHNRLz5CC998Dg9+t+DkgS\nARI3NyGtahV9IdcuEvNkIB3sanz33x2s3QJIEgESNzchPWO+Rv9G8XkykJ6m3/2Z2i2AJBEg\ncXMT0gj5N7uQgXQF/e5X1m4BJIkAiZubkN6ml+Q2Eq8tLgOpE/32nWq3AJJEgMTNTUg5/xe9\nJI+TmCcDqTeF1Kd2CyBJBEjcXL3XblNGa+2MF2TmyUDKbmM4On517RZAkgiQuLn8gOzB7XLz\npO7+/vq6Fi2uX1hnAyBJBEjcgv3MhpycmE8BSSJA4hZsSPUCJIkAiRsgSQRIEgESL0CSCJAk\nAiRugCQeIAESC5DEAyRAYgGSeIAESCxAEg+QAIkFSOIBEiCxAEk8QAIkFiCJB0iAxAIk8QAJ\nkFiAJB4gARILkMQDJEBiAZJ4gARILEASD5AAiQVI4gESILEASTxAAiQWIIkHSIDEAiTxAAmQ\nWIAkHiABEguQxAMkQGIBkniABEgsQBIPkACJBUjiARIgsQBJPEACJBYgiQdIgMQCJPEACZBY\ngCQeIAESC5DEAyQHkPY8cVffcfo/iEwdkDkxXHsISPECJIkCDSk8cOyW5cOHETK53/LsrPG1\nh4AUL0CSKNCQNqUVE7I2rbys12JCVnUvsA4BKW6AJFGgIVWXk/Ltkx4mG9JKCKlKX20dGr+s\nluttK0hUcVXCXZxURsqUzosUqZxWSNSe2vJypeOqSKHKcUURldOUn7VVxUrHxTtr7UPSG5F2\n126ytLtxNGOedah/yOus94aNAQgFswg7ZgdS0cH37i5b0sM4mjHXOtQ/lL2st6QsUeXVCXdx\nUphUKp1XXa50HFF8asNKx0WI0lOr+KytJGpPrftnbYl9SDuz9Q81PZdvSCvTAaZnW4fW13Eb\nKSbcRpIo0LeRvu6r//oqSc8u7bmckHXd8qxDQIobIEkUaEiFGRO2/Dj6/gry+qCt24ZMIOwQ\nkOIFSBIFGhLZNPLOe58/qF+dm9w/c1K49hCQ4gVIEgUbUoIAKSZAkgiQuAGSRIAkESDxAiSJ\nAEkiQOIGSOIBEiCxAEk8QAIkFiCJB0iAxAIk8QAJkFhehrR1yJVXPrRV4UBAkgiQuHkY0vZz\nNL0O29VNBCSJAImbhyE9oEUbpG4iIEkESNw8DOkSCumX6iYCkkSAxM3DkC6lkC5WNxGQJAIk\nbh6GNIRCGqxuIiBJBEjcPAxp53mGo5/uVDcRkCQCJG4ehhTaObzLtcN3KBwISBIBEjcvQ8ID\nslIBEi9AkgiQJAIkboAkHiABEguQxAMkQGIBkniABEgsQBIPkACJBUjiARIgsQDJdvveefqf\ne+tuACRAYgGS3b7poGna2QvrbAEkQGIBks0OXBh95t95+2o3ARIgsQDJZrPpU2i1GbWbAAmQ\nWIBks/dMSG/WbgIkQGIBks2WmJAW1G4CJEBiAZLdukUd/b7OFkACJBYg2W1736O0o+6q+6JG\ngARILECy365vYv/IEJAAiQVI4gESILEASTxAAiQWIIkHSIDEAiTxAAmQWIAkHiABEguQxAMk\nQGIBkniABEgsQBIPkACJJQlp3fcxnwKSRIDELciQPuqgaR2n1dkASBIBErcAQ/pfS+NpnC3r\nPB8akCRKaUiV4URV1STcxUkRElE6r6ZK/N92p39YcEedTUTxqVV8YonSceHgnrVxinPWVqiD\nVJaXqIKqhLs4qZSUKp1XVSD+b39OIXWq3ZJPwvJLqlNp4p+vk8IkX+W4gojKaXklHjpr40Qa\nXpAPq4OU0lftrqGQutRuwVU7iVL6ql1KQ3qVQppUuwWQJAIkbgGGFBpoOLq/zgZAkgiQuAUZ\nUmjh2LGL6n4OSBIBErdAQ6ofIEkESNwASSJAkgiQeAGSRIAkESBxAyTxAAmQWIAkHiABEguQ\nxAMkQGIBkniABEgsQBIPkACJBUjiARIgsQBJPEACJBYgiQdIgMQCJPEACZBYMZBynjuv5fnP\n54hPAySJAImbnyA9GP37oj+LTwMkiQCJm48graZ/qNd0jfA0QJIIkLj5CNI/zbdJfU94GiBJ\nBEjcfATpIxPSJ8LTAEkiQOLmI0g72kYdtdslPA2QJAIkbsmFtGfe/L3cHWLubPjAeMXHlh+J\nrIsGSBIBErekQpqk/4456Q3eHrGPI2U/3HPYaqF10QBJIkDilkxIn0evqrX4grNLKj8gu2na\ntE0y4wApVSDdTO88+B1nlxSGNPYYTWv9rMQ4QEoVSBdSSD/n7JK6kKbRH86/xccBUqpAup5e\nVm7k7JK6kMxf1zeIjwOkVIFkPsL6L84uqQvpIvrD+Zn4OHchbXrl0XcPyMxLMqTIzM8KgwIp\nNKqFprV4jLdH6kK6JfGv6wS5CunDE/TFXfh9o3snLnmQSrLOJ+Q2Teu4KyiQQmvfemsdd4fU\nhfQfCulj8XFuQvrxRNkrnkmENEzrTZZqWZ+3HRgYSAlLXUih8W007bgXJMa5CWmC+XSt9eLz\nkgepw22EjGpRQO7rCEhi+QpSaPuMT7fJjHMT0mMmpEWN75+o5EFq+RQhXa4lZFxLQBLLX5Bk\ncxPSO9RR853i85IH6dw7yJ6mjxNy75mAJBYgSRQD6cBlUUh/lZiXPEgjjn7o0qN+LB3fug8g\niQVIEsXea/d92lHaMX+V+LP/JEIqur1Jk6fIRu2czYAkFiBJVP8B2d3ZMoyS+zhSYREhBfNK\nbDvyO6T53S66/pVcdfMASaKAPbMhUA/IJog+3WyAuoGAJFFwIAXwAVluB0+ndw3NUTYRkCQK\nDqRUe0B2pflYxWhlEwFJouBASrUHZFeZkB5XNhGQJAoOpHgPyOaPv7fPmB36TaepAzInhmsP\ngwAp92wKaZ6yiYAkUXAgxXtA9tEh6zaNzcgjk/stz84aT9hhECCFPmtuOHpA3UBAkig4kOI8\nIHsobYP+WyhjTlmvxYSs6l5gHQYDUmjJ3Vfe9rbCeYAkUXAgxXlANvcD/WpcRc/ZG9JKCKlK\nX20d6l+q3quXm5eogqqEuziplJQqnRfJVzktn4RVjssrVXtiw0TtqY2onJZXovisrSpQOo40\nvCAfbgRSIw/IVoztX7S0u3EsY551qH/I66z3BkHeKWdZTrKXkFJF2LH6kGp2zJuzvTp20/z+\nIwvIkh7G8Yy51qH+oWSE3pyKRFVWJ9zFSVWkSum8mkql44jiU+vkxO6/Q7/B12MfZ49qIrug\nmCprlI5TfdZWu37WljcG6b+/NG58d/pvnU0Fj9y3oIaQDWllOsD0bOvQ+rLPbyMF6rl29G/H\nb+Q84wm3kSRycBtpZbPTn/zPjGfOaMackJo/P1UaBdNzOSHruuVZh4AUt2RCWmg+Kja/8V0A\nSSIHkG45+5BxcLjDrWzTmvQFa/RC5PVBW7cNmUDYISDFK5mQ3jUhTW18F0CSyAGkU0bRw0fb\ns02fpkWbRSKT+2dOMh6QNQ8BKV7JhPSlCWlW47sAkkQOIJ1sQTqF2A2QYkompJzoDVztIs5r\nwQGSRE6u2nWIXrXL63grsRsgxZTUOxuWG6/QfMG3nD0ASSIHkFY0O/3pGTOePbPZCkASK7nP\nbDjwnwnTua9NCkgSOfnDvrmdoq85/6VtR4AUG54iJFGAIJHqbXPnbKkm9gOkmABJvO/HDBgt\n/n7yccKL6PMCJIm8DGnGsfqVodYSbzvaoORAuiYmQBILkETbeyp9I+wd6kYCEi9AksjDkD4z\nHwX7UN1IXLXjBUgSeRjSByakN9WNBCRegCSRhyGtbUohLVc3EpB4AZJEHoYUGhJ1NFDhREDi\nBUgSeRnSgac6NDlr9H6FEwGJFyBJ5GVIoVBRJDgPyAKSXIAkUZCe2QBIcjWAtEduHiBJBEi8\n/AQpd8LZ2vH3bZWYB0gSARIvP0EaF70j6rqD4vMASSJA4uUjSPuOow+N/Et8HiBJBEi8fATp\nW/PB+lHi8wBJIkDi5SNI601Iz4nPAySJAImXjyCFrow6apUtPg+QJAIkXn6CtMp4R8AWr0rM\nAySJAImXnyCFdr+Y9egKmXmAJBEg8fI2pMXvTtukch4gSQRIvLwMKedO/apcm9fUDQQkmQCJ\nl5chjYzeudByofCAVeOfnBmzAZAkAiReXoZEX3ZAu1/03z/XQv/XN+2rswWQJAIkXl6GdDSF\n1E3wn8+h//zBOpsASSJA4uVlSB2phD8L/vMs+s/b1dkESBIBEi8vQ3olCuHEtYL/vAeF1LTO\nO4UBkkSAxMvLkEJjWmvaeZy3VeH3VwrpZ3U2AZJEgMTL05BCu1f8kCP8jzeeEoX0dp1NgCQR\nIPHyNiS5PzVfcLl+C+nlulsASSJA4hVkSKHQ5uzYt1IGJIkAiVewIdUPkFj7n7jml/esdjIP\nkHilMqQVvX926cjdEuN8DOngDcbtx2OXOZgHSLxSGNLSY4zL0pXc9+Tj52NIE+ldml0czAMk\nXikM6Xp6WXpBfJyPIfWlJ76Zg1eSASReKQypJb0s9RQf539IzQFJUSkM6Rh6WeojPs7HkN6g\nJ/4GB/MAiVcKQ/o9vSy9Lj7Ox5BybzZO+/ErHcyrByn3B4kXGAwBUoJ8BGltO+OydGtuo3sn\nzMeQQgee/83lA793Mi8G0t4hrbVWg3aJrSyau5AqElZZnXgfB1WRKqXzqiuVjiOKT23Mid3/\n1xu6vV4mMa7axvnloMoapeNcPWvpU+vvkhgX56wtVweprCBRRVUJd3FSmY1v6aRIkcpphUTt\nqS0vVzquihSqHFcYUTlN+VlbVeesXWO+xuAS8XHxzlp1kHDVLiY8s0EiN58iZL1F7WTxcbiN\nxA2QxPMRpNkmpGni4wCJGyCJ5yNI++lfL58h8QQrQOIGSOL5CFJovvHHXifNlhgHSNwASTw/\nQQrtePnhCTLv+QZI/ABJPF9Bkg6QuAGSeIAESCxAEg+QAImVypA+HjpU4u5gQAKkOqUupNxu\nxh3CaRLP5AQkQGKlLiT6nuvas+LjAAmQWKkL6WoK6QrxcYAESKzUhfQLCulC8XGABEis1IV0\nJ4V0h/g4QAIkVupCWnVc9AWuJN7jFpAAiZW6kELzrmp29P99JTEOkACJlcKQQqF9+xrZ0U45\n4y4/5ybht96IEyBxAySJvPzMhrujN7HeVTcQkLgBkkQehmT+qdxPJF74tV6AxA2QJPIwpMfN\nvzl18mrc/ACJGyBJ5GFIT5uQViqbCEjcAEkiD0NaRB11lHiZvXoBEjdAksjDkOhb3raU+ePt\negESN0CSyMuQQtN63zgwW+E8QOIGSBJ5GhIekAWk2gBJPEACJBYgiQdIgMQCJPEACZBYgCQe\nIAESC5DECzSkFa/8I/ZJG4DELUiQ9rz1+Ft7eDsAku2GN9c07YG6WwCJW4AgLTpTP+/PXMDZ\nox6kOWNGfym2LlqAIb1Dn7fxap1NgMQtOJD2nx8973/K+RujWEj3GrvfLfEUnwBDuolCurLO\nJkDiFhxIX5jPIp3R+C4xkF6mu48XX12AIV1Cfzjn1tkESNyCA+ldE9Jbje8SA+lauvv/ia8u\nwJB60B/OTXU2ARK34EBaYkJa1Pgu8V6O6wLx1QUY0tctoz+cuk/JBSRuwYFk/i+aztkjBpL5\nny5v/wQFGFLo/TM1rX3ML3dA4hYgSDsyj9aa3rOds0cMpMWtDEetOL/AEhVkSKGDK77NidkA\nSNwCBCkU2ruU+zBSvXvtPvu5pl34qcCqrAINqUGAxC1QkBJV/wHZLZulxgESILFSGpJkgARI\nLEASD5AcQarKKNI/RqYOyJwYrj0EpHgBkkTBhlS5dlyaAWlyv+XZWeNrDwEpXoAkUbAhTe/f\n14BU1msxIau6F1iHgBQ3QJIo2JAI2WJA2pBWol/JS19tHerbawr18g8lKj+ccBcnlZBipfOq\n8lROO0zUntrSUqXjKslhlePyqlROO1RMShztv2sn/+vhxBdOJ8U7a51DWtrdOJoxzzrUP+R1\n1nvDxgCE1Lf08ibaJQuSu4YIO2Yb0pIextGMudah/qH4j3qfhxNWk3gXB0VIROk8tasLE8Wn\nVvGJJUrHJfOs3RB9m7TW33F2cf+srXAOaUNamQ4wPds6tL6I20gx4TaSRI5uI/WhzxS8jbOL\nJ28jlfZcTsi6bnnWISDFDZAkcgTJ/IOh8zi7eBISeX3Q1m1DJtQeAlK8AEkiR5Cup5Au5+zi\nTUiRyf0zJ4VrDwEpXoAkkSNIr1BIz3F28RokboAUEyBJ5OxxpHsMRz15rzgBSLwASaIgQQrN\n/tvfZnJ3ACRegCRRoCAlDJB4AZJEgCQRIHEDJPEACZBYgCQeIAESC5DEAyRAYgGSeIAESCxA\nEk8O0sbh6QO+qLsBkLgBkkQBhvTNCcYjqo/W2QJI3ABJogBDupg+x2dh7RZA4gZIEgUX0g/m\nK5U/VrsJkLgBkkTBhbTahDSidhMgcQMkiYILKedkCumj2k2AxA2QJAoupNDbUUe/q7MFkLgB\nkkQBhhSa9n8nnj+y7ov+AxI3QJIoyJAaBEjcAEkiQJIIkHgBkkSAJBEgcQMk8QAJkFiAJB4g\nARILkMQDJEBiAZJ4gARILEASD5AAiQVI4gESILEASbxAQdr8r0mLuTsAEjdAEi9IkKYaf2fY\n+wBnD0DiBkjiBQjS0lbR59AO4+wCSNwASbwAQXqI/lFHW84ugMQNkMQLECTzfcu0/Y3vAkjc\nAEm8AEEaQR2dztkFkLgBkngBgvR92yik5zm7ABI3QBIvQJBCs8/TtJYjeXsAEjdAEi9IkEIH\nln65nbsDIHEDJPHchfRxnxv+sFpmHp7ZwAuQJPITpEeM2ywtZ0vMAyRegCSRjyB9Q+9G68h7\nk9gEARIvQJLIR5CeMR/YWSE+D5B4AZJEPoL0hAlpifg8QOIFSBL5CNJs6qgd72mkCQIkXoAk\nkY8ghe6OQvqnxDxA4gVIEvkJUs64y8688XOZeb6CVFGWqPLqhLs4KUwqlc6rLlc6jig+tWGl\n4yJE6alVfNZWErWn1v2ztkQdpPKiRJVEEu7ipHIb39JJkWKV04pJlcpxRRUVSsdVEbWn1uNn\nbUnMp3k/HpYaRxqe2kJ1kHDVLiZctZPI1acI7R7UXGt+3w6JcbiNxA2QxPMTpIzofRe3S4wD\nJG6AJJ6PIK00H5b6WnwcIHEDJPF8BOkDE9Jk8XGAxA2QxPMRJPPxXe1j8XGAxA2QxPMRpAPn\nRR2dtafx3RMFSNwASTwfQQotOF13dMpciXGAxA2QxPMTpNDu10ZMlLn3G5D4AZJ4voIkHSBx\nAyTxAAmQWIAkHiABEguQxAMkQGIBknjuQlr4UJ/RW2TmARIvQJLIT5BeaG68qP03EvMAiRcg\nSeQjSCtbRh9B/YXEPEDiBUgS+QjSWPM5PRIvEQlIvABJIh9BGmNCkrhuB0i8AEkiH0H6hDpq\ns098HiDxAiSJfAQp9LsopJck5gESL0CSyE+Qdg09s1kniT8fAiR+gCSRnyDJB0i8AEkiQJII\nkLgBkniABEgsQBIPkACJBUjiARIgsQBJPEACJBYgiQdIgMQCJPEACZBYgCQeIAESC5DEAyRA\nYgGSeIAESCxAEg+QjiSkQ2p/OqumZCudV6z0onVgitSbNTaooEDpuJlTJN7buGGKz9rvpqxU\nOk/tWXtwyqcNN6qDdMSb3vmzZC+BU0nnB5K9BF5DOhcm3ilpzeo8LdlL4BTuPJDzVUBSGyBJ\nBEhHMkCSCJDEA6QjGSBJBEhHssrCymQvgVNNYWmyl8CrtLAm2UvgFPb0WUu4Z63/ICHkwQAJ\nIQUBEkIKAiSEFOQnSB+n6XUjJDJ1QObEcLJX07B5f+796F6PLm9JWrSXvLk6QvLH35MxLuTR\nHx7JHXd3/3+UclfnJ0gvPZGdnb2akMn9lmdnjU/2aho0r9dXax+9v9qby8vXf3TZy/os9ebq\nCBk5fNmKUUM8et6WD3xy47rhj3JX5ydIwz+PHpT1WkzIqu4FSV5N/WoGzSIkNPagR5dnNGmy\nV394lenfEbIhLd+by1t6R4V+1qbt5K3OT5Aynux31xN79Z93CSFV6auTvZx67U47XGP8iD26\nPL3vBoY9u7qR4/YeGP+gR5f31Z01+q+l9IW81fkIUmHaU+vXjupXurS78VnGvGSvp17fdZve\nOy1zCfHo8gipHqz/f+rV1RVkpKXdGfLo8g72fKf08Itpn/FW5yNIkUP6/wsldyxY0sP4LGNu\nstdTr4Vpzxws/bj7bo8uT78Np98EIR5dXfngF3fufnVQsUeXt7J/Wo9/3fU1b3U+gkT70ycb\n0sp0VenZyV5Jvdak5ekfB3zm0eUR8tBsYlzx9OTqFveO6LcyM+d7dHmE5FVVpK/jrc5HkFYM\nLtL/6+r1bWnP5YSs65aX7PXUK5S+W/8p953n0eWRDT2Mp4p5dHULelXp1z3vmePN5RU8v0df\nYt8q3up8BKk0c8x3P4wZHCGvD9q6bciEZC+nQeOGrtnyQmaRV5c3dWT0wJurK8p8dtOmF+/O\n8+jyHvrruiV3T+f+8HwEiex87M57x+fr/+1P7p85yXMP2pHKif0zntrn2eX96V/RA4+ubu+z\nfTOe2OnV5R0c03uw8cc7nNX5CRJCng2QEFIQICGkIEBCSEGAhJCCAAkhBQESQgoCJIQUBEgI\nKQiQPFzXy/hff0Hz2F/ApXCA5OEAyT8BkocDJP8ESB6uEUhlK80jgOSdAMnDdb1s+23t2g8w\ntLx/xQnHXTLF2NZz1nEdCPng6jadJxqQunbbc/Mx7Qcar42/vffZbbp8oR8peuSnrTr+paTO\nEeR2gOThup52xuApPbQsQqZrVz47/Bfax/q2S0/sPVH/XXThqEGtzzEgXd3lkx2TmtxHyJo2\np414/KImbxLS7eg7nvy98c/YEeR2gOThumqTCam5uCMh3c+oJKSizR+MbW8REjruslJCljYx\nIGlfGXueRch1Zx0mJHz9ccWFTR7SN/U+n7AjyPUAycN1PTaif7y3PSGHjD9vDh3TV992QjUh\nn2ifGl//nQGprXFsQDuSpz1tHPtEm1fU5NK90X/PjiDXAyQP1/Ui42M/HRLZ8u7D17XQDEid\n9M+e03YYX3nEgPQr41hWO7JMM/uQPHlU0+tGLdM3syPI7QDJw9F77QxILzdr2/f11Wf2Nbf9\nnUJ61IAU3UeHlK2NXBDtACHrx1zTQkuL1DmCXA6QPByDVNIi08BwsgVpujbD+Eq3upAKtVHG\nsf0Lygs26jeg8rO0mexIsk5ACgVIHo5B+l57RT8yR8swtx1uc0UZId81rQuJ3Ngul5Dqm9pH\n5mnG67x/rn3GjiTvJKRMgOThGKTKM04d/c8/nXLGyW+b217UOo0Z2uaaGEirjz111GOXau+R\nknNaZz4/4KRzCtmRZJ6IFAmQPFztbaR1v21z1l07l3XJsp7t8MFVx13y8re/LTE/v/88/cOm\n7mcc/+tZxpHep7XokLWrzhHkdoCEkIIACSEFARJCCgIkhBQESAgpCJAQUhAgIaQgQEJIQYCE\nkIIACSEFARJCCgIkhBQESAgpCJAQUhAgIaQgQEJIQYCEkIIACSEFSUMqy7NXRXWBzT0dVlHo\nztyy6iKXBhe7M7ekusSlwaXuzC2qdmtwuTtzC6obDFYHqTRkrwpy2OaeDqsocGduGXFrcJE7\nc4tJsUuD7Z7HDiu0feFxWEG5O3PzSEX9TYCUMEBigwGJBkgiARIbDEg0QBIJkNhgQKIBkkiA\nxAYDEg2QRAIkNhiQaIAkEiCxwYBEAySRAIkNBiQaIIkESGwwINEASSRAYoMBiQZIIgESGwxI\nNEASSQ2k+Y/88bUD9QYDEg2QAMluo4x3R+60NXYwINEACZBsNpu+z3if2MGARAMkQLLZIAqp\nVW7MYECiARIg2exuCqnJvpjBgEQDJECy2TMU0vmxgwGJBkiAZLPd50chfRg7GJBogARIdltz\ne6smF7xTbzAg0QAJkOx3cE+DwYBEAyRAkhoMSDRAAiSpwYBEAyRAkhoMSDRAAiSpwYBEAyRA\nkhoMSDRAAiSpwYBEAyRAksRblGMAAB0ESURBVBoMSDRAAiSpwYBEAyRAkhoMSDRAAiSpwYBE\nAyRAkhoMSDRAAiSpwYBEAyRAkhoMSDRAAiSpwYBEAyRAkhoMSDRAqlt5kb2qSLHNPR1WVerO\n3Eri1mC7PzKHlds+L5wOrnBnbhmpdGdwadiduSWk/uBCdZAqyuwVIeU293RYxO4KHFZl+6Q5\nHVzpztwwcWtw2J25laTKncEVEXfmlpP6g0vUQcJVO8eDcdWOhqt2gCQ1GJBogARIUoMBiQZI\ngGSn7xfujj8YkGiAFChIm2ctzYm3XRLSsis1rfnQeJMByQyQAgQpZ1AzTfv5/DhfkYO089zo\nS0M+Em8wINEAKUCQRkQv76dtbvgVOUj/oC9WfMy+hl8CJDNACg6kA8fRC/yzDb8kB2konaut\njjMYkGiAFBxIG8zL+/0NvyQHyXz5/KN3xhkMSDRACg6kfS3pBf7xhl+Sg7Tu+OjcO+N8CZDM\nACk4kEIDo5f3E9Y1/IrkvXYfttXnXrs9zlcAyQyQAgRpz2365f2Uj+N8RfZxpK1vPz8r7hcA\nyQyQAgQpFFo0adqOeNvxzAY2GJBogCQSILHBgEQDJJEAiQ0GJBogiQRIbDAg0QBJJEBigwGJ\nBkgiARIbDEg0QBIJkNhgQKIBkkiAxAYDEg2QRAIkNhiQaIAkEiCxwYBEAySRAIkNBiQaIIkE\nSGwwINEASSRAYoMBiQZIIgESGwxINEASCZDYYECiAZJIgMQGAxINkEQCJDYYkGiAJBIgscGA\nRAMkkQCJDQYkGiCJBEhsMCDRAEkkQGKDAYkGSCIBEhsMSDRAEgmQ2GBAogGSSIDEBgMSDZBE\nAiQ2GJBogCQSILHBgEQDJJEAiQ0GJBogiQRIbDAg0QBJJEBigwGJBkgiARIbDEg0OUh7nrir\n7zj9H0SmDsicGK49BCTRwYBESy1I4YFjtywfPoyQyf2WZ2eNrz0EJNHBgERLLUib0ooJWZtW\nXtZrMSGruhdYh4AkPBiQaKkFqbqclG+f9DDZkFZCSFX6autQ/1LldL11xfaqIiU293RYVZk7\ncyuJS4PD5e7MrSAVLg2udGduGXFrcNiduaWkqt6WIvuQ9Eak3bWbLO1uHM2YZx3qH/I6671h\nYwBCwSzCjtmBVHTwvbvLlvQwjmbMtQ71D5Vf6W0oslcVKba5p8PCpe7MrSRuDS5zZ65+1cGl\nwRXuzNV/I7kzuDTsztwSUn9woX1IO7P1DzU9l29IK9MBpmdbh9bXcRvJ8WDcRqKl1m2kr/vq\nv75K0rNLey4nZF23POsQkIQHAxIttSAVZkzY8uPo+yvI64O2bhsygbBDQBIdDEi01IJENo28\n897nD+pX5yb3z5wUrj0EJNHBgERLMUgJAiTHgwGJBkiAJDUYkGiABEhSgwGJBkiAJDUYkGiA\nBEhSgwGJBkiAJDUYkGiAFFhI679Yw44DEhsMSDRAsteW2zVNu2m9+RkgscGARAMke92mGV19\nkH4GSGwwINEAyVbfarRZ9FNAYoMBiQZItppmQnqVfgpIbDAg0QDJVgtNSJ/QTwGJDQYkGiDZ\nKveKqKML9tFPAYkNBiQaINnru1/pjn62xPwMkNhgQKIBks0Oznzl0wPWJ4DEBgMSDZBEAiQ2\nGJBogCQSILHBgEQDJJEAiQ0GJBogiQRIbDAg0QBJJEBigwGJBkgiARIbDEg0QGrQzn0JdwEk\nNhiQaIBUr+mdmjS9ekGCnQCJDQYkGiDF9lUL46lAbdfx9wIkNhiQaIAU2w30yalZ/L0AiQ0G\nJBogxXYahXQNfy9AYoMBiQZIsV1AIf2evxcgscGARAOk2EZRSFP5ewESGwxINECK7cDNhqOB\nCfYCJDYYkGiAVL9pw0bNSbQPILHBgEQDJJEAiQ0GJBogiQRIbDAg0QBJJEBigwGJBkgiARIb\nDEg0QBIJkNhgQKIBkkiAxAYDEg2QRAIkNhiQaIAkEiCxwYBEAySR5CHtb2QwINEACZASt7n/\niU3Pnxx3MCDRAAmQEpZzVfSJsZPiDQYkGiDVraLcXhHbezosUunO3CoiM/h9+gzzdiVxBocl\n5nIKE7cGV7kzt5K4NTjiztwKUn9wqTpIZQX2CpMim3s6LFziztwKIjP4L+Z7xPwQZ3CpxFxO\nZbbPC6eDy92ZW0IqXBpc6c7cIhKuv0kdJFy1i9tj1FGTLXEG46odDVftAClhi6MvsKJdF28w\nINEACZAS91xz3dGZ38UbDEg0QAIkGy36631/3x13MCDRAAmQpAYDEg2QAElqMCDRAAmQpAYD\nEg2QAElqMCDRAAmQpAYDEg2QAElqMCDRAAmQpAYDEg2QAElqMCDRAAmQpAYDEg2QAElqMCDR\nAAmQpAYDEg2QAElqcKAg5axJ/L7WjQRIgCQ1OECQ9g1pqTXru1VsLiABktTgAEHKiv7V1c25\nQnMBCZCkBgcH0vdH0T8E/kJoLiABktTg4ECaYb4yxUtCcwEJkKQGBwfSQhPSO0JzAQmQpAYH\nB1LuL6KOTt0hNBeQAElqcHAghZaeZbx438y6m9a++Ld/HbQ1F5AASWpwgCCF9k4e8UrMvd+v\nt9Zp/WqznbmAlDqQ/jcsa/xe1YODBKl+K1tHr+zdbmdfQEoZSGOMC0WHHxQPDjKk0fTeh6a7\nbOwLSKkCaR69VHRVPDjIkB4078f73sa+gJQqkIaZ/73uUTs4yJBepj+yE3Ns7Bt0SJGZnxX6\nANKPf7isy+h9rkK63/zv1dZNZ/uDgwxpz4XRn9hYO/sGGFJJ1vmE3KZpHXd5HtK6k4wz7IoD\nbkJ6hTo6Q+ypZI0ODjKk0OpbjtJOfNrWTyzAkIZpvclSLevztgM9D+l2eiF/zk1I+y+Jfo9/\nKh4caEih0O41Nv/jCTCkDrcRMqpFAbmvo+chnUQh/d7Ve+029j2h6YWKHQUeku0CDKnlU4R0\nuZaQcS0ByXxAVvjP1hofDEi0AEM69w6yp+njhNx7puchpR+Bq3YuDQYkWoAhjTj6oUuP+rF0\nfOs+noe0tq3h6HJX72xwaTAg0QIMqej2Jk2eIhu1czZ7HlLoh4GXXPPYXjzXjgVIVsmHREhh\nESEF80psOwr0A7IuDQYkmkcg5S75eLW9PYP5gCwbDEg0QLJyBGnl5fpNhPRtdnYN5AOytYMB\niQZIVk4g7bsoeqdVdzv7BvIB2drBgEQDJCsnkD40nwa2xsa+gXxAtnYwINEAycoJpOdNSLNs\n7BvIB2RrBwMSDZCsnEB6z4Rk5/4GuQdk88ff22fMDkIiUwdkTgzXHgKS6GBAonkC0p7zoo5u\nsbOv3AOyjw5Zt2lsRh6Z3G95dtZ4wg4BSXQwINE8ASn0jfGHH9dtsrOr1AOyh9I26L+FMuaU\n9VpMyKruBdYhIAkPBiSaNyCFDsx+c4G9PaUekM39QL8aV9Fz9oY0fWNV+mrrUP9SyQi9ORX2\nqiaVNvd0WHXYnbkR4tbgKnfmVhG3BkfcmRsmbg12aW5lgwWXNwqpZse8Odur622sGNu/aGl3\n41jGPOtQ/5DXWe8Nglys8sXfXvXQ/mSvAjVShB2rB+m/vzRueXX6b91tNfP7jywgS3oYxzPm\nWof6h+q9erl59qokBTb3dFhlkTtzy4lbg0vs73uoi3F+tF1nZ99S4mCwk0rL3JlbTMrdGVxU\n4c7cAlJZb8vhRiCtbHb6k/+Z8cwZzbJrtxU8ct+CGkI2pJXpANOzrUPry7iN5Hiwg9tIE+h9\ns7fa2Re3kayS/6TVW84+ZBwc7nBr7e+jPz9VGgXTczkh67rlWYeAJDzYAaQ7KKRj7ewLSFbJ\nh3TKKHr4aHu2aU36gjV6IfL6oK3bhkwg7BCQRAcDEi3AkE62IJ3CNn2aFm0WiUzunznJeEDW\nPAQk0cEOII138BKVgGSVfEi3dIhetcvreCuxGyA5HuwAUs6vDUcn2Pp7GUCySj6kFc1Of3rG\njGfPbLYCkDwBKbTv8Wt/lWXn9X8BqbbkQyJzOxn/A/78S9uOAMn5YDyzgRZkSKR629w5W+o/\nIAtIClv/6WdbE+8lECBZeQGS4wDJWcOb67d5JroxGZCskgvpmpgAyR1IE6P3wrWY68JoQLIC\nJJH8BelX9P7sPi6MBiQrXLUTyV+QTqaQrnVhNCBZAZJI/oJ0KYWU4cJoQLICJJH8BWlK1FHL\n/7kwGpCsAEkkf0EKjW6lae2mujEZkKwASSSfQQpt/uK/dt7y23mAZAVIIvkNEp7ZYAVIgCQ1\nGJBogARIUoMBiZZsSLlvXN7+8tccvNM2IIkESGxwQCH9LXoP6iP25wKSSIDEBgcT0vpmUUhH\nr7U9F5BEAiQ2OJiQ3jdf9dv+u9YDkkiAxAYHE9I0E9IHtucCkkiAxAYHE9L246OO2tj/6zBA\nEgmQ2OBgQgq91Vx31HyK/bmAJBIgscEBhRRalPXbAYsczAUkkVID0txXp+9LODiokJwGSCKl\nAqRNxqt9dfgq0WBAogGSSKkA6dboze0zdyQY7M7lff97z725xZXJgCQSIJk5h/S9eQfw5ASD\nXYG02njPybb/cWM0IIkESGbOIc0zIT2ZYLArkK6Mfut2m10YDUgiAZKZc0ibjqKQ3kkw2A1I\nq0zEr7kwG5BEAiQzgdtI90Qvy50S3G/nCqT/mpCedWE2IIkESGYCkHbfpV+Ur1qZaLAbkDYf\nTSFNc2E2IIkESGZCjyNtmLkq4R/kuHMbaXDU0a8PujAakEQCJDOfPbNh/9BWWtM7NrgxGpBE\nAiQzn0EKhQ5vO+TOYEASyU1I2wd3OP7qT1UPBiRaUp/ZIHKVEpBE0iHR98tTfXsYkMySByl3\n/LlHnTZsj9O5gCSSDukNesfSOYoHAxIteZCeiJ6v3Z3OBSSRdEh/NB/rUPuMMEAySxqk7S3o\n+Trb4VxAEkmHNJT+vI/arXYwINGSBmmu+R/kWIdzAUkkHdJs+vO+XvFgQKIlDdJiE9KrDuf6\nC9Kid+c7ukfFzXvtor+STlmteDAg0ZIGKff8qKNj1zuc6ydIP3bRT+GlKx2cOlcfR5qR1WP0\nNtWDAYmWvDsb/tfWeP/RBH890jBXIVVG7FVDbO12c/T/ikvKbU41Blfb39dJNcStwS7NrXZr\nwdWuLbjGpcEJ5+aOGzB6o/PBDRYcVgep9JC9Kkmejb2Wmtde/2Nzql5Fgf19nVRGCl0aXOzO\n3BLi1mC757HDCkmZS4Mr3JmbTxoMVgjJ5q9Fe1ftPjYhvWL/9y2eImSGq3ZWeIpQaJkJycHT\ncgDJDJCsACkU+i29jXTA/qkDJDNAsgKkUGjjb3RHV2Q7OHWAZAZIVoBktOz9RQ7e+wmQWIBk\nBUgiAZIZIFkBkkiAZAZIVoAkEiCZAZIVIIkESGaAZJUY0puXtPnZ4/udzgUkkQCJDQ4cpOej\nD7Lc7XQuIIkESGxw0CDtak0f9k/0Nhz1AySRAIkNDhok61VexzmcC0giARIbHDRIi0xILzuc\nC0giARIbHDRIB8+OOmq1xuFcQBIJkNjgoEEKzT7GgPSS07mAJFLQIS2YNG2nzcGBgxRaN+z2\nPy5wPBeQRAo2pN3Gu16eau/98gIISSxAEinYkPpFbyW0tfXyH4BkBkgiBRrS3ub0fqunbQ0G\nJBogiRRoSOvMO4AH2xoMSDRAEinQkPaZj+3/3dZgQKIBkkiBhhQaFnV05lZbgwGJBkgiBRvS\ngT8crWm/sHcPMCCZAZJIwYYUCm2euczmS0MDkhkgiRR0SA4GAxINkEQCJDYYkGiAJBIgscGA\nRAMkkQCJDQYkmh8h5ay3+yZJgGQGSFaAZLV9QAvtmKH77A0GJJrvIIU2OXpHOfsBklVa9OHC\nLHuDAYnmM0j7HmyuNUn/0Y3RaiAdWPpNvdcZ8h2kr+jzV5rY+gtGQDLzGST6lvFX5bgwWgmk\nd0/TtJ/Evqmf7yBNNJ9SOc3WYECi+QvSpqb0PP7IhdkqIM1vEV3ezLrbfAfpfc3ByyUBkpm/\nIFmv5POMC7NVQOpOl3dj3W2+g7Tz1OiJ+JmtX/uAZOYvSNkmpDdcmK0C0iV0eefW3eY7SKGZ\nxrtOn7bI3mBAovkLUujq6AX15C0ujFYBib4ruHZV3W3+gxTa8tKwV3fZHAxINJ9B+u4C/XL6\nkxlujFYB6T0K6bW623wIyUGAZOYzSKED//7729tcmazkXruRzTWt+ZCYTYAkEiCxwan5zIbs\nSa+siN0CSCIBEhucmpAaBkgiARIbDEg0WUhVGUX6x8jUAZkTw7WHgCQ6GJBoKQapcu24NAPS\n5H7Ls7PG1x4CkuhgQKKlGKTp/fsakMp6LSZkVfcC6xCQhAcDEi3FIBGyxYC0Ia1Ev5KXvto6\nBCThwYBES01IS7sbRzPmWYf6h/wb9P5ZYy9CbO7oNNfmyg3OH3pas0umK1qLvdz6SbiW/xZc\nf8VVziEt6REFNNc61D8Uput9ELFXDbG5o9Nqql2aS2QGh2+IPjAe76fj1oKrpRbMG1zj0lzi\n1mCX5kYaLLj2TjcHV+3KCImkZ1uH1hdx1S5ub9NnmPwkzvNucdXOLDWv2pX2XE7Ium551iEg\ncXvYfHJznL9NBCSz1IREXh+0dduQCbWHgMTrUfOPfDfHGQxItBSFFJncP3NSuPYQkHgton9g\n+et4gwGJlnKQuAFS/J4yHLXPjjcYkGiABEg2mj8045kdcQcDEg2QAElqMCDRAAmQpAYDEg2Q\nAElqMCDRAAmQpAYDEg2QAElqMCDRAAmQpAYDEg2QAElqsB8hbR927Q2j96qdC0iAJDXYh5C2\nnm08vnyxvTfasRsgAZLUYB9CGkCfOjhK6VxAAiSpwT6EdB6FdI3SuYAESFKDfQjpXArpaqVz\nAQmQpAb7EFJfCukvSucCEiBJDfYhpI3to++zs1vpXEACJKnBPoQU2viHizsP3a52LiABktRg\nP0JyI0ACJKnBgEQDJECSGgxINEACJKnBgEQDJECSGgxINEACJKnBgEQDJECSGgxINEACJKnB\ngEQDJECSGgxINEByBOmjPjc+sBaQ6gwGJBogOYEUfTX5Y/4HSLWDAYkGSA4gzadPG74QkGoH\nAxINkBxAesx8e5N1gMQGAxINkBxAesSEtBqQ2GBAogGSA0jTqaPTDgISGwxINEByACnUPQrp\nQ9xGqh0MSDRAcgJp/1OXnP7bL3D3d53BgEQDJCeQrACJDQYkmguQ1i7cEwIksQDJDJAWX65p\nLf5yEJCEAiSzlIe0PfpysNpoQBIKkMxSHtLz9B7iNgcASSRAMkt5SH8yH7NcD0giAZJZykMa\nTR013wtIIgGSWcpDWn1cFFJf3EYSCpDMUh5S6L0TdUfX7/QppEP2qiR5Nvd0WEWBO3PLSKFL\ng4vdmVtC3Bps9zx2WCEpUzxxy1svzDYGVyiea5ZPGgxWB6myyl41xOaOTquJuDO3mrg12LUF\nV7s02KW5EbcWHHFtwTX1tlSqg4Srdo4H46odTeSq3YahN935bqKd/HnVzuYiAIkNBiSaAKTl\nJxh3J2Ql2AuQRAIks1SAdDW9h/tT/l6AJBIgmaUApN1NKKQH+bsBkkiAZJYCkLabz10YxN8N\nkEQCJLMUgBS6gEJ6i78XIIkESGapAGlW1NFvcvl7AZJIgGSWCpBCX91y+i9H7U2wEyCJBEhm\nKQHJVoAkEiCZAZIVIIkESGaAZAVIIgGSGSBZAZJIgETbP+ONzw+4MhmQzABJJJ9BWnye8V4F\nK9wYDUhmgCSSvyDtvzD6CMvFOS7MBiQzQBLJX5A+M588M9eF2YBkBkgi+QvSmyak91yYDUhm\ngCSSvyB9ZUJa7MJsQDIDJJH8BSn3N1FHt7kwGpCsAEkkf0EKbUzTtCZ3bHFjNCCZAZJIPoMU\nCu34dpc7gwHJDJBE8h0kPLPBCpBEAiQzQLICJJEAyQyQrABJJEAyAyQrQBIJkMwAyQqQRAIk\nM0CyAiSRAMkMkKwASSRAMgMkK0ASCZDMAMkKkEQCJDNAsgIkkQDJDJCsAg3pyyl7FJ8ss6I8\nd+YumrLJncGF+e7MXTVltTuD8136H+XHKW780Ydenkv/Ve2c0uDvJtVBstuIzjlH6lup6eXO\n2clegrM+6fx5spfgrKWd30j2Epy1q/NjjX4NkBoLkNwOkEQCJLcDJLcDJJEAye0ASaSywuoj\n9a3UVFEYSfYSnFVZGE72EpxVVViR7CU4q7qwrNGvHTFICAU5QEJIQYCEkIIACSEFuQ6pKqNI\n/5g//t4+Y3YQEpk6IHOit28T0wV/nKbXzUcLzh9/T8a4kB8WbF0WrJV7fsW1Cybr04saWbDL\nkCrXjkszfliPDlm3aWxGHpncb3l21nh3v6dU1oJfeiI7O3s18c+CRw5ftmLUEB8smF0WrJV7\nfsXWggkpHWCsOO6CXYY0vX9f41sfStugQ86YU9ZrMSGruhe4+01lMhdMhtPHZHyz4Mr07wjZ\nkJbv/QVblwVr5Z5fMVswIX9/WF9x/AW7ftVui/HDyv1A/01Y0XP2hrQS/Rd6+mq3v6lM0QWT\njCf73fXEXuKfBY8ct/fA+Ad9sGDrskDMlXt+xbUL/vr+7/UVx1/wkYFkVDG2f9HS7saxjHlu\nf1OZogsuTHtq/dpR/Up9s2BSkJGWdmeI+GHB9LJAzJX7YsXRBedkbDZWHH/BRwpSzfz+IwvI\nkh7RNcx1+5vKFF1w5FANISV3LPDNgssHv7hz96uDiv2wYHpZIObKfbBiuuDqv34UXXH8BR8h\nSAWP3Legxvg1XqZfSNM9/SQ29iuUkD994psFL+4d0c/wzPk+WLB5WSDsqp3XV2wu+NNBu/Yu\nSduYF3/BRwZSzZ+fKjU+Ke25nJB13fLc/qYyRRe8YrDxv3yvb32z4AW9qgipvmeO9xdsXRaI\nuXLPr9ha8KS0aC/FX/CRgbQmfcEavRB5fdDWbUMmuP09paLnbuaY734YMzjimwUXZT67adOL\nd+d5f8HssmD98vf6iussmK447oKPDKRPKeZZJDK5f+YkLz/4Zp27Ox+7897x+cQ/C977bN+M\nJ3b6YMHssmCt3OsrrrNg8wZ0vAXjKUIIKQiQEFIQICGkIEBCSEGAhJCCAAkhBQESQgoCJIQU\nBEgIKQiQglnXy5K9ghQLkPzSC9ohBzsC0hEOkPwSIHk6QPJLFqSylXZ2BKQjHCB5vKJHftqq\n419KyPWapvUlXXvOOq4DIdt7n92myxf6V7t223PzMe0HFupHv7zu+Cve+Pux1o6Xbb+tXfsB\n3n1JkaAFSB6v29F3PPl7LYus+aP22QbS9dITe08ka9qcNuLxi5q8qUO6ussnOyY1uY+Qfx91\n8RODWpx+rLXjaWcMntJD/3foyARI3q6wyUP6x97nW9fYtLf0T6876zAh4euPK9Y//0r/vOtZ\npPKsy8sJ+Vw7lu04mZCaizsmefmpEyB5u6Iml+6lx6iPE6oJydOeNjZ8os0jXdsaxwa0Iwu1\nD41jFzBIxxpvSnNv+yQtO/UCJI/35FFNrxu1jFg+OunHlmlmH5KuvzL2yWpHpmprjWM9GKSL\njE/7AdKRCpC83vox17TQ0iJ17ozL1kYuiHbAvHNOh/QahdTr2Jh77QDpiAVI3q5gYykh+Vna\nzDo+CrVRxpf2LyivhTRP+8g49ktASlKA5O3macartX+ufab7yLV83NhOP1p9U/tILaTin1xV\naewdhZQLSEc8QPJ2Jee0znx+wEnnFJJ/aI98Y/pYfeypox67VHuP1ELSbyRd9sxDJ1x3Eqm7\nIyAdsQDJ423qfVqLDlm7CNl5Q+sHrCcsbOp+xvG/Nl4fin5+/3n6h0+ubHP9//7285gdAemI\nBUjBKHKo3Di464ZkLyRVA6RgVNL8fv1jTutnkr2QVA2QAtIfmgx4/9Vz2uQmex2pGiAFpMqn\nz291Vvq2ZC8jZQMkhBQESAgpCJAQUhAgIaQgQEJIQYCEkIIACSEFARJCCgIkhBQESAgp6P8B\nrULxyAmRBSMAAAAASUVORK5CYII=",
 293       "text/plain": [
 294        "plot without title"
 295       ]
 296      },
 297      "metadata": {},
 298      "output_type": "display_data"
 299     }
 300    ],
 301    "source": [
 302     "hardloss <- ggplot(rubber, aes(x=hardness, y=loss)) + geom_point()\n",
 303     "strloss <-  ggplot(rubber, aes(x=strength, y=loss)) + geom_point()\n",
 304     "\n",
 305     "# plot_grid(hardloss, strloss, labels = \"AUTO\")\n",
 306     "\n",
 307     "multiplot(hardloss, strloss)"
 308    ]
 309   },
 310   {
 311    "cell_type": "markdown",
 312    "metadata": {},
 313    "source": [
 314     "In exercise 3.5(a) you obtained the following output for the regression of abrasion loss on hardness."
 315    ]
 316   },
 317   {
 318    "cell_type": "code",
 319    "execution_count": 6,
 320    "metadata": {},
 321    "outputs": [
 322     {
 323      "data": {
 324       "text/plain": [
 325        "\n",
 326        "Call:\n",
 327        "lm(formula = loss ~ hardness, data = rubber)\n",
 328        "\n",
 329        "Residuals:\n",
 330        "   Min     1Q Median     3Q    Max \n",
 331        "-86.15 -46.77 -19.49  54.27 111.49 \n",
 332        "\n",
 333        "Coefficients:\n",
 334        "            Estimate Std. Error t value Pr(>|t|)    \n",
 335        "(Intercept) 550.4151    65.7867   8.367 4.22e-09 ***\n",
 336        "hardness     -5.3366     0.9229  -5.782 3.29e-06 ***\n",
 337        "---\n",
 338        "Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
 339        "\n",
 340        "Residual standard error: 60.52 on 28 degrees of freedom\n",
 341        "Multiple R-squared:  0.5442,\tAdjusted R-squared:  0.5279 \n",
 342        "F-statistic: 33.43 on 1 and 28 DF,  p-value: 3.294e-06\n"
 343       ]
 344      },
 345      "metadata": {},
 346      "output_type": "display_data"
 347     },
 348     {
 349      "data": {
 350       "text/html": [
 351        "<table>\n",
 352        "<thead><tr><th></th><th scope=col>Df</th><th scope=col>Sum Sq</th><th scope=col>Mean Sq</th><th scope=col>F value</th><th scope=col>Pr(&gt;F)</th></tr></thead>\n",
 353        "<tbody>\n",
 354        "\t<tr><th scope=row>hardness</th><td> 1          </td><td>122455.0    </td><td>122455.037  </td><td>33.43276    </td><td>3.294489e-06</td></tr>\n",
 355        "\t<tr><th scope=row>Residuals</th><td>28          </td><td>102556.3    </td><td>  3662.726  </td><td>      NA    </td><td>          NA</td></tr>\n",
 356        "</tbody>\n",
 357        "</table>\n"
 358       ],
 359       "text/latex": [
 360        "\\begin{tabular}{r|lllll}\n",
 361        "  & Df & Sum Sq & Mean Sq & F value & Pr(>F)\\\\\n",
 362        "\\hline\n",
 363        "\thardness &  1           & 122455.0     & 122455.037   & 33.43276     & 3.294489e-06\\\\\n",
 364        "\tResiduals & 28           & 102556.3     &   3662.726   &       NA     &           NA\\\\\n",
 365        "\\end{tabular}\n"
 366       ],
 367       "text/markdown": [
 368        "\n",
 369        "| <!--/--> | Df | Sum Sq | Mean Sq | F value | Pr(>F) | \n",
 370        "|---|---|\n",
 371        "| hardness |  1           | 122455.0     | 122455.037   | 33.43276     | 3.294489e-06 | \n",
 372        "| Residuals | 28           | 102556.3     |   3662.726   |       NA     |           NA | \n",
 373        "\n",
 374        "\n"
 375       ],
 376       "text/plain": [
 377        "          Df Sum Sq   Mean Sq    F value  Pr(>F)      \n",
 378        "hardness   1 122455.0 122455.037 33.43276 3.294489e-06\n",
 379        "Residuals 28 102556.3   3662.726       NA           NA"
 380       ]
 381      },
 382      "metadata": {},
 383      "output_type": "display_data"
 384     },
 385     {
 386      "name": "stdout",
 387      "output_type": "stream",
 388      "text": [
 389       "          Df   Sum Sq    Mean Sq  F value       Pr(>F)   PctExp\n",
 390       "hardness   1 122455.0 122455.037 33.43276 3.294489e-06 54.42171\n",
 391       "Residuals 28 102556.3   3662.726       NA           NA 45.57829\n"
 392      ]
 393     }
 394    ],
 395    "source": [
 396     "fit <- lm(loss ~ hardness, data = rubber)\n",
 397     "summary(fit)\n",
 398     "anova(fit)\n",
 399     "af <- anova(fit)\n",
 400     "afss <- af$\"Sum Sq\"\n",
 401     "print(cbind(af,PctExp=afss/sum(afss)*100))"
 402    ]
 403   },
 404   {
 405    "cell_type": "markdown",
 406    "metadata": {
 407     "solution": "shown",
 408     "solution2": "hidden",
 409     "solution2_first": true,
 410     "solution_first": true
 411    },
 412    "source": [
 413     "## Exercise\n",
 414     "Now repeat the for the regression of abrasion loss on tensile strength."
 415    ]
 416   },
 417   {
 418    "cell_type": "markdown",
 419    "metadata": {},
 420    "source": [
 421     "### Solution"
 422    ]
 423   },
 424   {
 425    "cell_type": "code",
 426    "execution_count": 15,
 427    "metadata": {
 428     "solution2": "hidden"
 429    },
 430    "outputs": [
 431     {
 432      "data": {
 433       "text/plain": [
 434        "\n",
 435        "Call:\n",
 436        "lm(formula = loss ~ strength, data = rubber)\n",
 437        "\n",
 438        "Residuals:\n",
 439        "     Min       1Q   Median       3Q      Max \n",
 440        "-155.640  -59.919    2.795   61.221  183.285 \n",
 441        "\n",
 442        "Coefficients:\n",
 443        "            Estimate Std. Error t value Pr(>|t|)    \n",
 444        "(Intercept) 305.2248    79.9962   3.815 0.000688 ***\n",
 445        "strength     -0.7192     0.4347  -1.654 0.109232    \n",
 446        "---\n",
 447        "Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
 448        "\n",
 449        "Residual standard error: 85.56 on 28 degrees of freedom\n",
 450        "Multiple R-squared:  0.08904,\tAdjusted R-squared:  0.0565 \n",
 451        "F-statistic: 2.737 on 1 and 28 DF,  p-value: 0.1092\n"
 452       ]
 453      },
 454      "metadata": {},
 455      "output_type": "display_data"
 456     },
 457     {
 458      "data": {
 459       "text/html": [
 460        "<table>\n",
 461        "<thead><tr><th></th><th scope=col>Df</th><th scope=col>Sum Sq</th><th scope=col>Mean Sq</th><th scope=col>F value</th><th scope=col>Pr(&gt;F)</th></tr></thead>\n",
 462        "<tbody>\n",
 463        "\t<tr><th scope=row>strength</th><td> 1       </td><td> 20034.77</td><td>20034.772</td><td>2.736769 </td><td>0.1092317</td></tr>\n",
 464        "\t<tr><th scope=row>Residuals</th><td>28       </td><td>204976.59</td><td> 7320.593</td><td>      NA </td><td>       NA</td></tr>\n",
 465        "</tbody>\n",
 466        "</table>\n"
 467       ],
 468       "text/latex": [
 469        "\\begin{tabular}{r|lllll}\n",
 470        "  & Df & Sum Sq & Mean Sq & F value & Pr(>F)\\\\\n",
 471        "\\hline\n",
 472        "\tstrength &  1        &  20034.77 & 20034.772 & 2.736769  & 0.1092317\\\\\n",
 473        "\tResiduals & 28        & 204976.59 &  7320.593 &       NA  &        NA\\\\\n",
 474        "\\end{tabular}\n"
 475       ],
 476       "text/markdown": [
 477        "\n",
 478        "| <!--/--> | Df | Sum Sq | Mean Sq | F value | Pr(>F) | \n",
 479        "|---|---|\n",
 480        "| strength |  1        |  20034.77 | 20034.772 | 2.736769  | 0.1092317 | \n",
 481        "| Residuals | 28        | 204976.59 |  7320.593 |       NA  |        NA | \n",
 482        "\n",
 483        "\n"
 484       ],
 485       "text/plain": [
 486        "          Df Sum Sq    Mean Sq   F value  Pr(>F)   \n",
 487        "strength   1  20034.77 20034.772 2.736769 0.1092317\n",
 488        "Residuals 28 204976.59  7320.593       NA        NA"
 489       ]
 490      },
 491      "metadata": {},
 492      "output_type": "display_data"
 493     },
 494     {
 495      "name": "stdout",
 496      "output_type": "stream",
 497      "text": [
 498       "          Df    Sum Sq   Mean Sq  F value    Pr(>F)    PctExp\n",
 499       "strength   1  20034.77 20034.772 2.736769 0.1092317  8.903893\n",
 500       "Residuals 28 204976.59  7320.593       NA        NA 91.096107\n"
 501      ]
 502     }
 503    ],
 504    "source": [
 505     "fit <- lm(loss ~ strength, data = rubber)\n",
 506     "summary(fit)\n",
 507     "anova(fit)\n",
 508     "af <- anova(fit)\n",
 509     "afss <- af$\"Sum Sq\"\n",
 510     "print(cbind(af,PctExp=afss/sum(afss)*100))"
 511    ]
 512   },
 513   {
 514    "cell_type": "markdown",
 515    "metadata": {
 516     "solution2": "hidden"
 517    },
 518    "source": [
 519     "Individually, then, regression of abrasion loss on hardness seems satisfactory, albeit accounting for a disappointingly small percentage of the variance in the response; regression of abrasion loss on tensile strength,  on the other hand, seems to have very little explanatory power. "
 520    ]
 521   },
 522   {
 523    "cell_type": "markdown",
 524    "metadata": {
 525     "solution": "hidden"
 526    },
 527    "source": [
 528     "However, back in Exercise 3.19 it was shown that the residuals from the former regression showed a clear relationship with the tensile strength values, and this suggested that a regression model involving both variables is necessary. \n",
 529     "\n",
 530     "Let us try it. The only change is to include `strength` in the equations being fitted in the `lm()` function."
 531    ]
 532   },
 533   {
 534    "cell_type": "code",
 535    "execution_count": 17,
 536    "metadata": {},
 537    "outputs": [
 538     {
 539      "data": {
 540       "text/plain": [
 541        "\n",
 542        "Call:\n",
 543        "lm(formula = loss ~ hardness + strength, data = rubber)\n",
 544        "\n",
 545        "Residuals:\n",
 546        "    Min      1Q  Median      3Q     Max \n",
 547        "-79.385 -14.608   3.816  19.755  65.981 \n",
 548        "\n",
 549        "Coefficients:\n",
 550        "            Estimate Std. Error t value Pr(>|t|)    \n",
 551        "(Intercept) 885.1611    61.7516  14.334 3.84e-14 ***\n",
 552        "hardness     -6.5708     0.5832 -11.267 1.03e-11 ***\n",
 553        "strength     -1.3743     0.1943  -7.073 1.32e-07 ***\n",
 554        "---\n",
 555        "Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
 556        "\n",
 557        "Residual standard error: 36.49 on 27 degrees of freedom\n",
 558        "Multiple R-squared:  0.8402,\tAdjusted R-squared:  0.8284 \n",
 559        "F-statistic:    71 on 2 and 27 DF,  p-value: 1.767e-11\n"
 560       ]
 561      },
 562      "metadata": {},
 563      "output_type": "display_data"
 564     },
 565     {
 566      "data": {
 567       "text/html": [
 568        "<table>\n",
 569        "<thead><tr><th></th><th scope=col>Df</th><th scope=col>Sum Sq</th><th scope=col>Mean Sq</th><th scope=col>F value</th><th scope=col>Pr(&gt;F)</th></tr></thead>\n",
 570        "<tbody>\n",
 571        "\t<tr><th scope=row>hardness</th><td> 1          </td><td>122455.04   </td><td>122455.037  </td><td>91.96967    </td><td>3.458255e-10</td></tr>\n",
 572        "\t<tr><th scope=row>strength</th><td> 1          </td><td> 66606.59   </td><td> 66606.586  </td><td>50.02477    </td><td>1.324645e-07</td></tr>\n",
 573        "\t<tr><th scope=row>Residuals</th><td>27          </td><td> 35949.74   </td><td>  1331.472  </td><td>      NA    </td><td>          NA</td></tr>\n",
 574        "</tbody>\n",
 575        "</table>\n"
 576       ],
 577       "text/latex": [
 578        "\\begin{tabular}{r|lllll}\n",
 579        "  & Df & Sum Sq & Mean Sq & F value & Pr(>F)\\\\\n",
 580        "\\hline\n",
 581        "\thardness &  1           & 122455.04    & 122455.037   & 91.96967     & 3.458255e-10\\\\\n",
 582        "\tstrength &  1           &  66606.59    &  66606.586   & 50.02477     & 1.324645e-07\\\\\n",
 583        "\tResiduals & 27           &  35949.74    &   1331.472   &       NA     &           NA\\\\\n",
 584        "\\end{tabular}\n"
 585       ],
 586       "text/markdown": [
 587        "\n",
 588        "| <!--/--> | Df | Sum Sq | Mean Sq | F value | Pr(>F) | \n",
 589        "|---|---|---|\n",
 590        "| hardness |  1           | 122455.04    | 122455.037   | 91.96967     | 3.458255e-10 | \n",
 591        "| strength |  1           |  66606.59    |  66606.586   | 50.02477     | 1.324645e-07 | \n",
 592        "| Residuals | 27           |  35949.74    |   1331.472   |       NA     |           NA | \n",
 593        "\n",
 594        "\n"
 595       ],
 596       "text/plain": [
 597        "          Df Sum Sq    Mean Sq    F value  Pr(>F)      \n",
 598        "hardness   1 122455.04 122455.037 91.96967 3.458255e-10\n",
 599        "strength   1  66606.59  66606.586 50.02477 1.324645e-07\n",
 600        "Residuals 27  35949.74   1331.472       NA           NA"
 601       ]
 602      },
 603      "metadata": {},
 604      "output_type": "display_data"
 605     },
 606     {
 607      "name": "stdout",
 608      "output_type": "stream",
 609      "text": [
 610       "          Df    Sum Sq    Mean Sq  F value       Pr(>F)   PctExp\n",
 611       "hardness   1 122455.04 122455.037 91.96967 3.458255e-10 54.42171\n",
 612       "strength   1  66606.59  66606.586 50.02477 1.324645e-07 29.60143\n",
 613       "Residuals 27  35949.74   1331.472       NA           NA 15.97686\n"
 614      ]
 615     }
 616    ],
 617    "source": [
 618     "fit <- lm(loss ~ hardness + strength, data = rubber)\n",
 619     "summary(fit)\n",
 620     "anova(fit)\n",
 621     "af <- anova(fit)\n",
 622     "afss <- af$\"Sum Sq\"\n",
 623     "print(cbind(af,PctExp=afss/sum(afss)*100))"
 624    ]
 625   },
 626   {
 627    "cell_type": "markdown",
 628    "metadata": {},
 629    "source": [
 630     "Looking at the regression coefficient output next, the estimated model for the mean response is\n",
 631     "\n",
 632     "$$ \\hat{y} = \\hat{\\alpha} + \\hat{\\beta}_1 x_1 + \\hat{\\beta}_2 x_2 = 885.2 − 6.571 x_1 − 1.374 x_2 $$ \n",
 633     "\n",
 634     "where $x_1$ and $x_2$ stand for values of hardness and tensile strength, respectively. \n",
 635     "\n",
 636     "Associated with each estimated parameter, GenStat gives a standard error (details of the calculation of which need not concern us now) and hence a _t_-statistic (estimate divided by standard error) to be compared with the distribution on `d.f. (Residual)`=27 degrees of freedom. GenStat makes the comparison and gives _p_ values, which in this case are all very small, suggesting strong evidence for the non-zeroness (and hence presence) of each parameter, $\\alpha$, $\\beta_1$ and $\\beta_2$, individually."
 637    ]
 638   },
 639   {
 640    "cell_type": "markdown",
 641    "metadata": {},
 642    "source": [
 643     "So, even though the single regression of abrasion loss on tensile strength did not suggest a close relationship, when taken in conjunction with the effect of hardness the effect of tensile strength is also a considerable one. A key idea here is that $\\beta_1$ and $\\beta_2$ (and more generally $\\beta_1, \\beta_2, \\ldots , \\beta_k$ in model (5.1)) are partial regression coefficients. That is, $\\beta_1$ measures the effect of an increase of one unit in $x_1$ treating the value of the other variable $x_2$ as fixed. Contrast this with the single regression model $E(Y) = \\alpha + \\beta_1 x_1$ in which $\\beta_1$ represents an increase of one unit in $x_1$ treating $x_2$ as zero. The meaning of $\\beta_1$ in the regression models with one and two explanatory variables is not the same.\n",
 644     "\n",
 645     "You will notice that the percentage of variance accounted for has increased dramatically in the two-explanatory-variable model to a much more respectable 82.8%. This statistic has the same interpretation — or difficulty of interpretation — as for the case of one explanatory variable (see [Unit 3](unit3.ipynb)).\n",
 646     "\n",
 647     "Other items in the GenStat output are as for the case of one explanatory variable: the `Standard error of observations` is $\\hat{\\sigma}$, and is the square root of `m.s. (Residual)`; the message about standardised residuals will be clarified below; and the message about leverage will be ignored until [Unit 10](unit10.ipynb).\n",
 648     "\n",
 649     "The simple residuals are, again, defined as the differences between the observed and predicted responses:\n",
 650     "\n",
 651     "$$ r_i = y_i - \\left( \\hat{\\alpha} + \\sum_{j=1}^{k} \\hat{\\beta}_j x_{i, j} \\right) , i = 1, 2, \\ldots, n. $$\n",
 652     "\n",
 653     "GenStat can obtain these for you and produce a plot of residuals against fitted values. The fitted values are simply $\\hat{Y}_i = \\hat{\\alpha} + \\sum_{j=1}^{k} \\hat{\\beta}_j x_{i, j}$. As in the regression models in [Unit 3](unit3.ipynb), the default in GenStat is to use deviance residuals which, in this context, are equivalent to standardised residuals (simple residuals divided by their estimated standard errors)."
 654    ]
 655   },
 656   {
 657    "cell_type": "markdown",
 658    "metadata": {},
 659    "source": [
 660     "## Example 5.2"
 661    ]
 662   },
 663   {
 664    "cell_type": "code",
 665    "execution_count": null,
 666    "metadata": {},
 667    "outputs": [],
 668    "source": []
 669   }
 670  ],
 671  "metadata": {
 672   "kernelspec": {
 673    "display_name": "R",
 674    "language": "R",
 675    "name": "ir"
 676   },
 677   "language_info": {
 678    "codemirror_mode": "r",
 679    "file_extension": ".r",
 680    "mimetype": "text/x-r-source",
 681    "name": "R",
 682    "pygments_lexer": "r",
 683    "version": "3.4.2"
 684   },
 685   "toc": {
 686    "nav_menu": {},
 687    "number_sections": true,
 688    "sideBar": true,
 689    "skip_h1_title": false,
 690    "title_cell": "Table of Contents",
 691    "title_sidebar": "Contents",
 692    "toc_cell": false,
 693    "toc_position": {
 694     "height": "calc(100% - 180px)",
 695     "left": "10px",
 696     "top": "150px",
 697     "width": "342px"
 698    },
 699    "toc_section_display": true,
 700    "toc_window_display": true
 701   }
 702  },
 703  "nbformat": 4,
 704  "nbformat_minor": 2
 705 }