── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.2.1 ✔ readr 2.2.0
✔ forcats 1.0.1 ✔ stringr 1.6.0
✔ ggplot2 4.0.3 ✔ tibble 3.3.1
✔ lubridate 1.9.5 ✔ tidyr 1.3.2
✔ purrr 1.2.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(cmdstanr)
This is cmdstanr version 0.9.0.9000
- CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
- CmdStan path: /home/ken/.cmdstan/cmdstan-2.39.0
- CmdStan version: 2.39.0
Introduction
I have some data where different observations may be more or less reliable (in a way that I think I can quantify), and so each observation comes with a weight that expresses how reliable it is. I want to incorporate the weights into the estimation. It seems that the way to do this is to use a “pseudo-likelihood”: when you work out the log-likelihood, you take the contribution from each observation and multiply it by its case weight.
How does this go in Stan? There are two steps.
The toyest of toy examples
To illustrate, let’s use the simplest example we can think of: estimating a population mean where the observations are normally distributed and the variance is known. I’m doing a simple example like this so that we can focus on what happens to the Stan code in a simple case and not get burdened by details.
I put the following code in n1.stan:
// estimate normal mean with known SD 1
data {
int<lower = 1> n;
array[n] real y;
}
parameters {
real mu;
}
model {
// prior
mu ~ normal(0, 10);
// likelihood
y ~ normal(mu, 1);
}
and compile it:
n1 <-cmdstan_model(stan_file ="n1.stan")
I don’t normally like diffuse priors, but I have one here so that we can easily compare the output with what we are expecting.
To run this, we need some data. In the spirit of this being a toy example, we’ll have some toy data: exactly two observations, 1 and 3:
stan_data <-list(n =2, y =c(1, 3))
We’d expect the posterior mean to be about 2, maybe a bit less because the prior mean is zero (but not too much less because the prior is diffuse). Here and below, please forgive the voluminous Stan output that I have not figured out how to suppress:
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -1.52 -1.23 0.72 0.30 -2.91 -1.02 1.00 1786 1968
mu 1.96 1.97 0.70 0.69 0.79 3.12 1.00 1585 1881
and that seems to check out.
Modification 1: an explicit loop
If you go back and look at the Stan code in n1.stan, you’ll see that the model part is quietly vectorized: y is a vector of length n, and the likelihood line of the model actually says that y is a vector of n iid normals, not just one normal. This makes coding in Stan a lot more fluent, because you don’t have to think about individual observations. But you can do that if you want to. It turns out that we do need to write the loop explicitly in modification 2, so we will do it here (even though it seems rather pointless at this stage):
// likelihood
for (i in 1:n) {
y[i] ~ normal(mu, 1);
}
I saved the edited Stan code in n2.stan.
This is the exact equivalent of the earlier vectorized code, and so it should give the same results on the same data:
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -1.50 -1.23 0.70 0.29 -2.94 -1.02 1.00 1959 1928
mu 2.02 2.02 0.69 0.68 0.88 3.13 1.00 1300 1624
and so it does.
If you use vectorization, Stan does some optimization behind the scenes, so the loop will run more slowly, but the results will be the same.
Modification 2: explicitly constructing the likelihood
The “tilde” notation in Stan is actually a notational convenience: when specifying models, it simplifies the thinking to say “this is distributed as that”. But what is actually happening in Stan is that it is constructing a log-posterior density by adding together the log prior and the log likelihood (or, more precisely, the log prior and a contribution to the log likelihood for each observation). It does this by using a special quantity target which is initialized to zero and which you then increment with the right things. Stan has special functions to calculate the log probability density for all the distributions it knows about. The one for the normal distribution is called normal_lpdf. So if we modify our model section to do this calculation explicitly, we get this (in n3.stan):
model {
// prior
target += normal_lpdf(mu | 0, 10);
// likelihood
for (i in 1:n) {
target += normal_lpdf(y[i] | mu, 1);
}
}
and once again, running this on the same data should give the same result:
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -6.57 -6.30 0.70 0.30 -7.95 -6.08 1.00 2084 2550
mu 2.00 2.00 0.70 0.69 0.84 3.13 1.00 1635 2312
and it looks as if we do.
My understanding is that this code using target += is exactly equivalent to the code with the loop, and so should run as quickly as that.
Incorporating case weights
So far, it looks as if all we have done is to make our code less clear. But there is method to this madness: we wanted to have access to each observation’s contribution to the likelihood, so that we could multiply it by the case weight before it went in. But the term normal_lpdf(y[i] | mu, 1) is precisely the contribution of the \(i\)-th observation to the likelihood, and multiplying it by the case weight before incrementing target will handle the case weight.
So now we need to add to our data a vector w of case weights, and incorporate w[i] into the calculation of the log-likelihood:
data {
int<lower = 1> n;
array[n] real y;
array[n] real w;
}
parameters {
real mu;
}
model {
// prior
target += normal_lpdf(mu | 0, 10);
// likelihood
for (i in 1:n) {
target += normal_lpdf(y[i] | mu, 1) * w[i];
}
}
n4 <-cmdstan_model(stan_file ="n4.stan")
Then we set up the data. Let’s give the observation 3 twice the weight of the observation 1, so the posterior mean should be close to \((1 + 2(3)) / (1+2)\) = 2.33:
stan_data_4 <-list(n =2, y =c(1, 3), w =c(1, 2))fit.4<- n4$sample(data = stan_data_4)
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -7.86 -7.58 0.72 0.33 -9.29 -7.34 1.00 1630 2182
mu 2.34 2.33 0.59 0.60 1.38 3.32 1.00 1501 1411
and it looks as if we have successfully incorporated our case weights into the estimation.
Where I am going with this
I have some sports data collected over time for a number of players, and I want to estimate the players’ strengths at a certain time point, using the data I have up to that time point. It seems reasonable to expect that a player’s performance will fluctuate: a result from two years ago says less about the player’s level now than does a result from yesterday. I would expect this to be especially true for individual sports, as younger players gain experience or older players lose the ability to compete at the highest level. A natural thing to try is exponentially-decaying weights: for example, a result from a year ago has half the weight of a result today (which implies that a result from two years ago has a quarter the weight of a result today). Players’ strengths would be expected to change fairly slowly, so a gradual decay in the weight attached to a result seems better than (for example) simply removing all games from more than a year ago.
Of course, it would be better still to model the time-dependence of player strengths, but that is another story, one to tell another time.