# Pre- and post-treatment analysis with measurement errors

#1

I am looking at pre- and post-treatment data, x1 and x2, respectively, with paired measurements x1,obs and x2,obs that have a normal measurement error. Physically, x2 ≥ x1 must hold for each subject; however, because of measurement errors, some observed x2’s are less than the corresponding x1’s. Right now my model is just a simple linear regression with measurement errors, with the slope and intercept constrained to be positive. However, I’m not sure how to incorporate the x2 ≥ x1 constraint into the model. Any suggestions? Thanks!

``````data {
int<lower=0> N;	//	number of observations
vector[N] x_obs; // x-measurement, observed with error
vector[N] y_obs; // y-measurement, observed with error
real<lower=0> ub;	// upper bound of measurements
real<lower=0> tau_b; // beta parameter for gamma distribution

// Data for prediction
int<lower=0> N_hat;
vector[N_hat] x_hat;

int method;	// 0 for Normal likelihood, 1 for Student
}

parameters {
// Latent model
real<lower=0> tau;	// measurement error

vector<lower=0, upper=ub>[N] x;	// unknown true x
real<lower=0, upper=ub> mu_x;	// mean of unknown true x
real<lower=0> sigma_x;	// SD of unknown true x

// Difference model
real<lower=0> b0;
real<lower=0> b1;
real<lower=0, upper=ub> sigma_y;
real<lower=0> nu;	// Student likelihood
}

transformed parameters {
vector<lower=0>[N] mu;
mu = b0 + b1*x;
}

model {
// Latent model
tau ~ gamma(1, tau_b);
x ~ normal(mu_x, sigma_x);
x_obs ~ normal(x, tau);

// Regression model
b0 ~ normal(0, 3);
b1 ~ normal(1, 3);
sigma_y ~ gamma(2, 40);
nu ~ gamma(2, 0.1);

if (method==0)
y_obs ~ normal(mu, sigma_y);
else if (method==1)
y_obs ~ student_t(nu, mu, sigma_y);
}

generated quantities {
// PPC
vector[N] y_rep;

// New predictions
vector[N_hat] y_hat;
vector[N_hat] mu_hat;
mu_hat = b0 + b1*x_hat;

for (i in 1:N) {
if (method==0)
y_rep[i] = normal_rng(mu[i], sigma_y);
else if (method==1)
y_rep[i] = student_t_rng(nu, mu[i], sigma_y);
}

for (i in 1:N_hat) {
if (method==0)
y_hat[i] = normal_rng(mu_hat[i], sigma_y);
else if (method==1)
y_hat[i] = student_t_rng(nu, mu_hat[i], sigma_y);
}
}``````

#2

Sorry I missed this earlier. Great question!

You just want to use `ordered` as the variable type for the true value.

``````data {
vector[2] y_obs;
parameters {
ordered[2] y;  // ensures y[1] <= y[2]

model {
y_obs ~ normal(y, sigma);
``````

It’s also sensible to throw a prior on `y` in addition to the ordering constraint.