Where to locate the Bayesian analysis

There are at least three options for anchoring the Bayesian analysis. I could:

• use the most recent election result and anchor the state of the hidden markov model for that day to the election result;
• decide that one or more pollsters on average is/are absolutely correct (that is to say individually or collectively they have a house effect of zero), or at least has a known constant house effect; or
• locate the poll estimate such that we assume the polling industry collectively is unbiased (that is to say, the sum of the house effects is zero).

Recently, I have been using the last approach. But others would argue that the first approach is most accurate. Writing here, Simon Jackman noted that election anchored models did better in 2010 than  sum-to-zero models. It looks like a similar outcome in the 2013 election.

In the 2013 election, I ran a number of different models. One model simply used the polling data from Newspoll and Neilsen, the two pollsters I judged the most accurate in the past. In effect, this was the second approach I outlined above. It ended up being very close to the actual result in 2013. However, I wont have that option in 2016 with the churn we are seeing in the polling industry at the moment.

The TPP-HMM line in the next chart (the fatter, brown line) is the sum-to-zero result for all houses excluding Essential. In the eight months prior to the 2013 election, Essential is largely below this line. Since the 2013 election, its poll estimate has been above the line on average. For those who use Essential in their analysis, this repositioning, whether by chance or some change in methodology, will impact on their anchored aggregation, which assumes no change in practice or house effect on average.

The next two charts are the sum-to-zero TPP chart, followed by the election anchored TPP chart for the period since 1 January 2013. In these charts, I have only used the Morgan multi-mode polls. I have not used Morgan phone nor face-to-face polls. The difference between the two charts is half a percentage point, with the 2013 election anchored chart being more favorable to the Coalition.

In the next chart we take the medians from the previous charts. We can see that these analyses are very similar (if not the same), except they are located half a percentage point apart on the y axis.

Pollster preference flows

I was wondering whether the difference between my TPP aggregations from the TPP polling and the primary vote polling was an artifact of the preference flows that the pollsters were applying to the primary vote estimates they had derived.

This is a question for a simple multiple regression against the formula:

TPP_estimate = coalition_pv + α green_pv + β other_pv

In English, the Coalition's two-party preferred vote-share estimate comprises the Coalition primary vote, plus a proportion of the Greens' primary vote ( α ), and a proportion of the other parties' primary vote ( β ). In this equation, α and β are both values between 0 and 1 (on the continuum of no flow of preferences through to 100% flow). I decided to solve this regression using a simple Bayesian model, as follows.

model {

    ## -- preference flows
    for(poll in 1:NUMPOLLS) { # for each poll result - rows
        yhat[poll] <- pv_coalition[poll] +
            (alpha * pv_greens[poll]) +
            (beta * pv_other[poll])
        y[poll] ~ dnorm(yhat[poll], tau)
    }

    ## priors
    alpha ~ dunif(0.0, 1.0)
    beta ~ dunif(0.0, 1.0)
    sigma ~ dunif(0.001, 0.1)
    tau <- pow(sigma, -2)
}

I undertook the analysis for each polling house, using their polling data since the last Federal election, with the following results.

On both charts, I have marked with a vertical gray line the preference flow I use in my models (0.1697 for the Greens and 0.533 for other parties). I used Antony Green's earlier work to set my preference flows within my models.

The gray line falls within the 95% credibility interval for each of the polling houses. Therefore, I cannot argue that any of the polling houses are using different preference flows from the one I am using. If the pollsters are using different preference flows, this test did not demonstrate that.

Update

I have re-ran this analysis with uninformative, rather than weakly-informative priors. The charts and model have been updated. The width of the credibility intervals for Nielsen and Ipsos are unsurprising, as they have 7 and 6 observations respectively.

Four Bayesian model comparison

I now have four different Bayesian models for estimating the two-party preferred (TPP) vote share for the Coalition from multiple polling house sources.

• A dynamic linear model with input from TPP polling data
• A dynamic linear model with input from the primary voting intention polling data
• A Beta-walk model with input from TPP polling data; and
• A Dirichlet-walk model with input from the primary voting intention polling data

While each is a hidden Markov model, they use different input data (the TPP estimate or the primary vote estimates from the polling houses), and they use different distributions to propagate the Markov process in the hidden Markov model (the Normal, Dirichlet, or Beta distributions). The models based on the Normal distribution are dynamic linear models (DLM).

Because of size constraints, only the models which use the pollsters' TPP input data have a daily hidden Markov model. In these models, we estimate the hidden Coalition TPP vote share for each day in the period under analysis.

The models that use the primary vote data from the pollsters track a weekly hidden Markov model. In these models, we estimate the hidden Coalition TPP vote share for each week in the period under analysis. The weeks in these models begin on the first day of the earliest poll in the input.

All of the models are processed using the free JAGS software. The fastest model is the TPP dynamic linear model. The slowest model (by a country mile) is the primary vote share dynamic linear model. In the context of model speed, I am thinking of testing the Biips software, which claims to be faster than JAGS. But any move from JAGS will take some time.

Differences can arise between the models for a host of reasons. Comparing the models can be useful for better understanding the underlying trends in voting preferences (and whether there are issues that need to be factored into future analysis).

At the moment there is a substantial divergence between the TPP models and the primary vote share models. My suspicion is that this has come about because of the collapse of the Palmer United Party, and the way in which pollsters allocate preferences based on preference-flows at the last election and the stated voting intention of poll respondents, but I will need to do some work to test this hypothesis.

Charting the collapse of the Palmer United Party

There was a time when the Palmer United Party (PUP) looked like a force to be reckoned with. It commanded 6.5 per cent of the primary vote, a seat in the House of Representatives, and Senate seats for Queensland, Tasmania and Western Australia. But since then, the PUP share of the primary vote has fallen dramatically.

For those who are interested in these things, I aggregated the Palmer primary vote polls using a Beta-walk model. The Beta-walk model follows.

model {

    #### -- observational model
    for(poll in 1:NUMPOLLS) { # for each poll result - rows
        adjusted_poll[poll] <- walk[pollDay[poll]] + houseEffect[house[poll]]
        palmerVotes[poll] ~ dbin(adjusted_poll[poll], n[poll])
    }

    #### -- temporal model (a daily walk where this today is much like yesterday)
    tightness <- 50000 # tightness of fit parameter
    for(day in 2:PERIOD) { # rows
        binomial[day] <- walk[day-1] * tightness
        walk[day] ~ dbeta(binomial[day], tightness - binomial[day])
    }

    ## -- weakly informative priors for first day in the temporal model
    alpha ~ dunif(1, 1500)
    walk[1] ~ dbeta(alpha, 10000-alpha)

    #### -- sum-to-zero constraints on house effects
    for(house in 2:HOUSECOUNT) { # for each house ...
        houseEffect[house] ~ dnorm(0, pow(0.1, -2))
    }
    houseEffect[1] <- -sum(houseEffect[2:HOUSECOUNT])
}

I did not include Palmer United in my primary vote models because Newspoll does not publish a primary vote estimate for Palmer United. Newspoll includes Palmer United in its other parties count.

Dirichlet-walk hidden Markov model

Forgive me, but I am going to talk technical for a bit. I have been playing with a latent Dirichlet process hidden Markov model for aggregating opinion polls of primary voting intention. But before we get there, let's locate this approach.

The broad description for the method of poll aggregation I use on this site is the hidden Markov model (HMM). These models are analogous to the state-space models developed in the 1960s (beginning with the Kalman Filter). Hidden Markov models are one subset of Bayesian hierarchical models.

Using a HMM, I model the population voting intention (which cannot be observed directly - it is "hidden") as a series of states (either daily or weekly, depending on the model). Each state is directly dependent on the previous state and a probability distribution linking the states. Collectively, these links form a Markov chain or process. The model is informed by irregular and noisy data from the selected polling houses.

Solving the model necessitates integration over a series of complex multidimensional probability distributions. The definite integral is typically impossible to solve algebraically. But it can be solved using a numerical method based on Markov chains and random numbers known as Markov Chain Monte Carlo (MCMC) integration. I use a free software product called JAGS to solve the model.

The two poll aggregation models I had developed previously were dynamic linear models. In both models, the probability distribution linking the hidden states in the Markov model was the normal or Gaussian distribution: the bell-curve of high school statistics. Similarly, the probability distribution linking the hidden states with the polling observations was also the normal distribution. It is this dual independent use of the normal distribution that makes these models: dynamic linear models.

For the model based on the univariate Coalition two-party preferred poll results, the dynamic linear model works a treat. However, the multivariate primary vote model was far more difficult to construct. Quite some effort went into constraining the primary vote shares so that they always summed to one. It resulted in a large (and very slow) directed acyclic graph.

To address this complexity problem, I wondered whether a Dirichlet distribution could be used. Named for Peter Gustav Lejeune Dirichlet, the distribution is pronounced dirik-lay. The Dirichlet distribution has a number of useful features: First it is a multivariate distribution. Second, the output from the distribution always sums to one. It sounded ideal for modelling a time series of dynamically changing, continuous primary vote proportions (all on the unit interval).

In 2013, Emil Aas Stoltenberg wrote a paper on Bayesian Forecasting of Election Results in Multiparty Systems. In that paper he explored a Dirichlet-Multinomial process for the hidden Markov model. His model was coded in Python.

The Coalition TPP estimate from the this model over 6-months, and over the period since the previous election is not dissimilar to the output from the two dynamic linear models.

Turning to the model, rather than estimate the poll result, I use the multinomial distribution to estimate the number of people in each poll sample who expressed a preference for each of the parties. This is a very different approach to my other models. So that you can see it, I will include the R-code where I set up the input data for the model.

It should go without saying that the usual caveats apply. This is new code, and may include bugs. You will note that I am testing a few options in different places (see comments).

df <- read.csv(args[6], header=TRUE)
df <- df[order(df$Week), ]
NUMPOLLS <- nrow(df)
PERIOD <- max(df$Week)
HOUSECOUNT <- length(levels(df$House)) # a factor
HOUSENAMES <- levels(df$House)
PARTYNAMES <- c('Coalition', 'Labor', 'Greens', 'Other')
PARTIES <- length(PARTYNAMES)
primaryVotes <- df[PARTYNAMES] * df$Sample
primaryVotes <- sapply(primaryVotes, function(x) round(x,0))
day0 <- min(as.Date(df$Date)) - 1

## Assume Coalition gets preferences as follows:
## - 16.97% of the Green vote [was 16.97 in 2013 and 21.16 in 2010]
## - 53.3% of the Other vote [was 53.3 in 2013 and 58.26 in 2010]
## See: Antony Green - http://blogs.abc.net.au/antonygreen/2013/11/
##   preference-flows-at-the-2013-federal-election.html
preference_flows <- c(1.0, 0.0, 0.1697, 0.533)


data = list(PERIOD = PERIOD,
        HOUSECOUNT = HOUSECOUNT,
        NUMPOLLS = NUMPOLLS,
        PARTIES = PARTIES,
        primaryVotes = primaryVotes,
        pollWeek = df$Week,
        house = as.integer(df$House),
        # manage rounding issues with df$Sample ...
        n = rowSums(primaryVotes),
        preference_flows = preference_flows
    )
print(data)


# ----- JAGS model ...
library(rjags)
model <- "
model {

    #### -- observational model
    for(poll in 1:NUMPOLLS) { # for each poll result - rows
        adjusted_poll[poll, 1:PARTIES] <- walk[pollWeek[poll], 1:PARTIES] +
            houseEffect[house[poll], 1:PARTIES]
        primaryVotes[poll, 1:PARTIES] ~ dmulti(adjusted_poll[poll, 1:PARTIES], n[poll])
    }

    #### -- temporal model (a weekly walk where this week is much like last week)
    tightness <- 10000 # KLUDGE: tightness of fit parameter selected by trial and error
    for(week in 2:PERIOD) {
        # Note: use math not a distribution to generate the multinomial ...
        multinomial[week, 1:PARTIES] <- walk[week-1,  1:PARTIES] * tightness
        walk[week, 1:PARTIES] ~ ddirch(multinomial[week, 1:PARTIES])
    }

    ## -- weakly informative priors for first week in the temporal model
    for (party in 1:2) { # for each major party
        alpha[party] ~ dunif(250, 600) # majors between 25% and 60%
    }
    for (party in 3:PARTIES) { # for each minor party
        alpha[party] ~ dunif(10, 250) # minors between 1% and 25%
    }
    walk[1, 1:PARTIES] ~ ddirch(alpha[])

    ## -- estimate a Coalition TPP from the primary votes
    for(week in 1:PERIOD) {
        CoalitionTPP[week] <- sum(walk[week, 1:PARTIES] *
            preference_flows[1:PARTIES])
    }

    #### -- sum-to-zero constraints on house effects
    for (party in 2:PARTIES) { # for each party ...
        # house effects across houses sum to zero
        # NOTE: ALL MUST SUM TO ZERO
        houseEffect[1, party] <- -sum( houseEffect[2:HOUSECOUNT, party] )
    }
    for(house in 1:HOUSECOUNT) { # for each house ...
        # house effects across the parties sum to zero
        houseEffect[house, 1] <- -sum( houseEffect[house, 2:PARTIES] )
    }
    # but note, we do not apply a double constraint to houseEffect[1, 1]
    monitorHouseEffectOneSumParties <- sum(houseEffect[1, 1:PARTIES])
    monitorHouseEffectOneSumHouses <- sum(houseEffect[1:HOUSECOUNT, 1])

    ## -- vague normal priors for house effects - centred on zero
    for (party in 2:PARTIES) { # for each party (cols)
        for(house in 2:HOUSECOUNT) { #  (rows)
            houseEffect[house, party] ~ dnorm(0, pow(0.1, -2))
       }
    }
}
"

jags <- jags.model(textConnection(model),
        data = data,
        n.chains=4,
        n.adapt=n_adapt
    )

The input for the 6-month model was as follows:

$PERIOD
[1] 26

$HOUSECOUNT
[1] 5

$NUMPOLLS
[1] 35

$PARTIES
[1] 4

$primaryVotes
      Coalition Labor Greens Other
 [1,]       532   574    154   140
 [2,]       560   518    168   154
 [3,]       350   410    115   125
 [4,]       439   450    139   127
 [5,]       385   385     95   135
 [6,]       375   395    120   110
 [7,]      1465  1483    417   325
 [8,]       504   602    154   140
 [9,]       532   560    154   154
[10,]       504   602    154   140
[11,]       355   415    120   110
[12,]       412   483    141   141
[13,]      1345  1450    392   312
[14,]       375   405    100   120
[15,]       448   448    142   142
[16,]       588   504    168   140
[17,]       390   380    115   115
[18,]       441   453    139   128
[19,]       380   400    110   110
[20,]       471   425    126   126
[21,]       957   979    278   205
[22,]       405   360    125   110
[23,]       546   532    182   126
[24,]       471   413    126   138
[25,]       385   380    120   115
[26,]      1008   995    301   228
[27,]       400   375    115   110
[28,]       457   410    141   164
[29,]       690   656    185   151
[30,]       603   491    182   126
[31,]       415   355    125   105
[32,]       464   429    139   128
[33,]      1307  1218    385   273
[34,]       410   370    130    90
[35,]       479   433    152   105

$pollWeek
 [1]  1  1  2  2  6  8  8  9  9 10 10 10 10 12 12 13 14 14 16 16 17 18 19 19 20
[26] 21 22 22 24 24 24 24 24 26 26

$house
 [1] 1 2 3 4 3 3 5 1 2 1 3 4 5 3 4 2 3 4 3 4 5 3 2 4 3 5 3 4 1 2 3 4 5 3 4

$n
 [1] 1400 1400 1000 1155 1000 1000 3690 1400 1400 1400 1000 1177 3499 1000 1180
[16] 1400 1000 1161 1000 1148 2419 1000 1386 1148 1000 2532 1000 1172 1682 1402
[31] 1000 1160 3183 1000 1169

$preference_flows
[1] 1.0000 0.0000 0.1697 0.5330

Update

This page has been updated. Updates were the result of further analysis, which identified glitches with the original model. The tightness of fit parameter in the temporal model was not resolving. I have returned to specifying the tightness of fit parameter in the model.

Second, I have removed the dmulti() - the multinomial distribution step - in the temporal model and replaced it with simple arithmetic to calculate the multinomial. This makes the model a simpler Dirichlet-walk.

Updated Bayesian analysis of opinion polls

There have been two polls in the past fortnight.

• Morgan Poll: which had a two-party preferred (TPP) estimate of 48 for the Coalition, with preferences distributed by how electors voted at the 2013 Federal Election (-0.5 on the previous Morgan Poll); and
• Newspoll: which had a TPP estimate of 48 for the Coalition (+1 on the previous Newspoll).

Before going to the Bayesian model, I noticed in this tweet from @PhantomTrend, that I was an outlier among those who aggregate polls (my estimate being the most favourable to Labor). I looked at my methodology and concluded there were two possible sources of difference that might contribute to this outcome.

First, I was using the Morgan series that is based on preferences distributed by how poll respondents say they will vote. In the following charts, we explore this difference a little further. These charts show that the self-reported preference allocations have typically been more favourable to Labor (and correspondingly less favourable to the Coalition). I have changed my processing arrangement to use the Morgan TPP series, with preferences distributed by how electors voted at the 2013 Federal Election.

I commented on the Morgan poll following the Budget, that its house effect appears to be changing (with a reduced lean away from the Coalition for five polls now). In the moving average in the first chart above (which looks past individual poll noise), we can see the Morgan series aligning with the Bayesian estimate over recent months. I have been concerned for sometime that Morgan is over-distributed (too bouncy) for the reported sample sizes. For modelling purposes, I already reduce the sample size to 1000. More recently, with five polls now closer to the aggregation, I am coming to the view there may be a methodology change in the way Morgan processes its raw data collection. This would represent a significant dis-continuity, which perhaps should be factored into the model. I am reviewing whether I continue to use Morgan in the poll aggregation.

The second reason my aggregation will differ from others is that I do not use Essential polls in the model. My concern with the Essential poll series is that it appears under-dispersed for the reported sample size (it is not bouncy enough). Looking at the reported weekly sample sizes (typically over 1000) and the combined fortnightly sample sizes (typically around 1800), I have wondered whether the weekly samples include some people from the previous week (or earlier). Anyway, (as noted here) because I cannot explain the under-dispersion, I do not use the Essential poll in my aggregation. Please note: I am not suggesting there is anything wrong with the Essential poll series. All I am saying is that I do not understand it fully.

In the next chart, we can see that the Essential poll (the red line) has typically been more favourable to the Coalition in comparison my aggregation (the fat brown line) which ignores it. If I included the Essential poll, the sum-to-zero constraint would move my aggregation towards the Coalition.



In terms of the aggregation, and not withstanding the changes noted above, the Coalition continues to improve following its February troubles, but it has some way to go before it is in winnable territory.

If we focus on the most recent six months, and run the model in respect of that data alone, we can see the following.

Refactored Bayesian model

Today, I updated the technicals page for the 2016 election. This is the page where I explain the Bayesian model I use. The page includes additional commentary on the JAGS model, which explains how the various elements of the model interact.

In the process, I re-examined the model I had been using. I became concerned that the prior I had been using for the house effect was not sufficiently uninformative.  I had been using a normal prior, centred on zero, with a standard deviation of five percentage points. I have changed this to a uniform prior between -15 and +15 percentage points. I would expect house effects to be in the range of -2 to +2 percentage points, so the 30 point uniform range should be uninformative.

The change has had little impact on the analysis, so my fears were probably unfounded. Nonetheless, I have retained the uniform prior, as it is clearly less informative than the normal prior.

The other change I have made to the charts is that I now include an extra band in the Bayesian output to indicate the middle 99 per cent of samples. Previously, I had only indicated the median sample with a line, and the ranges for the middle 95, 80 and 50 percent of the samples with increasingly darker shading.

Let's look at the outcome. For the three month analysis, the median estimate of voting intention is unchanged. The first chart is the revised chart, the second chart is the earlier analysis (from here).

Turning now to house effects. The first chart is the updated analysis. The second chart is the earlier analysis. The only significant difference is the extra band in the first chart, that indicates the range for the middle 99 per cent of samples.

These changes in the second decimal place for the medians are very small, and should be ignored.

Coalition TPP by Polling House

The Henderson Moving Average (HMA) has some significant limitations. Technically, it should only be used when the data points come in equally spaced periods of time, and it has no mechanism for dealing with missing data.

In what is not quite kosher analysis, I have applied a HMA to the houses that poll fairly regularly. Essential usually produces a weekly estimate. Morgan and Newspoll typically produce a fortnightly result. And ReachTEL yields a monthly estimate.

To obtain a rough six month moving average across these houses, I have applied a 25-term HMA to Essential, a 13-term HMA to Morgan and Newspoll, and a 7-term HMA to ReachTEL. The results follow. Only Morgan saw the post 2014 Budget as the more significant slump for the Coalition. The other three pollsters had the late January come early February 2015 slump as the more significant. All agree the Coalition has been improving since the early 2015 slump, but an election winning position in the polls will require further improvement.

Data and code for election 2016

I have made most of my data and code base available on Google Drive. Please note, this is my live code base, which I play with quite a bit. So, there will be times when it is broken or in some stage of being edited.

What I have not made available is the Excel spreadsheets into which I initially place my data. These live in the (hidden) raw-data directory. However, the collated data for the Bayesian model lives in the intermediate directory, visible from the above link.

The program that collates and organises the spreadsheets into a single CSV input file for the Bayesian analysis is TPP-step1.py. There are two intermediate input files (at the moment): TPP-3-stage1.csv and TPP-all-stage1.csv. The first of these input files covers the most recent three months. The second of these input files is for all polls since the 2013 election.

The Bayesian model itself lives in the file TPP-step2.R.

The code for producing the plots lives in TPP-step3.py.

The files that begin with the letter 'z' are bash shell scripts.

The most recent set of charts live in the graphs directory. I don't keep historical charts.

There are a handful of helper programs that live in the bin directory.

There are a couple files that I am working on in respect of a primary votes model. This is still a long way from finished.

If you see an error in my code or data, please drop me a line (comments below or email address in right hand column), and let me know. I can only improve with your help.

Morgan poll and Bayesian aggregation

Morgan was the final poll out of the blocks in the post Budget tsunami of opinion polls. But before we look at the Morgan poll, I want to reflect a little on house effects (the systemic polling bias for each pollster). 

The Bayesian model I use includes the assumption that the individual polling houses do not change their methodology (and consequently their systemic house bias) throughout the period under analysis. The model also assumes this bias is constant. Each house, on average, leans to Labor or the Coalition by fixed number of percentage points.

I have two problems with these assumptions:

• The first assumption is probably not true. It would be more reasonable to assume that polling houses are continually reviewing and from time-to-time improving their statistical practice. Unfortunately, it is rare for polling houses to expose their methodology changes to the public. So I cannot readily introduce discontinuities into the Bayesian model when polling practices change.

• Second, even if the polling houses did not change their statistical practice, there is no guarantee that the systemic house bias is constant. It might vary, for example, depending on the vote share of the parties.

I am reflecting on these assumptions because the past four Morgan polls have been a little more favourable to the Coalition than the earlier polls where on average. This may reflect a change of polling practice at Morgan. And it may be nothing more than the random noise associated with opinion polling. I simply do not know. However, it is a trend worthy of monitoring further.

One way to reduce the potential erroneous impact of the above assumptions in the Bayesian model is to reduce the time period under analysis. A shorter window of analysis is less likely to include methodology changes from polling houses. When changes do happen, they will pass through the window of analysis quickly. The voting intention for the period is also more likely to be in a narrower range. With the voting intention in a narrower range, the non-constant rate biases will be better modeled by a constant. However, you do not get something for nothing. With fewer data points under analysis, the precision of the model is reduced.

To help you judge how things stand at the end of the Budget period, I have run the model using the polling data for the past three months, as well as the polling data since the last election. You should note the differences in the median and precisions of the relative house effects (also noting that the longer period includes ACNielsen, which ceased polling in the middle of 2014).

Of some comfort, both models yield a very similar end-point in terms of the Coalition's two-party preferred vote share (47.8 or 47.9 per cent).

Bayesian model over three months

Just a reminder in respect of the above charts. The shading indicates the proportion of samples in the model. The Markov Chain Monte Carlo model is run 100,000 times. In each iteration, for each node in the model, a sample is drawn. There are nodes in the model for each day under analysis and the house effect for each polling organisation. In the charts:

• the palest shaded area represents the middle 95 per cent of the samples;
• the next shaded area represents the middle 80 per cent of the samples;
• the darkest shaded area represents the middle 50 per cent of sample;
• the dark line and white triangle is the median (middle) sample; and
• the white line on the above charts is a 61-term Henderson moving average.

Bayesian model since the last election