EIV model of temporal trend with Berkson-like error

[roberto_r]

Hello,

I am trying to build an error-in-variables model in which the error is Berkson-like rather than a classical measurement error.

My goal is to recover the long-term trend of a continuous variable. To begin with, I am working with a simple linear regression (Github readme) so that I can understand the concepts before moving to more complex models.

For each simulated sample, I have:

a. A continuous measurement (for example, the length of a bone).
b. A time interval within which the sample could be dated (as is common in archaeology). For the purposes of this example, I assume that the error reflects the duration of the archaeological phase (rather than also thinking about the ‘real’ duration of each sample).

To recreate a realistic archaeological scenario (Github R Code: Script 1 and Script 2), I (1) simulate a true date by drawing it from a uniform distribution over the full study period (let’s say between 100-900). I then (2) partition the timeline into archaeological phases using a broken-stick Dirichlet process.

Each simulated sample is (3) assigned to the phase containing its true date, and its observed dating interval is taken to be the start and end dates of that phase.

I know this might seem unnecessary, but I then need to evaluate many possible phase configurations to see where the model performs best or worse. This phase assignment reflects the fact that a sample in archaeology is usually assigned to a chronological phase based on the context rather than an exact date (unless an object has been directly dated using scientific methods, for which the error is usually Gaussian). For instance, if the context has material from 300 to 500 and a bone is found in that layer, the bone will be dated uniformly between 300 and 500.

I fit two models: one with the EIV, and one where the true date is compared against the midpoint. In theory what I want to obtain is that the sample’s estimated date shifts from a uniform probability to a more skewed distribution. For instance if the linear regression shows that the value around year 500 is 5 and my sample (dated 300-500) is 4.8 I would think that the true date is closer to 500 than to 300. The midpoint date in this case would be 400.

However, somehow the intercept and slope for both models seem to be almost identical, even in terms of the uncertainty of the credible intervals. Am I designing my simulation wrong, the model, or this symmetry is to be expected and a midpoint model will work the same as the EIV model because the uniform distribution will push values towards the mean?

I am really sorry for the long question, I have tried to sort it out by doing some reading, AI checking/polishing the code, etc. but I am stuck on why mathematically the two models recover the parameters in a similar way, while I was expecting the midpoint one to perform worse.

I apologise in advance if I’ve misunderstood or mistated any statistical concept here — I don’t have formal training in statistics. I try to compensate working with manuals, exercises, forums like this one, and (alas!) occasionally LLMs. I am not sure ‘Berkson error’ is the right title here, or I should say structured non random error

[Bob_Carpenter]

Thanks for the careful question roberto_r.

I would have found the models much less confusing without all the transformations to start and end norms and mid norms. One issue your choice of normalization will run into is that you’re providing positive-only covariates. With linear regression, this induces a lot of correlation between the slope and intercept parameters. I’d suggest using [-1, 1] instead.

If you have a standard linear regression and then generalize it to a measurement error model where you assume the errors are symmetric around zero, you shouldn’t expect your posterior means to shift much. You can see this with a much simpler simulation than you’re running. The reason that uncertainties don’t shift much is a coincidence—you get more uncertainty from the unknown measurements, but then the posteriors for those can be more consistent with a linear regression and thus you can get less uncertainty in the regression itself. At least I think I’m thinking about this the right way.

If the measurements have bias to the low side, then your model isn’t really identified in how to deal with that. If they’re all shifted by the same absolute amount, then you can just soak this up in the intercept.

Did you intend the measurement error model to be uniform within that interval? That’s not going to let you model an underlyingly biased measurement error (e.g., all dates reported too early). It will let you infer them, but you don’t get the benefit of any partial pooling.

I didn’t trace through the simulation code as I have a hard time reading R, but did you simulate the data where the measurement error was not symmetric?

@roberto_r

Each simulated sample is (3) assigned to the phase containing its true date,

I would think this would also be a noisy process and you shouldn’t necessarily assume that it’s assigned the right phase. This is assuming that the phases are rather arbitrary cutoffs on continuously changing processes. You can adjust for this by letting the measurement-error boundaries go beyond the phase to which it’s assigned.

P.S. You can also make the Stan code faster by vectorizing.

parameters {
  vector<lower=0, upper=1>[N] true_date_raw;
}
transformed parameters {
  vector[N] true_date_norm = start_norm + true_date_raw .* (end_norm - start_norm);
}
model {
  alpha ~ normal(0, 10);
  beta  ~ normal(0, 10);
  sigma ~ exponential(1);
  y ~ normal(alpha + beta * true_date_norm, sigma);
}

If you don’t need true_date_norm in the output, you can define it as a local variable in the model before you define y, i.e.,

model {
  ...
  vector[N] true_date_norm = start_norm + true_date_raw .* (end_norm - start_norm);
  y ~ normal(alpha + beta * true_date_norm, sigma);

I would also advise against using “true” in a parameter name, since we assume parameters represent “true” parameters by default and the “true” will just clutter things up. If you have a variable alpha, just call it alpha and if you have a noisy observation of it, you can call that alpha_obs.

[roberto_r]

Hi Bob,

Thank you so much for this. I chose that year normalization mainly for two reasons:

1. In this example my timeline is only of 800 years, but in real life case studies the timeline could span many thousand years, and I thought that this would complicate the MCMC. I may be very wrong of course in this.
2. I was initially afraid, but I think it does not change much, that feeding negative years (years BC/BCE, prehistoric phases etc) would break this normalisation. But I will test the [-1;1] centering. In reality I have only tested this on positive years

@Bob_Carpenter

Did you intend the measurement error model to be uniform within that interval? That’s not going to let you model an underlyingly biased measurement error (e.g., all dates reported too early). It will let you infer them, but you don’t get the benefit of any partial pooling.

Yes, I believe that in absence of other information the ME can only be uniform, as in: we know that the sample belongs to a phase with Start_date and End_date, but we cannot assume a priori any shape. This brings me to the arbitrary phases: in reality they are arbitrary (and this is what I am trying to show). They can be historical periods (e.g. the Bronze Age), kingdoms (e.g. the reign of Augustus), pottery-phases, etc. and the result is that if you plot the start date and end date of these samples you get ‘blocks’ of samples with a hard boundary. I think a way ‘around’ this (even though it introduces a lot of complexity) would be to model the phase boundaries themselves, but what I am interested is mostly how what I am trying to model changes in calendar years, even though I only know phases (that can be very coarse at times).

Tomorrow I will:

1. Try to normalize the years with [-1;1]
2. Vectorize the Stan code (hopefully I can manage in someway)
3. Change the parameter’s name to something like date and date_raw ?

I forgot to mention something maybe important. The only difference (but maybe by construction) between the two models is sigma. In the midpoint model is always positively inflated, while the EIV model should capture the uncertainty already somehow and sigma remains low.

I was expecting however that, since a linear regression is just a line between points, if midpoints are ‘skewed’ when compared to the estimated dates the slope should at least change a bit. Is this the right way to think about it?

Again apologies for my statistical knowledge!

P.S. I use the true date in the output for the diagnostic plots that I show in the readme.

[Bob_Carpenter]

@roberto_r

feeding negative years

I wasn’t talking about BC/BCE, but rather just replacing (t - start) / (start - end), which ranges through (0, 1) with 2 * (t - start) / (start - end) + 1, which ranges through (-1, 1). It’d probably be even better if you just mapped x[n] to (x[n] - mean(x)) / sd(x).

@roberto_r

we cannot assume a priori any shape

It’s typical to fit a measurement error model, but these are often normal, as in the examples in User’s Guide section on measurement error. You could fit something that allowed skew if you anticipate that the measurements have asymmetric error.

@roberto_r

the result is that if you plot the start date and end date of these samples you get ‘blocks’ of samples with a hard boundary

I was suggesting widening the boundaries because you’re not going to be sure that a sample was classified into the right period, especially if arose somewhere near the boundary date. I don’t know this area, but am supposing that the end of the Bronze Age is fuzzy rather than being defined by a particular calendar date.

@roberto_r

The only difference (but maybe by construction) between the two models is sigma. In the midpoint model is always positively inflated, while the EIV model should capture the uncertainty already somehow and sigma remains low.

This isn’t surprising. The sigma goes down (but not necessarily the uncertainty in it), because the measurement error is encouraging the inferred data point from the measurement to be lower error.

@roberto_r

if midpoints are ‘skewed’ when compared to the estimated dates the slope should at least change a bit.

Not necessarily given the uniformity of the error you allow. If you simulate where the midpoint is known to be skewed low, you should be able to recover that.

Also keep in mind that the middle points in a regression have much less “leverage”.

[roberto_r]

Hi Bob,

Following up your advice:

1. I normalised the dates mapping them to [-1,1] and the model block is now vectored. Nothing really changed, but the sampling is more efficient (Neff is bigger) and faster. Also, the alpha/beta correlation dropped in both models (EIV and midpoint).
2. I tried more simulations varying how the timeline is divided, which phases get sampled, and where within a phase the dates actually fall.

When deposition is uniform within a phase, the two models end up with an almost identical slope regardless of how coarse the phases are or how many samples are in each phase. Sigma is still better for the EIV model, closer to the simulated value and wider credible intervals. In one example, EIV’s 90% sigma interval is 1.27–1.60 (Simulated value = 1.5), while midpoint’s is 1.98–2.27. The midpoint model slope’s CI is only very very modestly wider.

I then tried deposition that’s skewed within a phase, with more dates closer to one edge. On its own this does not shift the slope. The slope is only biased when I force wider phases towards one side of the timeline (at the beginning or at the end). The issue is that while I know in advance (in real life) the width of the phases on my timeline, I can’t know within a phase how samples are distributed. I guess this is why the EIV model is estimating them on average to be in the centre of the phase, which makes it equal to a midpoint model?

At this point I am not sure whether I should give up on my original idea, right now: the EIV model doesn’t give me a better trend estimate when deposition might be skewed. It only gives me a (loosely) ‘calibrated’ date for my sample, but - at least with a linear model (I’ll try changepoint after) - I am not getting what I was expecting (much wider credible intervals for the midpoint model).

One curiosity, when you suggested letting “the measurement-error boundaries go beyond the phase to which it’s assigned” to account for misclassification near a boundary, did you mean:

1. treating start_date/end_date themselves as uncertain, putting a prior on the phase boundary rather than treating it as a fixed year, or
2. widening the fixed window used in the uniform prior/likelihood by some fixed margin (e.g. Uniform(start - δ, end + δ)) without modelling the boundary as a parameter at all?

Thanks again for the detailed replies, this thread is helping me a lot!!