<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/were-there-more-notable-deaths-than-expected-in-2016/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://masalmon.github.io/2017/02/12/wikideaths1_ts/"> Maëlle</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>After exploring my <a target="_blank" rel="nofollow" href="http://www.masalmon.eu/2017/02/12/wikideaths2_population/">study population</a> of Wikipedia deaths, I want to analyse the time series of monthly counts of notable deaths. This is not a random interest of mine, <a target="_blank" rel="nofollow" href="https://edoc.ub.uni-muenchen.de/19877/">my PhD thesis</a> was about monitoring time series of count, the application being weekly number of reported cases of various diseases.</p>
<p><span></span></p>
<h1>Number of deaths reported over time</h1>
<div>
<pre><code><span>library</span><span>(</span><span>&quot;ggplot2&quot;</span><span>)</span><span>
</span><span>library</span><span>(</span><span>&quot;viridis&quot;</span><span>)</span><span>
</span><span>library</span><span>(</span><span>&quot;dplyr&quot;</span><span>)</span><span>
</span><span>library</span><span>(</span><span>&quot;lubridate&quot;</span><span>)</span><span>
</span><span>deaths</span><span> </span><span>&lt;-</span><span> </span><span>readr</span><span>::</span><span>read_csv</span><span>(</span><span>&quot;data/deaths_with_demonyms.csv&quot;</span><span>)</span><span>
</span><span>deaths</span><span> </span><span>%>%</span><span>
</span><span>group_by</span><span>(</span><span>date</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>n_deaths</span><span> </span><span>=</span><span> </span><span>n</span><span>())</span><span> </span><span>%>%</span><span>
</span><span>ggplot</span><span>(</span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>date</span><span>,</span><span> </span><span>xend</span><span> </span><span>=</span><span> </span><span>date</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>0</span><span>,</span><span> </span><span>yend</span><span> </span><span>=</span><span> </span><span>n_deaths</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_segment</span><span>()</span><span> </span><span>+</span><span>
</span><span>xlab</span><span>(</span><span>&quot;Time (days)&quot;</span><span>)</span><span> </span><span>+</span><span> </span><span>ylab</span><span>(</span><span>&quot;Number of reported deaths&quot;</span><span>)</span><span>
</span></code></pre>
</div>
<p><img alt="plot of chunk unnamed-chunk-1" src="https://i1.wp.com/masalmon.github.io/figure/source/2017-02-12-wikideaths1_ts/unnamed-chunk-1-1.png?w=456" data-recalc-dims="1"></p>
<p>I have two remarks about this figure:</p>
<ul>
<li>
<p>There is an increase in the number of reported deaths over time, which is not surprising given the growth of Wikipedia.</p>
</li>
<li>
<p>There are three major peaks. They are due to the <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/2010_Haiti_earthquake">2010 Haiti earthquake</a> of the 12th of January; the <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/2010_Polish_Air_Force_Tu-154_crash">2010 Polish Air Force Tu-154 crash</a> on the 10th of April of 2010; and the <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/2011_Lokomotiv_Yaroslavl_air_disaster">2011 Lokomotiv Yaroslavl air disaster</a> on the 7th of September 2011. In the following, I will remove the fatalities from these events.</p>
</li>
</ul>
<h1>Modeling the time series</h1>
<p>When I think of modeling time series of count in R, I either think of <a target="_blank" rel="nofollow" href="https://www.jstatsoft.org/article/view/v070i10"><code>surveillance</code></a>, which was the implementational repository of my methodological developments <a target="_blank" rel="nofollow" href="http://staff.math.su.se/hoehle/blog/2016/12/12/hackinthedark.html">as Michael would say</a>; or of <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/tscount/vignettes/tsglm.pdf"><code>tscount</code></a>, which was the package of my PhD writing buddy Tobias Liboschik (after meeting at a conference we exchanged feedback on our thesis chapters which was awesome). If you have an interest in prospectively detecting outbreaks and a time series with seasonality, like the time series of number of cases of Salmonella, then <code>surveillance</code> is the best choice. Now in the case of Wikipedia deaths, I felt more inclined towards using <code>tscount</code>. I had never properly used it apart from a comparison section in my thesis so it was also more exciting. So <code>tscount</code> it was!</p>
<p>In <code>tscount</code>, you have <code>tsglm</code> models which take into account both the count nature of the data and its time series nature, i.e. you can treat auto-correlation. The <a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/tscount/vignettes/tsglm.pdf">vignette of <code>tscount</code></a> is a really great manuscript.</p>
<p>I’ll use the time series of monthly counts of deaths in the rest of the post.</p>
<div>
<pre><code><span>library</span><span>(</span><span>&quot;tscount&quot;</span><span>)</span><span>
</span><span>monthly_deaths</span><span> </span><span>&lt;-</span><span> </span><span>deaths</span><span> </span><span>%>%</span><span>
</span><span>group_by</span><span>(</span><span>wiki_link</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>mutate</span><span>(</span><span>month</span><span> </span><span>=</span><span> </span><span>update</span><span>(</span><span>date</span><span>,</span><span> </span><span>day</span><span> </span><span>=</span><span> </span><span>1</span><span>))</span><span> </span><span>%>%</span><span>
</span><span>ungroup</span><span>()</span><span> </span><span>%>%</span><span>
</span><span>group_by</span><span>(</span><span>month</span><span>)</span><span> </span><span>%>%</span><span>
</span><span>summarize</span><span>(</span><span>n</span><span> </span><span>=</span><span> </span><span>n</span><span>())</span><span> </span><span>haiti_earthquake</span><span> </span><span>&lt;-</span><span> </span><span>update</span><span>(</span><span>ymd</span><span>(</span><span>&quot;2010-01-12&quot;</span><span>),</span><span> </span><span>day</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span><span>aircrash1</span><span> </span><span>&lt;-</span><span> </span><span>update</span><span>(</span><span>ymd</span><span>(</span><span>&quot;2010-01-12&quot;</span><span>),</span><span> </span><span>day</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span><span>aircrash2</span><span> </span><span>&lt;-</span><span> </span><span>update</span><span>(</span><span>ymd</span><span>(</span><span>&quot;2011-09-07&quot;</span><span>),</span><span> </span><span>day</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span><span>monthly_deaths</span><span> </span><span>&lt;-</span><span> </span><span>mutate</span><span>(</span><span>monthly_deaths</span><span>,</span><span> </span><span>n</span><span> </span><span>=</span><span> </span><span>ifelse</span><span>(</span><span>month</span><span> </span><span>==</span><span> </span><span>haiti_earthquake</span><span>,</span><span> </span><span>n</span><span> </span><span>-</span><span> </span><span>10</span><span>,</span><span> </span><span>n</span><span>))</span><span>
</span><span>monthly_deaths</span><span> </span><span>&lt;-</span><span> </span><span>mutate</span><span>(</span><span>monthly_deaths</span><span>,</span><span> </span><span>n</span><span> </span><span>=</span><span> </span><span>ifelse</span><span>(</span><span>month</span><span> </span><span>==</span><span> </span><span>aircrash1</span><span>,</span><span> </span><span>n</span><span> </span><span>-</span><span> </span><span>96</span><span>,</span><span> </span><span>n</span><span>))</span><span>
</span><span>monthly_deaths</span><span> </span><span>&lt;-</span><span> </span><span>mutate</span><span>(</span><span>monthly_deaths</span><span>,</span><span> </span><span>n</span><span> </span><span>=</span><span> </span><span>ifelse</span><span>(</span><span>month</span><span> </span><span>==</span><span> </span><span>aircrash2</span><span>,</span><span> </span><span>n</span><span> </span><span>-</span><span> </span><span>44</span><span>,</span><span> </span><span>n</span><span>))</span><span>
</span></code></pre>
</div>
<div>
<pre><code><span>ggplot</span><span>(</span><span>monthly_deaths</span><span>)</span><span> </span><span>+</span><span>
</span><span>geom_segment</span><span>(</span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>month</span><span>,</span><span> </span><span>xend</span><span> </span><span>=</span><span> </span><span>month</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>0</span><span>,</span><span> </span><span>yend</span><span> </span><span>=</span><span> </span><span>n</span><span>))</span><span> </span><span>+</span><span>
</span><span>ylab</span><span>(</span><span>&quot;Number of reported deaths&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>xlab</span><span>(</span><span>&quot;Time (months)&quot;</span><span>)</span><span>
</span></code></pre>
</div>
<p><img alt="plot of chunk unnamed-chunk-3" src="https://i0.wp.com/masalmon.github.io/figure/source/2017-02-12-wikideaths1_ts/unnamed-chunk-3-1.png?w=456" data-recalc-dims="1"></p>
<p>We still see the upward trend, and I have a feeling we’re also dealing with overdispersion. First step before doing any modeling: getting a ts object. I hold off the last 12 months which are 2016 months.</p>
<div>
<pre><code><span>ts_deaths</span><span> </span><span>&lt;-</span><span> </span><span>xts</span><span>::</span><span>xts</span><span>(</span><span>monthly_deaths</span><span>$</span><span>n</span><span>[</span><span>1</span><span>:</span><span>144</span><span>],</span><span> </span><span>monthly_deaths</span><span>$</span><span>month</span><span>[</span><span>1</span><span>:</span><span>144</span><span>])</span><span>
</span><span>ts_deaths</span><span> </span><span>=</span><span> </span><span>ts</span><span>(</span><span>ts_deaths</span><span>,</span><span> </span><span>freq</span><span>=</span><span>12</span><span>,</span><span> </span><span>start</span><span>=</span><span>c</span><span>(</span><span>2004</span><span>,</span><span> </span><span>1</span><span>))</span><span>
</span><span>plot</span><span>(</span><span>ts_deaths</span><span>)</span><span>
</span></code></pre>
</div>
<p><img alt="plot of chunk unnamed-chunk-4" src="https://i2.wp.com/masalmon.github.io/figure/source/2017-02-12-wikideaths1_ts/unnamed-chunk-4-1.png?w=456" data-recalc-dims="1"></p>
<p>Then I use a model with regression on the past observation, on the mean from one year ago, and with a temporal trend. I log-transform time because the regression model uses a log link and I don’t expect the time trend to be exponential on the response scale. I’ll fit a Poisson and a negative binomial models in order to compare them.</p>
<div>
<pre><code><span>time</span><span> </span><span>&lt;-</span><span> </span><span>log</span><span>(</span><span>1</span><span>:</span><span>nrow</span><span>(</span><span>monthly_deaths</span><span>))</span><span>
</span><span>fit_pois</span><span> </span><span>&lt;-</span><span> </span><span>tsglm</span><span>(</span><span>ts_deaths</span><span>[</span><span>1</span><span>:</span><span>144</span><span>],</span><span> </span><span>model</span><span> </span><span>=</span><span> </span><span>list</span><span>(</span><span>past_obs</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>past_mean</span><span> </span><span>=</span><span> </span><span>13</span><span>),</span><span> </span><span>xreg</span><span> </span><span>=</span><span> </span><span>time</span><span>[</span><span>1</span><span>:</span><span>144</span><span>],</span><span> </span><span>distr</span><span> </span><span>=</span><span> </span><span>&quot;poisson&quot;</span><span>)</span><span>
</span><span>fit_nb</span><span> </span><span>&lt;-</span><span> </span><span>tsglm</span><span>(</span><span>ts_deaths</span><span>[</span><span>1</span><span>:</span><span>144</span><span>],</span><span> </span><span>model</span><span> </span><span>=</span><span> </span><span>list</span><span>(</span><span>past_obs</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>past_mean</span><span> </span><span>=</span><span> </span><span>13</span><span>),</span><span> </span><span>xreg</span><span> </span><span>=</span><span> </span><span>time</span><span>[</span><span>1</span><span>:</span><span>144</span><span>],</span><span> </span><span>distr</span><span> </span><span>=</span><span> </span><span>&quot;nbinom&quot;</span><span>)</span><span>
</span></code></pre>
</div>
<p>For choosing between the two models I’ll use scoring rules.</p>
<div>
<pre><code><span>rbind</span><span>(</span><span>Poisson</span><span> </span><span>=</span><span> </span><span>scoring</span><span>(</span><span>fit_pois</span><span>),</span><span> </span><span>NegBin</span><span> </span><span>=</span><span> </span><span>scoring</span><span>(</span><span>fit_nb</span><span>))</span><span>
</span></code></pre>
</div>
<div>
<pre><code>##         logarithmic    quadratic  spherical rankprob    dawseb     normsq
## Poisson    10.23428  0.003338134 -0.0495852 40.83497 17.633417 11.8105665
## NegBin      5.68638 -0.004157600 -0.0655695 35.61633  9.417219  0.9722222
##          sqerror
## Poisson 3835.296
## NegBin  3835.296
</code></pre>
</div>
<p>The smaller the scoring rule the better, so we should use the negative binomial model which I’m happy about since I had assumed we’d be dealing with overdispersion. Note that I used scoring rules whose values were computed on the part of the data I used for fitting the model, which is what Tobias Liboschik did in his manuscript for the Campylobacter example.</p>
<h1>Predicting 2016 values</h1>
<p>I realized there’s no <code>broom</code> “adapter” for <code>tsglm</code> so I’ll write untidy code and hope you’ll all forgive me.</p>
<div>
<pre><code><span>set.seed</span><span>(</span><span>1</span><span>)</span><span>
</span><span>pred2016</span><span> </span><span>&lt;-</span><span> </span><span>predict</span><span>(</span><span>fit_nb</span><span>,</span><span> </span><span>n.ahead</span><span> </span><span>=</span><span> </span><span>12</span><span>,</span><span> </span><span>level</span><span> </span><span>=</span><span> </span><span>0.9</span><span>,</span><span> </span><span>global</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span> </span><span>B</span><span> </span><span>=</span><span> </span><span>3000</span><span>,</span><span> </span><span>newxreg</span><span> </span><span>=</span><span> </span><span>time</span><span>[</span><span>145</span><span>:</span><span>156</span><span>])</span><span>
</span><span>monthly_deaths</span><span> </span><span>&lt;-</span><span> </span><span>mutate</span><span>(</span><span>monthly_deaths</span><span>,</span><span> </span><span>lower</span><span> </span><span>=</span><span> </span><span>NA</span><span>,</span><span>
</span><span>upper</span><span> </span><span>=</span><span> </span><span>NA</span><span>,</span><span> </span><span>pred</span><span> </span><span>=</span><span> </span><span>NA</span><span>)</span><span>
</span><span>monthly_deaths</span><span>$</span><span>lower</span><span>[</span><span>145</span><span>:</span><span>156</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>pred2016</span><span>$</span><span>interval</span><span>[,</span><span>1</span><span>]</span><span>
</span><span>monthly_deaths</span><span>$</span><span>upper</span><span>[</span><span>145</span><span>:</span><span>156</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>pred2016</span><span>$</span><span>interval</span><span>[,</span><span>2</span><span>]</span><span>
</span><span>monthly_deaths</span><span>$</span><span>pred</span><span>[</span><span>145</span><span>:</span><span>156</span><span>]</span><span> </span><span>&lt;-</span><span> </span><span>pred2016</span><span>$</span><span>pred</span><span>
</span></code></pre>
</div>
<p>Let’s plot the result.</p>
<div>
<pre><code><span>ggplot</span><span>(</span><span>monthly_deaths</span><span>)</span><span> </span><span>+</span><span>
</span><span>geom_segment</span><span>(</span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>month</span><span>,</span><span> </span><span>xend</span><span> </span><span>=</span><span> </span><span>month</span><span>,</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>0</span><span>,</span><span> </span><span>yend</span><span> </span><span>=</span><span> </span><span>n</span><span>))</span><span> </span><span>+</span><span>
</span><span>ylab</span><span>(</span><span>&quot;Number of reported deaths&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>xlab</span><span>(</span><span>&quot;Time (months)&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>(</span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>month</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>pred</span><span>),</span><span> </span><span>col</span><span> </span><span>=</span><span> </span><span>&quot;blue&quot;</span><span>)</span><span> </span><span>+</span><span>
</span><span>geom_ribbon</span><span>(</span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>month</span><span>,</span><span> </span><span>ymin</span><span> </span><span>=</span><span> </span><span>lower</span><span>,</span><span> </span><span>ymax</span><span> </span><span>=</span><span> </span><span>upper</span><span>),</span><span>
</span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.5</span><span>)</span><span>
</span></code></pre>
</div>
<p><img alt="plot of chunk unnamed-chunk-8" src="https://i1.wp.com/masalmon.github.io/figure/source/2017-02-12-wikideaths1_ts/unnamed-chunk-8-1.png?w=456" data-recalc-dims="1"></p>
<p>So what? All number of deaths of 2016 are in the 90% coverage prediction interval computed by <code>tscount</code> which means that maybe 2016 wasn’t a bad year for notable deaths <em>in general</em>, but then one should make the same analysis on a subset of the data like notable dead from the arts by filtering them according to their role, or notable dead whose Wikipedia pages consistently had more than X% of total Wikipedia daily page views before their deaths. I’m actually planning on doing more with page views, but in the meantime, it was nice to play with time series of counts again!</p>
<p>I’d like to end this post with a note from my husband, who thinks having a blog makes me an influencer. If you too like Wikipedia, consider <a target="_blank" rel="nofollow" href="https://wikimediafoundation.org/wiki/Ways_to_Give">donating to the Wikimedia foundation</a>.</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://masalmon.github.io/2017/02/12/wikideaths1_ts/"> Maëlle</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/E8hAmd_vXnE"></a></div>