+++ /dev/null
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Section 5.1: Using the model"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Imports and defintions"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 1,
- "metadata": {},
- "outputs": [
- {
- "name": "stderr",
- "output_type": "stream",
- "text": [
- "── Attaching packages ─────────────────────────────────────── tidyverse 1.2.1 ──\n",
- "✔ ggplot2 2.2.1 ✔ purrr 0.2.4\n",
- "✔ tibble 1.4.2 ✔ dplyr 0.7.4\n",
- "✔ tidyr 0.8.0 ✔ stringr 1.2.0\n",
- "✔ readr 1.1.1 ✔ forcats 0.2.0\n",
- "── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──\n",
- "✖ dplyr::filter() masks stats::filter()\n",
- "✖ dplyr::lag() masks stats::lag()\n"
- ]
- }
- ],
- "source": [
- "library(tidyverse)\n",
- "# library(cowplot)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 2,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Multiple plot function\n",
- "#\n",
- "# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)\n",
- "# - cols: Number of columns in layout\n",
- "# - layout: A matrix specifying the layout. If present, 'cols' is ignored.\n",
- "#\n",
- "# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),\n",
- "# then plot 1 will go in the upper left, 2 will go in the upper right, and\n",
- "# 3 will go all the way across the bottom.\n",
- "#\n",
- "multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {\n",
- " library(grid)\n",
- "\n",
- " # Make a list from the ... arguments and plotlist\n",
- " plots <- c(list(...), plotlist)\n",
- "\n",
- " numPlots = length(plots)\n",
- "\n",
- " # If layout is NULL, then use 'cols' to determine layout\n",
- " if (is.null(layout)) {\n",
- " # Make the panel\n",
- " # ncol: Number of columns of plots\n",
- " # nrow: Number of rows needed, calculated from # of cols\n",
- " layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),\n",
- " ncol = cols, nrow = ceiling(numPlots/cols))\n",
- " }\n",
- "\n",
- " if (numPlots==1) {\n",
- " print(plots[[1]])\n",
- "\n",
- " } else {\n",
- " # Set up the page\n",
- " grid.newpage()\n",
- " pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))\n",
- "\n",
- " # Make each plot, in the correct location\n",
- " for (i in 1:numPlots) {\n",
- " # Get the i,j matrix positions of the regions that contain this subplot\n",
- " matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))\n",
- "\n",
- " print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,\n",
- " layout.pos.col = matchidx$col))\n",
- " }\n",
- " }\n",
- "}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Modelling abrasion loss\n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 3,
- "metadata": {
- "scrolled": true
- },
- "outputs": [
- {
- "data": {
- "text/html": [
- "<table>\n",
- "<thead><tr><th scope=col>loss</th><th scope=col>hardness</th><th scope=col>strength</th></tr></thead>\n",
- "<tbody>\n",
- "\t<tr><td>372</td><td>45 </td><td>162</td></tr>\n",
- "\t<tr><td>206</td><td>55 </td><td>233</td></tr>\n",
- "\t<tr><td>175</td><td>61 </td><td>232</td></tr>\n",
- "\t<tr><td>154</td><td>66 </td><td>231</td></tr>\n",
- "\t<tr><td>136</td><td>71 </td><td>231</td></tr>\n",
- "\t<tr><td>112</td><td>71 </td><td>237</td></tr>\n",
- "\t<tr><td> 55</td><td>81 </td><td>224</td></tr>\n",
- "\t<tr><td> 45</td><td>86 </td><td>219</td></tr>\n",
- "\t<tr><td>221</td><td>53 </td><td>203</td></tr>\n",
- "\t<tr><td>166</td><td>60 </td><td>189</td></tr>\n",
- "\t<tr><td>164</td><td>64 </td><td>210</td></tr>\n",
- "\t<tr><td>113</td><td>68 </td><td>210</td></tr>\n",
- "\t<tr><td> 82</td><td>79 </td><td>196</td></tr>\n",
- "\t<tr><td> 32</td><td>81 </td><td>180</td></tr>\n",
- "\t<tr><td>228</td><td>56 </td><td>200</td></tr>\n",
- "\t<tr><td>196</td><td>68 </td><td>173</td></tr>\n",
- "\t<tr><td>128</td><td>75 </td><td>188</td></tr>\n",
- "\t<tr><td> 97</td><td>83 </td><td>161</td></tr>\n",
- "\t<tr><td> 64</td><td>88 </td><td>119</td></tr>\n",
- "\t<tr><td>249</td><td>59 </td><td>161</td></tr>\n",
- "\t<tr><td>219</td><td>71 </td><td>151</td></tr>\n",
- "\t<tr><td>186</td><td>80 </td><td>165</td></tr>\n",
- "\t<tr><td>155</td><td>82 </td><td>151</td></tr>\n",
- "\t<tr><td>114</td><td>89 </td><td>128</td></tr>\n",
- "\t<tr><td>341</td><td>51 </td><td>161</td></tr>\n",
- "\t<tr><td>340</td><td>59 </td><td>146</td></tr>\n",
- "\t<tr><td>283</td><td>65 </td><td>148</td></tr>\n",
- "\t<tr><td>267</td><td>74 </td><td>144</td></tr>\n",
- "\t<tr><td>215</td><td>81 </td><td>134</td></tr>\n",
- "\t<tr><td>148</td><td>86 </td><td>127</td></tr>\n",
- "</tbody>\n",
- "</table>\n"
- ],
- "text/latex": [
- "\\begin{tabular}{r|lll}\n",
- " loss & hardness & strength\\\\\n",
- "\\hline\n",
- "\t 372 & 45 & 162\\\\\n",
- "\t 206 & 55 & 233\\\\\n",
- "\t 175 & 61 & 232\\\\\n",
- "\t 154 & 66 & 231\\\\\n",
- "\t 136 & 71 & 231\\\\\n",
- "\t 112 & 71 & 237\\\\\n",
- "\t 55 & 81 & 224\\\\\n",
- "\t 45 & 86 & 219\\\\\n",
- "\t 221 & 53 & 203\\\\\n",
- "\t 166 & 60 & 189\\\\\n",
- "\t 164 & 64 & 210\\\\\n",
- "\t 113 & 68 & 210\\\\\n",
- "\t 82 & 79 & 196\\\\\n",
- "\t 32 & 81 & 180\\\\\n",
- "\t 228 & 56 & 200\\\\\n",
- "\t 196 & 68 & 173\\\\\n",
- "\t 128 & 75 & 188\\\\\n",
- "\t 97 & 83 & 161\\\\\n",
- "\t 64 & 88 & 119\\\\\n",
- "\t 249 & 59 & 161\\\\\n",
- "\t 219 & 71 & 151\\\\\n",
- "\t 186 & 80 & 165\\\\\n",
- "\t 155 & 82 & 151\\\\\n",
- "\t 114 & 89 & 128\\\\\n",
- "\t 341 & 51 & 161\\\\\n",
- "\t 340 & 59 & 146\\\\\n",
- "\t 283 & 65 & 148\\\\\n",
- "\t 267 & 74 & 144\\\\\n",
- "\t 215 & 81 & 134\\\\\n",
- "\t 148 & 86 & 127\\\\\n",
- "\\end{tabular}\n"
- ],
- "text/markdown": [
- "\n",
- "loss | hardness | strength | \n",
- "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
- "| 372 | 45 | 162 | \n",
- "| 206 | 55 | 233 | \n",
- "| 175 | 61 | 232 | \n",
- "| 154 | 66 | 231 | \n",
- "| 136 | 71 | 231 | \n",
- "| 112 | 71 | 237 | \n",
- "| 55 | 81 | 224 | \n",
- "| 45 | 86 | 219 | \n",
- "| 221 | 53 | 203 | \n",
- "| 166 | 60 | 189 | \n",
- "| 164 | 64 | 210 | \n",
- "| 113 | 68 | 210 | \n",
- "| 82 | 79 | 196 | \n",
- "| 32 | 81 | 180 | \n",
- "| 228 | 56 | 200 | \n",
- "| 196 | 68 | 173 | \n",
- "| 128 | 75 | 188 | \n",
- "| 97 | 83 | 161 | \n",
- "| 64 | 88 | 119 | \n",
- "| 249 | 59 | 161 | \n",
- "| 219 | 71 | 151 | \n",
- "| 186 | 80 | 165 | \n",
- "| 155 | 82 | 151 | \n",
- "| 114 | 89 | 128 | \n",
- "| 341 | 51 | 161 | \n",
- "| 340 | 59 | 146 | \n",
- "| 283 | 65 | 148 | \n",
- "| 267 | 74 | 144 | \n",
- "| 215 | 81 | 134 | \n",
- "| 148 | 86 | 127 | \n",
- "\n",
- "\n"
- ],
- "text/plain": [
- " loss hardness strength\n",
- "1 372 45 162 \n",
- "2 206 55 233 \n",
- "3 175 61 232 \n",
- "4 154 66 231 \n",
- "5 136 71 231 \n",
- "6 112 71 237 \n",
- "7 55 81 224 \n",
- "8 45 86 219 \n",
- "9 221 53 203 \n",
- "10 166 60 189 \n",
- "11 164 64 210 \n",
- "12 113 68 210 \n",
- "13 82 79 196 \n",
- "14 32 81 180 \n",
- "15 228 56 200 \n",
- "16 196 68 173 \n",
- "17 128 75 188 \n",
- "18 97 83 161 \n",
- "19 64 88 119 \n",
- "20 249 59 161 \n",
- "21 219 71 151 \n",
- "22 186 80 165 \n",
- "23 155 82 151 \n",
- "24 114 89 128 \n",
- "25 341 51 161 \n",
- "26 340 59 146 \n",
- "27 283 65 148 \n",
- "28 267 74 144 \n",
- "29 215 81 134 \n",
- "30 148 86 127 "
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "# Load the data from the file.\n",
- "rubber <- read.csv('rubber.csv')\n",
- "rubber"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "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."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 4,
- "metadata": {},
- "outputs": [],
- "source": [
- "# ggplot(rubber, aes(x=hardness, y=loss)) + geom_point()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 5,
- "metadata": {
- "scrolled": false
- },
- "outputs": [
- {
- "data": {
- "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=",
- "text/plain": [
- "plot without title"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "hardloss <- ggplot(rubber, aes(x=hardness, y=loss)) + geom_point()\n",
- "strloss <- ggplot(rubber, aes(x=strength, y=loss)) + geom_point()\n",
- "\n",
- "# plot_grid(hardloss, strloss, labels = \"AUTO\")\n",
- "\n",
- "multiplot(hardloss, strloss)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "In exercise 3.5(a) you obtained the following output for the regression of abrasion loss on hardness."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 6,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "\n",
- "Call:\n",
- "lm(formula = loss ~ hardness, data = rubber)\n",
- "\n",
- "Residuals:\n",
- " Min 1Q Median 3Q Max \n",
- "-86.15 -46.77 -19.49 54.27 111.49 \n",
- "\n",
- "Coefficients:\n",
- " Estimate Std. Error t value Pr(>|t|) \n",
- "(Intercept) 550.4151 65.7867 8.367 4.22e-09 ***\n",
- "hardness -5.3366 0.9229 -5.782 3.29e-06 ***\n",
- "---\n",
- "Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
- "\n",
- "Residual standard error: 60.52 on 28 degrees of freedom\n",
- "Multiple R-squared: 0.5442,\tAdjusted R-squared: 0.5279 \n",
- "F-statistic: 33.43 on 1 and 28 DF, p-value: 3.294e-06\n"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- },
- {
- "data": {
- "text/html": [
- "<table>\n",
- "<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(>F)</th></tr></thead>\n",
- "<tbody>\n",
- "\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",
- "\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",
- "</tbody>\n",
- "</table>\n"
- ],
- "text/latex": [
- "\\begin{tabular}{r|lllll}\n",
- " & Df & Sum Sq & Mean Sq & F value & Pr(>F)\\\\\n",
- "\\hline\n",
- "\thardness & 1 & 122455.0 & 122455.037 & 33.43276 & 3.294489e-06\\\\\n",
- "\tResiduals & 28 & 102556.3 & 3662.726 & NA & NA\\\\\n",
- "\\end{tabular}\n"
- ],
- "text/markdown": [
- "\n",
- "| <!--/--> | Df | Sum Sq | Mean Sq | F value | Pr(>F) | \n",
- "|---|---|\n",
- "| hardness | 1 | 122455.0 | 122455.037 | 33.43276 | 3.294489e-06 | \n",
- "| Residuals | 28 | 102556.3 | 3662.726 | NA | NA | \n",
- "\n",
- "\n"
- ],
- "text/plain": [
- " Df Sum Sq Mean Sq F value Pr(>F) \n",
- "hardness 1 122455.0 122455.037 33.43276 3.294489e-06\n",
- "Residuals 28 102556.3 3662.726 NA NA"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- },
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- " Df Sum Sq Mean Sq F value Pr(>F) PctExp\n",
- "hardness 1 122455.0 122455.037 33.43276 3.294489e-06 54.42171\n",
- "Residuals 28 102556.3 3662.726 NA NA 45.57829\n"
- ]
- }
- ],
- "source": [
- "fit <- lm(loss ~ hardness, data = rubber)\n",
- "summary(fit)\n",
- "anova(fit)\n",
- "af <- anova(fit)\n",
- "afss <- af$\"Sum Sq\"\n",
- "print(cbind(af,PctExp=afss/sum(afss)*100))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "solution": "shown",
- "solution2": "hidden",
- "solution2_first": true,
- "solution_first": true
- },
- "source": [
- "## Exercise\n",
- "Now repeat the for the regression of abrasion loss on tensile strength."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Solution"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 15,
- "metadata": {
- "solution2": "hidden"
- },
- "outputs": [
- {
- "data": {
- "text/plain": [
- "\n",
- "Call:\n",
- "lm(formula = loss ~ strength, data = rubber)\n",
- "\n",
- "Residuals:\n",
- " Min 1Q Median 3Q Max \n",
- "-155.640 -59.919 2.795 61.221 183.285 \n",
- "\n",
- "Coefficients:\n",
- " Estimate Std. Error t value Pr(>|t|) \n",
- "(Intercept) 305.2248 79.9962 3.815 0.000688 ***\n",
- "strength -0.7192 0.4347 -1.654 0.109232 \n",
- "---\n",
- "Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
- "\n",
- "Residual standard error: 85.56 on 28 degrees of freedom\n",
- "Multiple R-squared: 0.08904,\tAdjusted R-squared: 0.0565 \n",
- "F-statistic: 2.737 on 1 and 28 DF, p-value: 0.1092\n"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- },
- {
- "data": {
- "text/html": [
- "<table>\n",
- "<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(>F)</th></tr></thead>\n",
- "<tbody>\n",
- "\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",
- "\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",
- "</tbody>\n",
- "</table>\n"
- ],
- "text/latex": [
- "\\begin{tabular}{r|lllll}\n",
- " & Df & Sum Sq & Mean Sq & F value & Pr(>F)\\\\\n",
- "\\hline\n",
- "\tstrength & 1 & 20034.77 & 20034.772 & 2.736769 & 0.1092317\\\\\n",
- "\tResiduals & 28 & 204976.59 & 7320.593 & NA & NA\\\\\n",
- "\\end{tabular}\n"
- ],
- "text/markdown": [
- "\n",
- "| <!--/--> | Df | Sum Sq | Mean Sq | F value | Pr(>F) | \n",
- "|---|---|\n",
- "| strength | 1 | 20034.77 | 20034.772 | 2.736769 | 0.1092317 | \n",
- "| Residuals | 28 | 204976.59 | 7320.593 | NA | NA | \n",
- "\n",
- "\n"
- ],
- "text/plain": [
- " Df Sum Sq Mean Sq F value Pr(>F) \n",
- "strength 1 20034.77 20034.772 2.736769 0.1092317\n",
- "Residuals 28 204976.59 7320.593 NA NA"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- },
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- " Df Sum Sq Mean Sq F value Pr(>F) PctExp\n",
- "strength 1 20034.77 20034.772 2.736769 0.1092317 8.903893\n",
- "Residuals 28 204976.59 7320.593 NA NA 91.096107\n"
- ]
- }
- ],
- "source": [
- "fit <- lm(loss ~ strength, data = rubber)\n",
- "summary(fit)\n",
- "anova(fit)\n",
- "af <- anova(fit)\n",
- "afss <- af$\"Sum Sq\"\n",
- "print(cbind(af,PctExp=afss/sum(afss)*100))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "solution2": "hidden"
- },
- "source": [
- "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. "
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "solution": "hidden"
- },
- "source": [
- "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",
- "\n",
- "Let us try it. The only change is to include `strength` in the equations being fitted in the `lm()` function."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 17,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "\n",
- "Call:\n",
- "lm(formula = loss ~ hardness + strength, data = rubber)\n",
- "\n",
- "Residuals:\n",
- " Min 1Q Median 3Q Max \n",
- "-79.385 -14.608 3.816 19.755 65.981 \n",
- "\n",
- "Coefficients:\n",
- " Estimate Std. Error t value Pr(>|t|) \n",
- "(Intercept) 885.1611 61.7516 14.334 3.84e-14 ***\n",
- "hardness -6.5708 0.5832 -11.267 1.03e-11 ***\n",
- "strength -1.3743 0.1943 -7.073 1.32e-07 ***\n",
- "---\n",
- "Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
- "\n",
- "Residual standard error: 36.49 on 27 degrees of freedom\n",
- "Multiple R-squared: 0.8402,\tAdjusted R-squared: 0.8284 \n",
- "F-statistic: 71 on 2 and 27 DF, p-value: 1.767e-11\n"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- },
- {
- "data": {
- "text/html": [
- "<table>\n",
- "<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(>F)</th></tr></thead>\n",
- "<tbody>\n",
- "\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",
- "\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",
- "\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",
- "</tbody>\n",
- "</table>\n"
- ],
- "text/latex": [
- "\\begin{tabular}{r|lllll}\n",
- " & Df & Sum Sq & Mean Sq & F value & Pr(>F)\\\\\n",
- "\\hline\n",
- "\thardness & 1 & 122455.04 & 122455.037 & 91.96967 & 3.458255e-10\\\\\n",
- "\tstrength & 1 & 66606.59 & 66606.586 & 50.02477 & 1.324645e-07\\\\\n",
- "\tResiduals & 27 & 35949.74 & 1331.472 & NA & NA\\\\\n",
- "\\end{tabular}\n"
- ],
- "text/markdown": [
- "\n",
- "| <!--/--> | Df | Sum Sq | Mean Sq | F value | Pr(>F) | \n",
- "|---|---|---|\n",
- "| hardness | 1 | 122455.04 | 122455.037 | 91.96967 | 3.458255e-10 | \n",
- "| strength | 1 | 66606.59 | 66606.586 | 50.02477 | 1.324645e-07 | \n",
- "| Residuals | 27 | 35949.74 | 1331.472 | NA | NA | \n",
- "\n",
- "\n"
- ],
- "text/plain": [
- " Df Sum Sq Mean Sq F value Pr(>F) \n",
- "hardness 1 122455.04 122455.037 91.96967 3.458255e-10\n",
- "strength 1 66606.59 66606.586 50.02477 1.324645e-07\n",
- "Residuals 27 35949.74 1331.472 NA NA"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- },
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- " Df Sum Sq Mean Sq F value Pr(>F) PctExp\n",
- "hardness 1 122455.04 122455.037 91.96967 3.458255e-10 54.42171\n",
- "strength 1 66606.59 66606.586 50.02477 1.324645e-07 29.60143\n",
- "Residuals 27 35949.74 1331.472 NA NA 15.97686\n"
- ]
- }
- ],
- "source": [
- "fit <- lm(loss ~ hardness + strength, data = rubber)\n",
- "summary(fit)\n",
- "anova(fit)\n",
- "af <- anova(fit)\n",
- "afss <- af$\"Sum Sq\"\n",
- "print(cbind(af,PctExp=afss/sum(afss)*100))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Looking at the regression coefficient output next, the estimated model for the mean response is\n",
- "\n",
- "$$ \\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",
- "\n",
- "where $x_1$ and $x_2$ stand for values of hardness and tensile strength, respectively. \n",
- "\n",
- "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."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "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",
- "\n",
- "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",
- "\n",
- "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",
- "\n",
- "The simple residuals are, again, defined as the differences between the observed and predicted responses:\n",
- "\n",
- "$$ r_i = y_i - \\left( \\hat{\\alpha} + \\sum_{j=1}^{k} \\hat{\\beta}_j x_{i, j} \\right) , i = 1, 2, \\ldots, n. $$\n",
- "\n",
- "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)."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Example 5.2"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "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": "3.4.2"
- },
- "toc": {
- "nav_menu": {},
- "number_sections": true,
- "sideBar": true,
- "skip_h1_title": false,
- "title_cell": "Table of Contents",
- "title_sidebar": "Contents",
- "toc_cell": false,
- "toc_position": {
- "height": "calc(100% - 180px)",
- "left": "10px",
- "top": "150px",
- "width": "342px"
- },
- "toc_section_display": true,
- "toc_window_display": true
- }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
projects / ou-jupyter-r-demo.git / blobdiff