{"id":824,"date":"2015-08-28T07:46:20","date_gmt":"2015-08-28T07:46:20","guid":{"rendered":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/?p=824"},"modified":"2025-02-26T13:21:37","modified_gmt":"2025-02-26T13:21:37","slug":"google-flu-trends","status":"publish","type":"post","link":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2015\/08\/28\/google-flu-trends\/","title":{"rendered":"Google Flu Trends"},"content":{"rendered":"

I am always on the lookout for interesting medical data sets that I can use\u00a0for teaching and the other day, quite by accident,\u00a0I came across “Google flu trends”. This is an internet project operated by Google that estimates the numbers of cases of influenza each week in 25 different countries. The records go back to the beginning of 2003.<\/p>\n

The novel aspect of\u00a0this project\u00a0is that influenza is monitored by tracking online Google searches related to influenza. The number of searches has been calibrated against surveys\u00a0of the number of medically diagnosed cases.\u00a0You can find\u00a0more information on Wikipedia<\/a>.<\/p>\n

There\u00a0is a fascinating on-going debate in the media and on the internet about the ability of such search engine data to predict cases of influenza. What is most interesting is the\u00a0emotion that this\u00a0debate generates. People who\u00a0are attracted by\u00a0big data get enthusiastic about\u00a0the project and play down its limitations,\u00a0while people who equate big data\u00a0with big brother stress the weakness in the flu estimates and overlook the potential. It is as though they are not really talking about Google flu trends but merely using\u00a0the debate\u00a0as a vehicle for justifying their wider prejudices.<\/p>\n

Anyway, I will just accept the data and see how we might model\u00a0them to detect influenza epidemics. After all, even if the search based predictions are\u00a0wrongly calibrated they might still rise when the true number of cases of flu rises, which could be enough to enable us to detect the start of an outbreak.\u00a0As an example, here is\u00a0the data for Chile.<\/p>\n

\"chile\"<\/a><\/p>\n

The dashed vertical lines denote periods of 52 weeks. Remember that Chile is in the southern hemisphere so the winter influenza peak occurs around the middle of the year. Just looking at the plot one can see a small seasonal trend of the type that is most obvious in year 4 (around 200 weeks). Then on top of that there is usually an annual epidemic that causes a sharp spike of variable size and\u00a0variable duration. For the purposes of monitoring influenza, we would like to analyse these data in real time and to decide as soon as possible when an epidemic is underway.<\/p>\n

I will try to derive a\u00a0Bayesian analysis of these data and in the process of developing and fitting\u00a0such a\u00a0model I will need to create a new sampler and integrate it into the general\u00a0structure described in my book.<\/p>\n

In my model the weekly counts\u00a0will be thought of as being dependent on a latent variable that denotes whether or not an\u00a0epidemic is in progress.\u00a0We will define\u00a0a variable Z that takes the value 0\u00a0when there is no current epidemic, and 1, when there is a current\u00a0epidemic. Of course, Z is not observed but has to be inferred from the pattern of counts.\u00a0Our aim is\u00a0to detect as soon as possible when Z switches from 0 to 1.<\/p>\n

It appears from the plot that it would be reasonable to start by assuming that there is exactly one flu epidemic every year\u00a0and that\u00a0the severity varies from the very minor (e.g. year 4) to the very extreme (year 7).<\/p>\n

Each\u00a0epidemic is characterised by the week that it starts, S, and the week that it ends, E. These parameters will take integer values and of course S must be less that E. It would be nice to be able to update\u00a0estimates of S and E\u00a0using a Metropolis-Hastings algorithm but the programs given in my book only work for parameters defined over a continuous range: mhsnorm<\/strong> is good for parameters defined over an infinite range (-\u221e,\u221e), mhslogn<\/strong> is good for (0,\u221e) and mhstrnc<\/strong> is good for (a,b), where a and b are finite. Our first job therefore is to design a new Metropolis-Hastings sampler that can deal with integer valued parameters.<\/p>\n

Let\u2019s suppose that \u03b8 is a parameter defined over the integers (1, 2, \u2026 , n) and and that we have a program that works out the log-posterior given \u03b8. If the current value of the parameter is \u03b80<\/sub>, we need to make a proposal say \u03b81<\/sub> and compare logpost(\u03b80<\/sub>) with logpost(\u03b81<\/sub>) to see if we accept a move to \u03b81<\/sub> or reject it and stay at \u03b80<\/sub>.<\/p>\n

We might propose a random move to any other integer in the range [1,n]. The disadvantage of completely random moves\u00a0is that\u00a0many of the proposed moves will take us to poorly supported values and such\u00a0moves\u00a0are very likely to be rejected. Lots of rejected moves will mean poor mixing.<\/p>\n

An alternative strategy would be to propose a move to a neighbouring value, i.e. from \u03b80<\/sub>=12 to \u03b81<\/sub>=11 or 13. Such local moves are more likely to be accepted but\u00a0the change is small\u00a0so the algorithm will\u00a0require lots of moves to cover the posterior and this might also be inefficient, especially early in the chain when we have not located the peak of the posterior.<\/p>\n

On balance my gut feeling\u00a0favours short moves to neighbouring points. So let\u2019s assume that when \u03b80<\/sub>=12 we are equally likely to suggest a move up to 13 or a move down to 11. Similarly if \u03b80<\/sub>=11 we are equally likely to propose a move to 10 or 12. The key feature is Pr(propose 11 | \u03b80<\/sub>=12 ) = Pr(propose 12 | \u03b80<\/sub>=11 ). This is the characteristic\u00a0that defines\u00a0the special case known as a\u00a0Metropolis algorithm in which the proposal probabilities cancel and we only need to worry about the comparative values of the log-posterior.<\/p>\n

There is a small problem at the end points of the range. Consider the case\u00a0when \u03b80<\/sub>=1, the\u00a0move to\u00a0\u03b81<\/sub>=0 is not allowed because the range is [1,n]. We can cope with this, either by\u00a0modifying the rule for making proposals, or by making the log-posterior\u00a0at \u03b81<\/sub>=0 so small that the move is never accepted. The latter is easier to program.<\/p>\n

The program for the new sampler will follow the same pattern as mhsnorm<\/strong>. So let\u2019s look at that program and see what needs to be changed. I have highlighted the few places that\u00a0need our attention\u00a0in green.<\/p>\n

program mhsnorm<\/strong><\/span>, rclass<\/span>
\n\u00a0\u00a0 syntax anything [fweight] [if] [in],\u00a0 <\/span>SD(string)<\/span>\u00a0\u00a0\u00a0<\/strong><\/span>\u00a0\u00a0 [\u00a0<\/span>\u00a0\u00a0 LOGP(real 1e20)<\/span>\u00a0\u00a0 DEBUG\u00a0\u00a0\u00a0\u00a0<\/span>]<\/span>
\n\u00a0\u00a0 tokenize “`anything'”<\/span>
\n\u00a0\u00a0 local logpost “`1′”<\/span>
\n\u00a0\u00a0 local b\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 “`2′”<\/span>
\nlocal ipar\u00a0\u00a0 = `3′<\/span>
\n\u00a0\u00a0 tempname logl newlogl<\/span>
\n\u00a0\u00a0 if “`debug'” == “” local qui “qui”<\/span>
\n\u00a0\u00a0 local oldvalue = `b'[1,`ipar’]<\/span>
\n\u00a0\u00a0 scalar `logl’ = .<\/span>
\n\u00a0\u00a0 if `logp’ < 1e19 scalar `logl’ = `logp’<\/span>
\n\u00a0 if `logl’ == . {<\/span>
\n\u00a0\u00a0\u00a0\u00a0 `qui’ `logpost’ `logl’ `b’ `ipar’ [`weight’`exp’] `if’ `in’<\/span>
\n\u00a0\u00a0\u00a0\u00a0 if `logl’ == . {<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 di as err “mhsnorm<\/strong><\/span>: cannot calculate log-posterior for parameter `ipar'”<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 matrix list `b’<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 exit(1400)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 }<\/span>
\n\u00a0\u00a0 }<\/span>
\n\u00a0\u00a0 matrix `b'[1,`ipar’] = `oldvalue’ + `sd’*rnormal()<\/strong><\/span><\/span>
\n\u00a0\u00a0 scalar `newlogl’ = .<\/span>
\n\u00a0\u00a0 `qui’ `logpost’ `newlogl’ `b’ `ipar’ [`weight’`exp’] `if’ `in’<\/span>
\n\u00a0\u00a0 if `newlogl’ == . {<\/span>
\n\u00a0\u00a0\u00a0\u00a0 di as err “mhsnorm<\/strong><\/span>: cannot calculate log-posterior for parameter `ipar'”<\/span>
\n\u00a0\u00a0\u00a0\u00a0 matrix list `b’<\/span>
\n\u00a0\u00a0\u00a0\u00a0 exit(1400)<\/span>
\n\u00a0\u00a0 }<\/span>
\n\u00a0\u00a0 if log(uniform()) < (`newlogl’-`logl’) local accept = 1<\/span>
\n\u00a0\u00a0 else {<\/span>
\n\u00a0\u00a0\u00a0\u00a0 local accept = 0<\/span>
\n\u00a0\u00a0\u00a0\u00a0 matrix `b'[1,`ipar’] = `oldvalue’<\/span>
\n\u00a0\u00a0\u00a0\u00a0 scalar `newlogl’ = `logl’<\/span>
\n\u00a0\u00a0 }<\/span>
\n\u00a0\u00a0 return scalar accept = `accept’<\/span>
\n\u00a0\u00a0 return scalar logp = `newlogl’<\/span>
\nend<\/span><\/p>\n

Most of the program is involved in housekeeping and error checking so the changes are minimal. First, since I\u00a0am going to call the new sampler mhsint<\/strong>,\u00a0we need to change all references to mhsnorm into mhsint, then we need to drop the redundant SD option and finally we change the line that makes the new proposal from<\/p>\n

matrix `b'[1,`ipar’] = `oldvalue’ + `sd’*rnormal()<\/strong><\/span><\/p>\n

to<\/p>\n

matrix `b'[1,`ipar’] = `oldvalue’ + 1 \u2013 2*(runiform()<0.5)<\/strong><\/span><\/p>\n

This ensures that the proposal has an equal chance of being 1 more or 1 less that the\u00a0previous value.<\/p>\n

Let\u2019s design a simple test to\u00a0show how mhsint can be used with mcmcrun in just the same way as the other samplers mentioned in my book. I will suppose that a variable can take the values 1 or 2 or 3 with posterior probabilities 0.2, 0.3 and 0.5. Here is a program for calculating the log-posterior. If we move outside the range [1,3] then the log-posterior is set to a very low number so as to effectively block the move.<\/p>\n

program logpost<\/span>
\n\u00a0\u00a0 args lnf b<\/span><\/p>\n

\u00a0 local k = `b'[1,1]<\/span>
\n\u00a0\u00a0 if `k’ == 1 scalar `lnf’ = log(0.2)<\/span>
\n\u00a0\u00a0 else if `k’ == 2 scalar `lnf’ = log(0.3)<\/span>
\n\u00a0\u00a0 else if `k’ == 3 scalar `lnf’ = log(0.5)<\/span>
\n\u00a0\u00a0 else scalar `lnf’ = -999.9<\/span>
\n end<\/span><\/p>\n

I\u2019ll start the run with the variable (that I\u2019m calling k) set to 2. This I use mcmcrun to create a short burnin and then 500 samples.<\/p>\n

matrix b = (2)<\/span>
\n mcmcrun logpost b using temp.csv, replace\u00a0 \/\/\/<\/span>
\n\u00a0\u00a0\u00a0\u00a0 param(k) burn(50) update(500) samp( mhsint )<\/span>
\n insheet using temp.csv , clear<\/span>
\n tabulate k<\/span><\/p>\n

          k |\u00a0\u00a0\u00a0\u00a0 Freq.\u00a0\u00a0\u00a0\u00a0 Percent\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 Cum.\r\n------------+-----------------------------------\r\n          1 |\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 102\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 20.40\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 20.40\r\n          2 |\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 145\u00a0\u00a0\u00a0\u00a0   29.00\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 49.40\r\n          3 |\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 253\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 50.60\u00a0\u00a0\u00a0\u00a0  100.00\r\n------------+-----------------------------------\r\n      Total |\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 500\u00a0\u00a0\u00a0\u00a0 100.00<\/pre>\n
Clearly the chain has settled to the correct proportions.<\/p>\n
Next let\u2019s design a test that is closer to the flu problem. Suppose that we generate at sequence of 50 values, i.e. the flu counts for each week over most of a year. At the start of the year the counts\u00a0will be Poisson with a mean of 2 but at some random point between 10 and 40 weeks they will switch to a Poisson with a mean of 4.<\/p>\n
clear<\/span>
\nset seed 626169<\/span>
\nset obs 50<\/span>
\nlocal start = int(10+runiform()*31)<\/span>
\ngen t = _n<\/span>
\ngen mu = 2 + 2*(t >= `start’)<\/span>
\ngen y = rpoisson(mu)<\/span>
\nscatter y t, c(l) xline(`start’)<\/span><\/p>\n
\"random50\"<\/a><\/p>\n
In this simulation the switching point is at 31. What we need to do is to estimate the time of the switch, which I\u2019ll again call k, from the 50 counts.<\/p>\n
Here is a program that calculates the log-posterior. If we step outside the range [10,40] then the\u00a0move is blocked by setting a very small log-posterior, otherwise we calculate a standard Poisson log-likelihood with known means 2 or 4 and a switch at the current k, which is stored as the first element of our parameter vector, b.<\/p>\n
program logpost<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 args lnf b<\/span><\/p>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if `b'[1,1] < 10 | `b'[1,1] > 40 scalar `lnf’ = -9999.9<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 else {<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 tempvar mu lnL<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 gen `mu’ = 2<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 qui replace `mu’ = 4 if _n >= `b'[1,1]<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 gen `lnL’ = y*log(`mu’) – `mu’<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 qui su `lnL’<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 scalar `lnf’ = r(sum)<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 }<\/span>
\nend<\/span><\/p>\n
We can run the analysis much as before<\/p>\n
matrix b = (25)<\/span>
\nmcmcrun logpost b using temp.csv, replace\u00a0 \/\/\/<\/span>
\n\u00a0\u00a0\u00a0\u00a0 param(k) burn(50) update(500) samp( mhsint )<\/span>
\ninsheet using temp.csv , clear<\/span>
\nhistogram k , start(9.5) width(1) xlabel(10(5)40)<\/span><\/p>\n
The estimated log-posterior does indeed peak at 31 but given the nature of the data that we saw in the time series plot it is not surprising that the switch looks more likely to have come later in the sequence rather than earlier.<\/p>\n
\"posthistogram\"<\/a><\/p>\n
If we inspect the trace plot of the chain then we can see that we only made steps of size one but,\u00a0despite this,\u00a0the mixing looks reasonable.<\/p>\n
\"flutrace1\"<\/a><\/p>\n
It would be a simple matter to change from steps of size 1 to steps of size 1 or 2 or steps of size 1 or 2 or 3. We just need to change the line of code\u00a0that creates the random proposal.<\/p>\n
As it stands, the program logpost does not include a prior on k, so the implicit assumption is that prior to seeing the data we thought that all switching points in the range [10,40] were equally likely.\u00a0As it happens, this\u00a0accords with the way that the data were generated but had the data been real counts of flu in a small population, we might well have had prior belief about when the switch would occur. For example, we\u00a0might believe\u00a0that 25 is the most likely switching point and that the prior probability is proportional to 1\/(1+abs(k-25))<\/p>\n
clear<\/span>
\nrange k 10 40 31<\/span>
\ngen prior = 1\/(1+abs(k-25))<\/span>
\ntwoway bar prior k<\/span><\/p>\n
\"priorhistogram\"<\/p>\n
We can incorporate this prior into the calculation of the log-posterior<\/p>\n
program logpost<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 args lnf b<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if `b'[1,1] < 10 | `b'[1,1] > 40 scalar `lnf’ = -9999.9<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 else {<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 tempvar mu lnL<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 gen `mu’ = 2<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 qui replace `mu’ = 4 if _n >= `b'[1,1]<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 gen `lnL’ = y*log(`mu’) – `mu’<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 qui su `lnL’<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 scalar `lnf’ = r(sum) – log(1+abs(`b'[1,1]-25))<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 }<\/span>
\nend<\/span><\/p>\n
and the posterior is pulled slightly to the left.\"posthistogram2\"<\/p>\n
Next time I will use the new sampler to help me build a model for the flu data.<\/p>\n","protected":false},"excerpt":{"rendered":"
I am always on the lookout for interesting medical data sets that I can use\u00a0for teaching and the other day, quite by accident,\u00a0I came across “Google flu trends”. This is an internet project operated by Google that estimates the numbers of cases of influenza each week in 25 different countries. The records go back to […]<\/p>\n","protected":false},"author":134,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[80,40],"class_list":["post-824","post","type-post","status-publish","format-standard","hentry","category-uncategorized","tag-google-flu-trends","tag-metropolis-hastings"],"_links":{"self":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/824","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/users\/134"}],"replies":[{"embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/comments?post=824"}],"version-history":[{"count":13,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/824\/revisions"}],"predecessor-version":[{"id":853,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/824\/revisions\/853"}],"wp:attachment":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/media?parent=824"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/categories?post=824"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/tags?post=824"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}