Dear all,
this is my first post here, so if I have a too complex or too little explained minimal example model or made any other mistake, please tell me and I will retry.
Essentially, what I am trying to do is apply Linzer2013s reverse-random walk election forecasting model (as implemented in STAN by Kremp2016 here(see Andrew Gelman’s explanation of the model here) to the German electoral context (at this stage only one national level, but multinomial (multi-party) dependent variable, with inspiration credits to Stoetzer et al. 2017.
The key set-up of the model is as follows: a reverse random-walk starting with a structural/random prior on election day is estimated, with the log party support vector of today being a function of tomorrows value and some random variance. Polls on a given day are treated as multinomial draws from the latent (unobserved) party support vector at that moment.
Attached, I have provided a full minimal example with data-csv, R and STAN-script. Instead of describing every step at length and making this post very long, I have opted for providing the code along with detailed comments at every step. Please tell me if you had prefered a more detailed description at this place.
After a week of trial and error, the model finally builds, but does fails to compile with the following error message:
SAMPLING FOR MODEL ‘minimal_model’ NOW (CHAIN 2).
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Random variable is nan, but must not be nan! (in ‘modele4b066842876_minimal_model’ at line 78)
polls.csv (6.4 KB)
minimal_model.stan (3.2 KB)
minimal_example.R (1.8 KB)
Thank you so much for your help and warm welcome to your community.
Marcel
.
STAN MODEL
data {
int<lower=0> nParties; // number of parties
int<lower=0> nPolls; // number of polls
int<lower=0> nPeriods; // number of
vector[nPolls] samplesize; //number of Observations in a given sample
int date[nPolls]; //days until election
matrix[nParties, 2] normal_priors; // normal-prior matrix (alternatively, had tried it with beta priors, but did not work either)
int y[nPolls, nParties]; // outcome matrix: absolute counts of votes per party (column) per poll (row)
}
transformed data{
}
parameters {
real alphanorm[nParties]; //vector of structural priors
simplex[nParties] theta[nPolls];// theta to be used to as input for multinomial, obtained with dirichlet of vote shares
real s[nParties-1]; //evolution variance
}
transformed parameters {
matrix[nPeriods,nParties] alphastar; //daily party support vector: determined by alphanorm for election day, otherwise value for tomorrow plus random evolution variance
matrix[nPeriods,nParties] ea; // see below
matrix[nPeriods,nParties] alpha; // see below
real t[nParties-1]; //transformed evolution variance
real alrforcast[nParties]; //final forecast
real ef[nParties];
real forcast[nParties];
vector[nParties] shares;
//log transformation of structural priors
for (j in 1:(nParties)) {
alphastar[nPeriods,j] = log(alphanorm[j]/alphanorm[nParties]);
}
//alphastar for [,nParties] (other parties) set to 0
for (i in 1:(nPeriods-1)) {
alphastar[i,nParties] = 0;
}
//convert log shares from alphastar back to vote shares
for (i in 1:(nPeriods)) {
for (j in 1:(nParties)) {
ea[i,j] = exp(alphastar[i,j]);
alpha[i,j] = ea[i,j]/sum(ea[i,]);
}
}
//transposition from row vector to vector
for (k in 1:nPolls){
shares = alpha[date[k],]';
}
//transformation of parameter for evolution variance s
for (j in 1:(nParties-1)) {
t[j] = pow(s[j], -2);
}
//final forecast
for (j in 1:(nParties)) {
alrforcast[j] = alphastar[nPeriods,j];
ef[j] = exp(alrforcast[j]);
forcast[j] = ef[j]/sum(ef[]);
}
}
model {
// Walk starts on election day with forcast based on structural prior
for (j in 1:(nParties)) {
target += normal_lpdf(alphanorm[j] | normal_priors[j,1],normal_priors[j,1]);
}
// Backward-Random-Walk
for (i in 1:(nPeriods-1)) {
for (j in 1:(nParties-1)) {
//alphastar realization of today is result of tomorrow plus random variance t
target += normal_lpdf(alphastar[i,j] | alphastar[i+1,j],t[j]);
//this is where the error seems to appear:
//("Rejecting initial value: Error evaluating the log probability at the initial value.Exception: normal_lpdf:
// Random variable is nan, but must not be nan!")
}
}
// Likelihood Model
for (k in 1:nPolls) {
theta[k] ~ dirichlet(shares); //shares is the party support vote share vector for the day in which poll k was conducted
y[k,1:nParties] ~ multinomial(theta[k]); //y[k,] is the outcome matrix with absolute counts for poll responses
}
// evolution variance with uniform prior (see transformed parameters: t[j] = pow(s[j], -2);)
for (j in 1:(nParties-1)) {
s[j] ~ uniform(0,2);
}
}
R-Script
load("stanData.rdata")
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
polls <- read.csv("polls.csv")
party_names <- grep("cdu|spd|fdp|lin|afd|gru|other|oth", names(polls),value = T)
Y <- round(as.matrix(polls[,grep("cdu|spd|fdp|lin|afd|gru|other|oth", names(polls))]/100) * polls$sample_size)
NObs <- apply(Y,1,sum)
nParties <-ncol(Y)
time_window <- c(365)
# Normal Priors (from external structural forecast)
normal_priors <-
data.frame(
prior_mean = c(0.09096958,
0.35017978,
0.06393053,
0.07916824,
0.08401857,
0.05717943,
0.27455386),
prior_sigma = c(0.02276087,
0.04404920,
0.02241541,
0.02136944,
0.02218496,
0.02231709,
0.02953791)
)
rownames(normal_priors) <-
c("afd","cdu","fdp","gru","lin","oth","spd")
stanData <- list(y = Y,
nParties = nParties,
nPeriods = time_window[1]+1,
nPolls = nrow(Y),
date = polls$t,
samplesize = NObs,
normal_priors = normal_priors[party_names,] # Order others last
)
fit1 <- stan(
file = "minimal_model.stan", # Stan program
data = stanData, # named list of data
chains = 4, # number of Markov chains
warmup = 1000, # number of warmup iterations per chain
iter = 2000, # total number of iterations per chain
cores = 4, # number of cores (using 2 just for the vignette)
refresh = 100 # show progress every 'refresh' iterations
)
summary(fit1)