forked from STAT545-UBC/STAT545-UBC-original-website
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathblock012_function-regress-lifeexp-on-year.html
More file actions
244 lines (220 loc) · 13.3 KB
/
Copy pathblock012_function-regress-lifeexp-on-year.html
File metadata and controls
244 lines (220 loc) · 13.3 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta charset="utf-8">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />
<title>Linear regression of life expectancy on year</title>
<script src="libs/jquery-1.11.0/jquery.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1.0" />
<link href="libs/bootstrap-2.3.2/css/united.min.css" rel="stylesheet" />
<link href="libs/bootstrap-2.3.2/css/bootstrap-responsive.min.css" rel="stylesheet" />
<script src="libs/bootstrap-2.3.2/js/bootstrap.min.js"></script>
<style type="text/css">code{white-space: pre;}</style>
<link rel="stylesheet"
href="libs/highlight/default.css"
type="text/css" />
<script src="libs/highlight/highlight.js"></script>
<style type="text/css">
pre:not([class]) {
background-color: white;
}
</style>
<script type="text/javascript">
if (window.hljs && document.readyState && document.readyState === "complete") {
window.setTimeout(function() {
hljs.initHighlighting();
}, 0);
}
</script>
<link rel="stylesheet" href="libs/local/nav.css" type="text/css" />
</head>
<body>
<style type = "text/css">
.main-container {
max-width: 940px;
margin-left: auto;
margin-right: auto;
}
</style>
<div class="container-fluid main-container">
<header>
<div class="nav">
<a class="nav-logo" href="index.html">
<img src="static/img/stat545-logo-s.png" width="70px" height="70px"/>
</a>
<ul>
<li class="home"><a href="index.html">Home</a></li>
<li class="faq"><a href="faq.html">FAQ</a></li>
<li class="syllabus"><a href="syllabus.html">Syllabus</a></li>
<li class="topics"><a href="topics.html">Topics</a></li>
<li class="people"><a href="people.html">People</a></li>
</ul>
</div>
</header>
<div id="header">
<h1 class="title">Linear regression of life expectancy on year</h1>
</div>
<div id="TOC">
<ul>
<li><a href="#overview">Overview</a></li>
<li><a href="#load-the-gapminder-data">Load the Gapminder data</a></li>
<li><a href="#get-data-to-practice-with">Get data to practice with</a></li>
<li><a href="#get-some-code-that-works">Get some code that works</a></li>
<li><a href="#turn-working-code-into-a-function">Turn working code into a function</a></li>
<li><a href="#test-on-other-data-and-in-a-clean-workspace">Test on other data and in a clean workspace</a></li>
<li><a href="#are-we-there-yet">Are we there yet?</a></li>
</ul>
</div>
<div id="overview" class="section level3">
<h3>Overview</h3>
<p>We recently learned how to write our own R functions (<a href="block011_write-your-own-function-01.html">Part 1</a>, <a href="block011_write-your-own-function-02.html">Part 2</a>, <a href="block011_write-your-own-function-03.html">Part 3</a>).</p>
<p>Now we use that knowledge to write another useful function, within the context of the Gapminder data:</p>
<ul>
<li>Input: a data.frame that contains (at least) a life expectancy variable <code>lifeExp</code> and a variable for year <code>year</code></li>
<li>Output: a vector of estimated intercept and slope, from a linear regression of <code>lifeExp</code> on <code>year</code></li>
</ul>
<p>The ultimate goal is to apply this function to the Gapminder data for a specific country. We will eventually scale up to <em>all</em> countries using external machinery, e.g., the <code>plyr</code> package.</p>
</div>
<div id="load-the-gapminder-data" class="section level3">
<h3>Load the Gapminder data</h3>
<p>As usual, load the Gapminder excerpt. Load <code>ggplot2</code> because we’ll make some plots.</p>
<pre class="r"><code>library(ggplot2)
gDat <- read.delim("gapminderDataFiveYear.txt")
str(gDat)
## 'data.frame': 1704 obs. of 6 variables:
## $ country : Factor w/ 142 levels "Afghanistan",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ year : int 1952 1957 1962 1967 1972 1977 1982 1987 1992 1997 ...
## $ pop : num 8425333 9240934 10267083 11537966 13079460 ...
## $ continent: Factor w/ 5 levels "Africa","Americas",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ lifeExp : num 28.8 30.3 32 34 36.1 ...
## $ gdpPercap: num 779 821 853 836 740 ...
## or do this if the file isn't lying around already
## gd_url <- "http://tiny.cc/gapminder"
## gDat <- read.delim(gd_url)</code></pre>
</div>
<div id="get-data-to-practice-with" class="section level3">
<h3>Get data to practice with</h3>
<p>I extract the data for one country in order to develop some working code interactively.</p>
<pre class="r"><code>j_country <- "France" # pick, but do not hard wire, an example
(j_dat <- subset(gDat, country == j_country))
## country year pop continent lifeExp gdpPercap
## 529 France 1952 42459667 Europe 67.410 7029.809
## 530 France 1957 44310863 Europe 68.930 8662.835
## 531 France 1962 47124000 Europe 70.510 10560.486
## 532 France 1967 49569000 Europe 71.550 12999.918
## 533 France 1972 51732000 Europe 72.380 16107.192
## 534 France 1977 53165019 Europe 73.830 18292.635
## 535 France 1982 54433565 Europe 74.890 20293.897
## 536 France 1987 55630100 Europe 76.340 22066.442
## 537 France 1992 57374179 Europe 77.460 24703.796
## 538 France 1997 58623428 Europe 78.640 25889.785
## 539 France 2002 59925035 Europe 79.590 28926.032
## 540 France 2007 61083916 Europe 80.657 30470.017</code></pre>
<p>Always always always plot the data. Yes, even now.</p>
<pre class="r"><code>p <- ggplot(j_dat, aes(x = year, y = lifeExp))
p + geom_point() + geom_smooth(method = "lm", se = FALSE)</code></pre>
<p><img src="block012_function-regress-lifeexp-on-year_files/figure-html/first-example-scatterplot.png" /></p>
</div>
<div id="get-some-code-that-works" class="section level3">
<h3>Get some code that works</h3>
<p>Fit the regression</p>
<pre class="r"><code>j_fit <- lm(lifeExp ~ year, j_dat)
coef(j_fit)
## (Intercept) year
## -397.7646019 0.2385014</code></pre>
<p>Whoa, check out that crazy intercept! Apparently the life expectancy in France around year 0 A.D. was minus 400 years! Never forget to sanity check a model. In this case, a reparametrization is in order. I think it makes more sense for the intercept to correspond to life expectancy in 1952, the earliest date in our dataset. Estimate the intercept eye-ball-o-metrically from the plot and confirm that we’ve got something sane and interpretable now.</p>
<pre class="r"><code>j_fit <- lm(lifeExp ~ I(year - 1952), j_dat)
coef(j_fit)
## (Intercept) I(year - 1952)
## 67.7901282 0.2385014</code></pre>
<div id="sidebar-regression-stuff" class="section level4">
<h4>Sidebar: regression stuff</h4>
<p>There are two things above that might prompt questions.</p>
<p>First, how did I know to get the estimated coefficients from a fitted model via <code>coef()</code>? Years of experience. But how might a novice learn such things? Read <a href="http://www.rdocumentation.org/packages/stats/functions/lm">the documentation for <code>lm()</code></a>, in this case. The “See also” section advises us about many functions that can operate on fitted linear model objects, including, but by no means limited to, <code>coef()</code>. Read <a href="http://www.rdocumentation.org/packages/stats/functions/coef">the documentation on <code>coef()</code></a> too.</p>
<p>Second, what am I doing here: <code>lm(lifeExp ~ I(year - 1952))</code>? I want the intercept to correspond to 1952 and an easy way to accomplish that is to create a new predictor on the fly: year minus 1952. The way I achieve that in the model formula, <code>I(year - 1952)</code>, uses the <code>I()</code> function which “inhibits interpretation/conversion of objects”. By protecting the expression <code>year - 1952</code>, I ensure it is interpreted in the obvious arithmetical way.</p>
</div>
</div>
<div id="turn-working-code-into-a-function" class="section level3">
<h3>Turn working code into a function</h3>
<p>Create the basic definition of a function and drop your working code inside. Add arguments and edit the inner code to match. Apply it to the practice data. Do you get the same result as before?</p>
<pre class="r"><code>le_lin_fit <- function(dat, offset = 1952) {
the_fit <- lm(lifeExp ~ I(year - offset), dat)
coef(the_fit)
}
le_lin_fit(j_dat)
## (Intercept) I(year - offset)
## 67.7901282 0.2385014</code></pre>
<p>I had to decide how to handle the offset. Given that I will scale this up to many countries, which, in theory, might have data for different dates, I chose to set a default of 1952. Strategies that compute the offset from data, either the main Gapminder dataset or the excerpt passed to this function, are also reasonable to consider.</p>
<p>I loathe the names on this return value. This is not my first rodeo and I know that, downstream, these will contaminate variable names and factor levels and show up in public places like plots and tables. Fix names early!</p>
<pre class="r"><code>le_lin_fit <- function(dat, offset = 1952) {
the_fit <- lm(lifeExp ~ I(year - offset), dat)
setNames(coef(the_fit), c("intercept", "slope"))
}
le_lin_fit(j_dat)
## intercept slope
## 67.7901282 0.2385014</code></pre>
<p>Much better!</p>
</div>
<div id="test-on-other-data-and-in-a-clean-workspace" class="section level3">
<h3>Test on other data and in a clean workspace</h3>
<p>It’s always good to rotate through examples during development. The most common error this will help you catch is when you accidentally hard-wire your example into your function. If you’re paying attention to your informal tests, you will find it creepy that your function returns <strong>exactly the same results</strong> regardless which input data you give it. This actually happened to me while I was writing this document, I kid you not! I had left <code>j_fit</code> inside the call to <code>coef()</code>, instead of switching it to <code>the_fit</code>. How did I catch that error? I saw the fitted line below, which clearly did not have an intercept in the late 60s and a positive slope, as my first example did. Figures are a mighty weapon in the fight against nonsense. I don’t trust analyses that have few/no figures.</p>
<pre class="r"><code>j_country <- "Zimbabwe"
(j_dat <- subset(gDat, country == j_country))
## country year pop continent lifeExp gdpPercap
## 1693 Zimbabwe 1952 3080907 Africa 48.451 406.8841
## 1694 Zimbabwe 1957 3646340 Africa 50.469 518.7643
## 1695 Zimbabwe 1962 4277736 Africa 52.358 527.2722
## 1696 Zimbabwe 1967 4995432 Africa 53.995 569.7951
## 1697 Zimbabwe 1972 5861135 Africa 55.635 799.3622
## 1698 Zimbabwe 1977 6642107 Africa 57.674 685.5877
## 1699 Zimbabwe 1982 7636524 Africa 60.363 788.8550
## 1700 Zimbabwe 1987 9216418 Africa 62.351 706.1573
## 1701 Zimbabwe 1992 10704340 Africa 60.377 693.4208
## 1702 Zimbabwe 1997 11404948 Africa 46.809 792.4500
## 1703 Zimbabwe 2002 11926563 Africa 39.989 672.0386
## 1704 Zimbabwe 2007 12311143 Africa 43.487 469.7093
p <- ggplot(j_dat, aes(x = year, y = lifeExp))
p + geom_point() + geom_smooth(method = "lm", se = FALSE)</code></pre>
<p><img src="block012_function-regress-lifeexp-on-year_files/figure-html/second-example-scatterplot.png" /></p>
<pre class="r"><code>le_lin_fit(j_dat)
## intercept slope
## 55.22124359 -0.09302098</code></pre>
<p>The linear fit is comically bad, but yes I believe the visual line and the regression results match up.</p>
<p>It’s also a good idea to clean out the workspace, rerun the minimum amount of code, and retest your function. This will help you catch another common mistake: accidentally relying on objects that were lying around in the workspace during development but that are not actually defined in your function nor passed as formal arguments.</p>
<pre class="r"><code>rm(list = ls())
gDat <- read.delim("gapminderDataFiveYear.txt")
le_lin_fit <- function(dat, offset = 1952) {
the_fit <- lm(lifeExp ~ I(year - offset), dat)
setNames(coef(the_fit), c("intercept", "slope"))
}
le_lin_fit(subset(gDat, country == "Zimbabwe"))
## intercept slope
## 55.22124359 -0.09302098</code></pre>
</div>
<div id="are-we-there-yet" class="section level3">
<h3>Are we there yet?</h3>
<p>Yes.</p>
<p>Given how I plan to use this function, I don’t feel the need to put it under formal unit tests or put in argument validity checks. Let’s move on to <a href="http://stat545-ubc.github.io/block013_plyr-ddply.html">the exciting part</a>, which is scaling this up to <strong>all</strong> countries.</p>
</div>
<div class="footer">
This work is licensed under the <a href="http://creativecommons.org/licenses/by-nc/3.0/">CC BY-NC 3.0 Creative Commons License</a>.
</div>
</div>
<script>
// add bootstrap table styles to pandoc tables
$(document).ready(function () {
$('tr.header').parent('thead').parent('table').addClass('table table-condensed');
});
</script>
<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
(function () {
var script = document.createElement("script");
script.type = "text/javascript";
script.src = "https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
document.getElementsByTagName("head")[0].appendChild(script);
})();
</script>
</body>
</html>