Following up on our post from the other day, we were trying to come up with parameters for the following measurement-error model:
– People have latent a distribution of attractiveness:
latent_attractiveness ~ normal(0, sigma_latent)
– Attractiveness for each person is measured twice, with different mean levels and different levels of added variance:
At time t, observed_attractiveness ~ normal(mu_obs_t + latent_attractiveness, sigma_obs_t), for t = 1, 2
– Attractiveness at each time is then recorded as a discrete variable, rounded to the nearest integer and constrained to fall between 1 and 5: That is, if observed attractiveness is less than 1.5, it is recorded as 1, if observed attractiveness is between 1.5 and 2.5, it is recorded as 2, etc.
– The only information easily accessible is recorded attractiveness binned into three categories: 1 or 2, 3 or 4, and 5. This summary is available for each person at both time points.
Here are the data, from the first two waves of the Add Health survey:
So, here’s the deal. We want to choose five parameters (sigma_latent, mu_obs_1, mu_obs_2, sigma_obs_1, sigma_obs_2) to fit 9 data points. Actually 8 data points because what we’re fitting are the proportions that sum to 1.
This shouldn’t be so hard, right?
Actually, it is kind of hard! At least, it’s hard to get a fit that gets real close to the data table.
1. Trial and error to find parameter values that match the data
I started with trial and error, first computing the mean and sd of the attractiveness ratings for each wave of the data (approximating the three binned levels as 1.5, 3.5, and 5 on the 1-to-5 scale), then setting mu_obs_1 and mu_obs_2 to the mean values and partitioning the variance so that sd_latent and sd_obs are roughly equal. In R:
discretize <- function(a) {
ifelse(a < 2.5, 1, ifelse(a < 4.5, 2, 3))
}
N_table <- rbind(c(208, 640, 40), c(501, 9879, 464), c(51, 1439, 1432))
N_data <- sum(N_table)
long_data <- array(NA, c(N_data, 2))
long_data[,1] <- rep(c(1,2,3, 1,2,3, 1,2,3), as.vector(N_table))
long_data[,2] <- rep(c(1,1,1, 2,2,2, 3,3,3), as.vector(N_table))
long_data_1to5 <- ifelse(long_data==1, 1.5, ifelse(long_data==2, 3.5, ifelse(long_data==3, 5, NA)))
means <- apply(long_data_1to5, 2, mean)
sds <- apply(long_data_1to5, 2, sd)
mu_obs <- means
sigma_latent <- sqrt(0.5*mean(sds^2))
sigma_obs <- sqrt(sds^2 - sigma_latent^2)
The result:
> print(round(mu_obs, 2))
[1] 3.68 3.59
> print(round(sigma_latent, 2))
[1] 0.54
> print(round(sigma_obs, 2))
[1] 0.61 0.45
I then estimate the probabilities in the 9 cells by simulating a million pairs of observations from this model, first simulating latent attractiveness, then the two observed attractiveness measures, then rounding:
N <- 1e6
y_latent <- rnorm(N, 0, sigma_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
y_obs[,j] <- rnorm(N, mu_obs[j] + y_latent, sigma_obs[j])
}
y_discrete <- discretize(y_obs)
And here's the result:
> print(table(y_discrete[,1], y_discrete[,2]))
1 2 3
1 18715 54792 393
2 41167 674479 55053
3 653 110776 43972
OK, it's hard to compare this to the data table because the numbers are so big, so let's rescale so that it sums to the number of people in the data:
> print(table(y_discrete[,1], y_discrete[,2])*N_data/N)
1 2 3
1 274.249610 802.921968 5.759022
2 603.261218 9883.815266 806.746662
3 9.569062 1623.311504 644.365688
This is hard to read so I'll round it:
> print(round(table(y_discrete[,1], y_discrete[,2])*N_data/N))
1 2 3
1 274 803 6
2 603 9884 807
3 10 1623 644
And, for convenience, let's look again at the data:
> print(N_table)
[,1] [,2] [,3]
[1,] 208 640 40
[2,] 501 9879 464
[3,] 51 1439 1432
Comparing the cells [1,3] and [3,1] in the simulated table to the observed data, we see that our simulations don't have enough variation between the two measurements. So let's rejigger by decreasing sigma_latent and increasing sigma_obs:
sigma_latent <- sqrt(0.25*mean(sds^2))
sigma_obs <- sqrt(sds^2 - sigma_latent^2)
Simulate again with these new parameter values:
y_latent <- rnorm(N, 0, sigma_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
y_obs[,j] <- rnorm(N, mu_obs[j] + y_latent, sigma_obs[j])
}
y_discrete <- discretize(y_obs)
And here's what we get:
> print(round(table(y_discrete[,1], y_discrete[,2])*N_data/N))
1 2 3
1 145 896 37
2 687 9592 1011
3 54 1827 404
Now the [1,3] and [3,1] entries are basically OK, but [1,1] and [3,3] are too small.
I tried some more tinkering with the five parameters but couldn't get the table to line up. And then I thought: Hey, why not cut out the middleman and just solve the optimization problem directly?
2. Stan program: first try
It goes like this:
data: the observed table, N_table
parameters: the five parameters, sigma_latent, mu_obs, sigma_obs
model: compute table(y_discrete) and compare to N_table
In this model, table(y_discrete)*N_data/N represents E(N_table|parameters), and we were computing this using simulation . . . but you can't simulate random data in the Stan model block! Or we could compute the expectations directly, but these are two-dimensional normal integrals, and there's no formula for that.
So I did a hack. You can't do simulations in the Stan model block, but you can simulate random numbers offline and pass them into the Stan model. Here's our Stan program, table_match.stan:
functions {
int discretize (real y) {
if (y < 2.5)
return 1;
else if (y > 2.5 && y < 4.5)
return 2;
else
return 3;
}
}
data {
array[3, 3] int N_table;
int N;
}
transformed data {
int N_data = 0;
for (j1 in 1:3) {
for (j2 in 1:3) {
N_data += N_table[j1,j2];
}
}
vector[N] z_latent;
matrix[N, 2] z_obs;
for (n in 1:N){
z_latent[n] = std_normal_rng();
for (j in 1:2) {
z_obs[n, j] = std_normal_rng();
}
}
real a = 1.*N_data/N;
}
parameters {
vector[2] mu_obs;
real s_latent;
vector[2] s_obs;
}
transformed parameters {
matrix[3,3] N_sim = rep_matrix(0, 3, 3);
for (n in 1:N){
real y_latent = s_latent * z_latent[n];
array[2] int y_discrete;
for (j in 1:2){
real y_obs = y_latent + mu_obs[j] + s_obs[j] * z_obs[n,j];
y_discrete[j] = discretize(y_obs);
}
N_sim[y_discrete[1], y_discrete[2]] += 1;
}
N_sim = a*N_sim;
}
model {
mu_obs ~ normal(3,3);
s_latent ~ exponential(1);
s_obs ~ exponential(1);
for (j1 in 1:3) {
for (j2 in 1:3) {
N_table[j1, j2] ~ poisson(N_sim[j1, j2]);
}
}
}
At the very end we're using the Poisson likelihood to compare the observed data to the model. And we include weak priors on the model parameters to keep everything stable.
We run the program from R:
library("cmdstanr")
table_match <- cmdstan_model("table_match.stan")
table_data <- list(N_table=N_table, N=1000)
fit <- table_match$sample(data=table_data, parallel_chains=4, refresh=0)
print(fit, max_rows=30)
Aaaand . . . it doesn't work:
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ 77677.57 77880.76 10570.22 13247.66 63077.85 91926.34 3.66 4 5
mu_obs[1] 0.24 0.10 0.61 0.51 -0.42 1.22 3.94 4 11
mu_obs[2] 0.47 0.40 1.64 2.24 -1.47 2.54 3.70 4 11
s_latent 2.78 1.82 2.52 1.37 0.49 7.00 3.42 4 12
s_obs[1] 1.84 1.60 1.18 1.08 0.46 3.68 3.47 4 13
s_obs[2] 1.65 1.70 1.36 1.93 0.18 3.03 3.43 4 11
N_sim[1,1] 8290.04 8865.67 1305.70 640.92 6096.06 9319.94 5.54 4 NA
N_sim[2,1] 890.16 593.49 669.06 358.48 351.70 2007.60 32.33 4 NA
N_sim[3,1] 468.94 80.60 711.48 65.18 14.65 1699.86 13.53 4 4
N_sim[1,2] 1759.81 1472.73 1368.41 1292.70 219.81 3883.31 5.95 4 4
N_sim[2,2] 852.42 622.80 655.06 369.34 219.81 1948.98 9.66 4 NA
N_sim[3,2] 322.47 307.73 251.19 336.75 29.31 644.78 8.22 4 4
N_sim[1,3] 626.65 256.44 779.88 260.71 14.65 1963.64 18.31 4 NA
N_sim[2,3] 408.84 300.41 339.12 271.58 58.62 952.51 16.39 4 NA
N_sim[3,3] 1034.67 381.00 1331.88 304.16 29.31 3326.46 16.98 4 NA
What about Stan in optimize setting?
fit_optimize <- table_match$optimize(data=table_data)
Nope:
Initial log joint probability = 69172.8
Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
1 69172.8 0 9.10807 0.001749 0.001 45
Chain 1 Optimization terminated with error:
Chain 1 Line search failed to achieve a sufficient decrease, no more progress can be made
Warning: Fitting finished unexpectedly! Use the $output() method for more information.
Damn.
3. Stan program: revision
What went wrong? The above code is mathematically correct, but the discreteness of all these simulations leads to difficulty with the Hamiltonian Monte Carlo. The problem is that, due to the use of fixed error terms and a finite N, continuous changes in the parameters result in jumps in the expected data and thus jumps in the likelihood.
But then I realized something: I can't compute two-dimensional normal integrals in closed form, but one-dimensional normal integrals are trivial. So, instead of passing in the vectors of random numbers z_latent, z_obs[,1], and z_obs[,2], I can just pass in z_latent, and then we can use the normal cdf function in Stan to integrate over z_obs[,1] and z_obs[,2]. This offers two advantages: first, by integrating over two of the three dimensions analytically rather than by simulation, we're computing the cell probabilities more accurately and we don't need such a large N; second, the target function will now be a smoothly continuous function of the parameters, so HMC should work better.
Here's our revised Stan program, table_match_2.stan:
data {
array[3, 3] int N_table;
int N;
}
transformed data {
int N_data = 0;
for (j1 in 1:3) {
for (j2 in 1:3) {
N_data += N_table[j1,j2];
}
}
vector[N] z_latent;
for (n in 1:N){
z_latent[n] = std_normal_rng();
}
real a = 1.*N_data/N;
}
parameters {
vector[2] mu_obs;
real s_latent;
vector[2] s_obs;
real sigma;
}
transformed parameters {
matrix[3,3] N_sim = rep_matrix(0, 3, 3);
for (n in 1:N){
real y_latent = s_latent * z_latent[n];
matrix[3,2] p;
for (j in 1:2) {
p[1,j] = normal_cdf(2.5 | y_latent + mu_obs[j], s_obs[j]);
p[2,j] = normal_cdf(4.5 | y_latent + mu_obs[j], s_obs[j]) - p[1,j];
p[3,j] = 1 - (p[1,j] + p[2,j]);
}
N_sim += p[,1] * p[,2]';
}
N_sim = a*N_sim;
}
model {
mu_obs ~ normal(3,3);
s_latent ~ exponential(1);
s_obs ~ exponential(1);
sigma ~ exponential(1);
for (j1 in 1:3) {
for (j2 in 1:3) {
sqrt(N_table[j1, j2]) ~ normal(sqrt(N_sim[j1, j2]), sigma);
}
}
}
You can see from the last line of code that I also changed the likelihood. The Poisson likelihood is appropriate if the model is correct, but it won't be correct; by fitting on the square root scale but including a free parameter, sigma, we're allowing for overdispersion.
We call it from R:
fit_2 <- table_match_2$sample(data=table_data, parallel_chains=4, refresh=0)
print(fit_2, max_rows=30)
Aaahhhh, much better:
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -25.68 -25.29 2.22 2.08 -29.93 -22.80 1.00 1009 1687
mu_obs[1] 3.79 3.79 0.08 0.07 3.67 3.92 1.00 2090 2047
mu_obs[2] 3.66 3.66 0.07 0.06 3.56 3.78 1.00 1762 2092
s_latent 0.68 0.68 0.06 0.05 0.59 0.77 1.00 2562 1803
s_obs[1] 0.51 0.51 0.08 0.07 0.39 0.64 1.00 2562 1894
s_obs[2] 0.27 0.28 0.11 0.10 0.08 0.44 1.00 1791 1139
sigma 4.22 4.07 1.01 0.93 2.84 6.14 1.00 1511 2172
N_sim[1,1] 338.69 326.65 114.39 103.68 173.55 539.20 1.00 2433 2009
N_sim[2,1] 434.37 418.15 164.48 152.14 199.93 722.66 1.00 2775 2271
N_sim[3,1] 1.63 0.52 3.59 0.70 0.01 6.71 1.00 2207 2454
N_sim[1,2] 588.17 570.14 193.09 180.59 304.30 926.40 1.00 2580 2098
N_sim[2,2] 9734.95 9745.35 520.62 488.97 8858.11 10563.91 1.00 3395 2533
N_sim[3,2] 1637.06 1621.71 322.01 300.94 1134.19 2184.29 1.00 3159 2436
N_sim[1,3] 0.66 0.13 1.92 0.18 0.00 2.98 1.00 2309 2220
N_sim[2,3] 595.11 575.31 193.51 175.52 319.79 954.80 1.00 2496 2634
N_sim[3,3] 1323.36 1298.78 300.80 278.37 863.93 1847.16 1.00 1967 1755
Although, if you look at the expected values in the table (N_sim), it's not so great! The [1,3] and [3,1] values are 0.66 and 1.63, much lower than the data values of 40 and 51.
I think the problem is with the Poisson or square-root likelihood, that it allows relatively large error for small-count cells.
4. Stan program: another revision
So I'll try a new model, this time with a lognormal likelihood. This will go into the file table_match_3.stan. I'll just show you the one line in the model block because that's the only part of the program that's changed:
log(N_table[j1, j2]) ~ normal(log(N_sim[j1, j2]), sigma);
We call it from R:
fit_3 <- table_match_3$sample(data=table_data, parallel_chains=4, refresh=0)
print(fit_3, max_rows=30)
Which yields:
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -6.98 -6.42 3.35 3.34 -13.32 -2.51 1.01 339 559
mu_obs[1] 4.12 4.00 0.59 0.37 3.48 5.25 1.00 819 721
mu_obs[2] 3.86 3.74 0.52 0.22 3.38 4.89 1.01 761 600
s_latent 0.82 0.69 0.48 0.22 0.40 1.69 1.01 434 306
s_obs[1] 0.98 0.95 0.37 0.28 0.46 1.58 1.00 1040 744
s_obs[2] 0.55 0.51 0.34 0.23 0.09 1.16 1.01 777 833
sigma 0.84 0.74 0.38 0.30 0.42 1.56 1.01 448 810
N_sim[1,1] 488.28 286.72 716.36 190.52 91.35 1536.07 1.00 1071 888
N_sim[2,1] 666.97 552.96 480.97 294.04 220.03 1474.32 1.00 1772 1300
N_sim[3,1] 115.33 64.10 308.25 42.04 19.43 302.69 1.01 1794 1046
N_sim[1,2] 1043.18 944.89 607.16 461.48 313.54 2133.92 1.00 1109 1153
N_sim[2,2] 5820.30 6119.15 2055.18 2028.03 1832.99 8729.56 1.01 432 610
N_sim[3,2] 2887.52 2773.01 1308.13 1259.89 941.96 5226.20 1.00 1142 1137
N_sim[1,3] 69.16 36.91 145.99 24.02 12.50 203.99 1.00 1730 1098
N_sim[2,3] 1139.58 969.93 709.41 479.43 390.12 2416.68 1.00 1478 1611
N_sim[3,3] 2423.68 1726.31 2045.37 1221.81 495.16 7118.74 1.01 526 635
Or, in optimize setting:
fit_optimize <- table_match$optimize(data=table_data)
print(fit_3_optimize, max_rows=30)
variable estimate
lp__ 1.62
mu_obs[1] 3.86
mu_obs[2] 3.68
s_latent 0.54
s_obs[1] 0.82
s_obs[2] 0.52
sigma 0.38
N_sim[1,1] 209.05
N_sim[2,1] 549.51
N_sim[3,1] 49.04
N_sim[1,2] 931.32
N_sim[2,2] 8175.10
N_sim[3,2] 2736.90
N_sim[1,3] 30.44
N_sim[2,3] 984.07
N_sim[3,3] 988.56
Not bad! Not perfect---you can see that the [3,3] cell is much too big. But it'll do.
5. Conclusion
The model doesn't quite fit! Maybe it could make sense to try longer-tailed distributions such as the t with low degrees of freedom. But for our original purpose, which is to understand the impact of measurement error on a nonlinear relationship, this will be just fine. Also we've learned something about statistical computing.