Pre- and post-treatment analysis with measurement errors

Modeling
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.

1 Like