In the comments to that post, I suggested that part of the problem may have been how I had defined the model.
Thinking about this some more, I think the problem was in this statement: << beta-for-pollster * (poll-result - minimum-poll-result) >>. Rather than the minimum, it should be a central tendency of some type (mid-point, mean, median, etc.).
I have now redefined the model to use the midpoint for each polling house (defined as the (min+max)/2 for that house). The use of house-specific mid-points acknowledges that the raw x scores include the alpha house effect I am trying to estimate. The revised model is:
house-effect-for-pollster = alpha-for-pollster + beta-for-pollster * (poll-result - pollster-mid-point)
The priors to the beta effect have been completely freed-up, so that they are uninformative. As a result, this model works better internally than the previous model.
However, I have added a constraint such that the aggregation assumes the "fizziness" factors also sum to zero. This may be a bollocks assumption; and will need some further analysis.
So - with the caveat that this is still very early exploratory analysis - the results of this new model follow.
On a plain reading, the above charts suggest that the current Rudd effect may actually be better for Labor than other aggregations have found. However, I am not convinced this is anything more than an artifact of the model. I am particularly concerned about the inclusion of data from polling houses that do not have data points across the highs and lows of the Gillard period, and across the Rudd and Gillard periods. Such data don't fully occupy the model, which can result in artifacts.
If we limit the analysis to data from Nielsen and Newspoll for a comparison using data that does fully occupy the model, we get the following.
Because I am committed to presenting poll analysis in a transparent and unbiased fashion, the next two code snippets show how I package up the data for the model (in R), and the model itself (in JAGS).
# - prepare the data ...
data$y <- data[ , y] / 100 # vote share between 0-1
data$variance <- data$y * (1 - data$y) / data$Sample #pq/n
data$standardError <- sqrt(data$variance)
data$samplePrecision <- 1 / data$variance
HOUSES <- levels(data$House)
HOUSECOUNT <- length(levels(data$House))
NUMPOLLS <- length(data$y)
cat('Number of polls: '); print(NUMPOLLS)
midpoints <- rep(NA, HOUSECOUNT)
for(i in seq_len(HOUSECOUNT)) {
midpoints[i] <- (max(data[data$House==HOUSES[i], 'y']) +
min(data[data$House==HOUSES[i], 'y'])) / 2
cat('Midpoints: '); cat(i); cat(' '); cat(HOUSES[i]); cat(': '); print(midpoints[i])
}
# - remember these dots
dotsFile <- paste('./files/', fPrefix, 'original.csv', sep='')
write.csv(data, file=dotsFile)
# - manage dates
day0 <- min(data$Date) - 1 # walk starts from earliest date
if(endDate == TODAY)
endDate <- max(data$Date)
PERIOD <- as.numeric(endDate - day0) # length of walk in days
DISCOUNTINUITYDAY <- as.numeric(as.Date(discontinuity) - day0)
cat('Discontinuity: '); print(DISCOUNTINUITYDAY)
tPrefix <- paste(format(day0+1, '%d-%b-%Y'), ' to ', format(endDate, '%d-%b-%Y'), sep='')
data$day <- as.numeric(data$Date - day0)
# - do the MCMC thing ...
parameters <- list(PERIOD = PERIOD,
HOUSECOUNT = HOUSECOUNT,
NUMPOLLS = NUMPOLLS,
DISCOUNTINUITYDAY = DISCOUNTINUITYDAY,
NEWSPOLL = which(levels(data$House) == 'Newspoll'),
y = data$y,
x = data$y,
day = data$day, house = as.integer(data$House),
samplePrecision = data$samplePrecision,
midpoints=midpoints
)
jags <- jags.model(textConnection(model),
data=parameters,
n.chains=4,
n.adapt=n.adapt
)
# - burn in
update(jags, n.iter=n.update) # burn-in the chains
# - capture results
jags.capture <- c('walk', 'alpha', 'beta', 'discontinuityValue')
coda.samples <- coda.samples(jags, jags.capture, n.iter=n.iter, thin=n.thin)
coda.matrix <- as.matrix(coda.samples)model {
## Derived from Simon Jackman's original model
## -- observational model
for(poll in 1:NUMPOLLS) {
# note: x and y are the original polling series
houseEffect[poll] <- alpha[house[poll]] +
beta[house[poll]]*(x[poll]-midpoints[house[poll]])
mu[poll] <- walk[day[poll]] + houseEffect[poll]
y[poll] ~ dnorm(mu[poll], samplePrecision[poll])
}
## -- temporal model
for(i in 2:PERIOD) { # for each day under analysis ...
day2DayAdj[i] <- ifelse(i==DISCOUNTINUITYDAY,
walk[i-1]+discontinuityValue, walk[i-1])
walk[i] ~ dnorm(day2DayAdj[i], walkPrecision)
}
sigmaWalk ~ dunif(0, 0.01) ## uniform prior on std. dev.
walkPrecision <- pow(sigmaWalk, -2) ## for the day-to-day random walk
walk[1] ~ dunif(0.4, 0.6) ## uninformative prior
discontinuityValue ~ dunif(-0.2, 0.2) ## uninformative prior
## -- house effects model
for(i in 2:HOUSECOUNT) { ## vague normal priors for house effects
alpha[i] ~ dunif(-0.1,0.1)
beta[i] ~ dunif(-1,1)
}
alpha[NEWSPOLL] <- -sum(alpha[2:HOUSECOUNT]) ## sum to zero
beta[NEWSPOLL] <- -sum(beta[2:HOUSECOUNT]) ## sum to zero
}
