# Bayesian estimates for Deming regression coinciding with least-squares estimates

Consider the following Deming model with independent replicates :

$$x_{i,j} \mid \theta_{i} \sim {\cal N}(\theta_{i}, \gamma_X^2), \quad y_{i,j} \mid \theta_{i} \sim {\cal N}(\alpha+\beta\theta_{i}, \gamma_Y^2), \\ i \in 1:I, \quad j \in 1:J,$$ all observations being independent, all parameters being unknown.

Let's fix the parameters and the design and write a function that simulates the data:

N <- 25  # number of groups
J <- 4  # number of repetitions within each group
theta0 <- 1:N  # true theta_i
gammaX0 <- sqrt(4);
lambda <- 1
gammaY0 <- sqrt(lambda)*gammaX0
alpha <- 0
beta <- 1
simulations <- function(){
X <- Y <- array(NA,dim=c(N,J))
for(i in 1:N){
for(j in 1:J){
X[i,j] <- rnorm(1,theta0[i],gammaX0)
Y[i,j] <- rnorm(1,alpha+beta*theta0[i],gammaY0)
}
}
return(list(X=X,Y=Y))
}
set.seed(666)
sims <- simulations()
X <- sims$X; Y <- sims$Y


When the ratio $\lambda=\frac{\gamma_Y^2}{\gamma_X^2}$ is known, one can get good estimates of $\alpha$ and $\beta$ by considering only the group means $$\bar{x}_{i\bullet} \sim {\cal N}(\theta_i, \gamma_X^2/J), \quad \bar{y}_{i\bullet} \sim {\cal N}(\alpha + \beta\theta_i, \gamma_Y^2/J)$$ and then by taking the maximum likelihood estimates:

MethComp::Deming(rowMeans(X), rowMeans(Y), vr=lambda)[1:2]
##  Intercept      Slope
## -0.3207636  1.0154527


On the other hand, the least-squares estimates:

coefficients(lm(rowMeans(Y)~rowMeans(X)))
## (Intercept) rowMeans(X)
##  0.04505284  0.98709755


are biased and unconsistent.

Now, using a Bayesian approach with the following naive and independent priors: $$\gamma_X^2 \sim {\cal IG}(0^+,0^+), \qquad \gamma_Y^2 \sim {\cal IG}(0^+,0^+), \\ \theta_i \sim {\cal N}(0, \infty), \\ \alpha \sim {\cal N}(0, \infty), \qquad \beta \sim {\cal N}(0, \infty),$$ then the Bayesian estimates (posterior means or medians) of $\alpha$ and $\beta$ are close to the least-squares estimates. And I wonder why.

I know there is no reason to get a good coincidence between maximum likelihood estimates and Bayesian estimates by using "naive" non-informative priors, but I am surprised by the coincidence between the least-squares estimates and the Bayesian estimates. Is there an intuitive reason to expect this result ?

This is what I discovered using JAGS:

library(rjags)
modelfile <- "gibbs1.txt"
jmodel <- function(){
for(i in 1:N){
for(j in 1:J){
X[i,j] ~ dnorm(theta[i], inv.gammaX2)
Y[i,j] ~ dnorm(nu[i], inv.gammaY2)
}
theta[i] ~ dnorm(0, 0.00001)
nu[i] <- alpha+beta*theta[i]
}
alpha ~ dnorm(0,0.001)
beta ~ dnorm(0,0.001)
#
inv.gammaX2 ~ dgamma(0.01,0.01)
inv.gammaY2 ~ dgamma(0.01,0.01)
gammaX2 <- 1/inv.gammaX2
gammaY2 <- 1/inv.gammaY2
lambda <- gammaY2/gammaX2
}
writeLines(paste("model", paste(body(jmodel), collapse="\n"), "}"), modelfile)
# run Gibbs sampler
data <- list(X=X, Y=Y, N=N, J=J)
inits1 <- list(alpha=0, beta=1, inv.gammaX2=1, inv.gammaY2=1, theta=rowMeans(X))
inits2 <- lapply(inits1, function(x) x*rnorm(length(x),1,.01))
inits <- list(inits1,inits2)
jags <- jags.model(modelfile,
data = data,
inits=inits,
n.chains = 2,
summary(samples)$statistics ## Mean SD Naive SE Time-series SE ## alpha -0.03753944 0.65340570 0.0046202760 0.016653968 ## beta 0.99340104 0.04402299 0.0003112895 0.001129354  Below are the results of$1000$simulations (I can provide the code for simulations) :   In fact, it is not difficult to derive the full conditional distributions in the Gibbs sampler. And one gets (I can provide details if needed): • the full conditional distribution of$\theta_i$is a Gaussian distribution and the expression of its mean as a function of$(\alpha, \beta)$is the same as the expression of the ML estimate of$\theta_i$as a function of the ML estimates$(\hat\alpha, \hat\beta)$assuming the known ratio$\lambda=\gamma^2_Y/\gamma_X^2$; • the full conditional distribution of$(\alpha, \beta)$is a Gaussian distribution centered around the ML ($=$least-squares) estimates for the ordinary regression of the$y_{ij}$vs the$\theta_i$(each$\theta_i$replicated$J$times). It is interesting to note that iterating the above procedure about the means, by assuming the true ratio at first step, yields the maximum likelihood estimates: alpha <- 0; beta <- 1 # initial values for(iter in 1:100){ theta <- rowMeans(X) + beta/(lambda+beta^2)*(rowMeans(Y)-(alpha+beta*rowMeans(X))) ab <- coefficients(lm(c(t(Y))~I(rep(theta, each=J)))) alpha <- ab; beta <- ab } ab ## (Intercept) I(rep(theta, each = J)) ## -0.3207636 1.0154527 MethComp::Deming(rowMeans(X), rowMeans(Y), vr=lambda)[1:2] ## Intercept Slope ## -0.3207636 1.0154527  In view of this fact I would be tempted to expect Bayesian estimates close to the maximum likelihood estimates, and not the least-squares estimates (corresponding to maximum likelihood estimates when$\lambda=\infty\$). What is wrong in this intuition ?