# Tolerance interval for Deming regression

I am trying to derive (one-sided) tolerance intervals related to the Deming regression model: $$x_i=x^*_i + \epsilon_i$$ $$y_i = (\alpha+\beta x^*_i) + \epsilon'_i$$ where the $x^*_i$'s are nonrandom fixed numbers, $\epsilon_i \sim {\cal N}(0,\sigma_x^2)$, $\epsilon'_i \sim {\cal N}(0,\sigma_y^2)$, and all variables $\epsilon_i, \epsilon'_i$ are mutually independent. The ratio $\sigma_x^2/\sigma_y^2$ is assumed to be known. Put $\tau^2=\sigma_x^2+\sigma_y^2$.

Given a new true $x$'' value $x^*_\text{new}$ I'm looking for an upper tolerance bound of the difference $$y_\text{new} - x_\text{new} \sim {\cal N}(\alpha+(\beta-1)x^*_\text{new}, \tau^2).$$ Actually this upper tolerance bound is nothing but an upper $100(1-\alpha)\%$-confidence bound about the $100p\%$-quantile of $y_\text{new} - x_\text{new}$. As for classical Gaussian models I'm looking for a tolerance bound having form $$B(x,y)=(\hat{y^*}_\text{new} - x^*_\text{new}) + K s$$ with $\hat{y^*}_\text{new} = \hat\alpha + \hat\beta x^*_\text{new}$ and $s^2= \hat\tau^2$. Denoting by $F_\text{new}$ the cumulative distribution function of the difference $y_\text{new} - x_\text{new}$, the tolerance factor $K$ must be derived in order that $$\Pr\left( F_\text{new}(B(x,y)) > p \right) = 1-\alpha$$ for some given $p\in]0,1[$ (the tolerance level) and confidence level $1-\alpha$ (of course here $\alpha$ is not the intercept parameter $\alpha$ of the model!)

Estimates of $\alpha$, $\beta$ and $\tau^2$ as well as their asymptotic variance matrix $V$ and an estimate $\hat V$ of this asymptotic variance matrix are available in some textbooks and/or research papers (see for instance Fuller's book or Iles & Gillard's report 1 and report2). In particular $\hat\tau^2$ is asymptotically independent of $(\hat\alpha,\hat\beta)$.

Now I'm trying to derive the tolerance bound by miming the way used for deriving the tolerance bound in a classical regression model. It is based on the following lemma.

Lemma Let $W \sim {\cal N}(0, \sigma^2)$ independent of $Q \sim \frac{\chi^2_{d}}{d}$. Then the number $K>0$ satisfying $$\Pr\left( \Phi(W+K\sqrt{Q}) > p \right) =1-\alpha$$ is given by a known formula $K=K(\sigma,d,p,\alpha)$ (see Krishnamoorthy & Mathew's book).

Now put $X_\text{new}=(1,x_\text{new})'$. One has the three asymptotically independent random variables $$\frac{(y_\text{new} - x_\text{new} ) - (\alpha+(\beta-1)x^*_\text{new})}{\sqrt{\tau^2}} \sim {\cal N}(0,1),$$ $$\frac{(\hat{y^*}_\text{new} - x^*_\text{new} ) - (\alpha+(\beta-1)x^*_\text{new})}{\sqrt{\tau^2}} \sim {\cal N}\left(0,\frac{X_\text{new}'VX_\text{new}}{\tau^2}\right),$$ $$\frac{s^2}{\tau^2} \sim \frac{\chi^2_{n-2}}{n-2}$$

The definition of $K$ becomes $$\Pr\left( \Phi(W+K\sqrt{Q}) > p \right) =1-\alpha$$ with $W \sim {\cal N}\left(0,\frac{X_\text{new}'VX_\text{new}}{\tau^2}\right)$ independent of $Q \sim \frac{\chi^2_{n-2}}{n-2}$.

Applying the lemma we get $K$ depending on the unknown values of $\tau$ and $V$. Substituting with the estimates $\hat\tau$ and $\hat V$ gives an approximate tolerance factor $\hat K$.

Well, I have done it and I have implemented everything in R. But assessing the performance of my tolerance bound with simulation shows that even with large sample sizes it does not work, except for the case $\beta=1$. I'm pretty sure of my formulas for the estimates and the variance matrices since I have also implemented confidence intervals and prediction intervals and they work well.

EDIT: In fact that works !! There was a silly error in the R code I used for simulations. Should I delete this post ? I think no because the method could interest others.

-
Hi Stephane: if you are able to post an answer to your own question with an illustrative example and/or code, I think that would make an excellent addition to the site! :-) –  cardinal Nov 26 '12 at 16:08
@cardinal Ok. I will do it whenever I will be less busy. –  Stéphane Laurent Dec 3 '12 at 19:43
Stephane: Sounds good. I know how that goes! :-) –  cardinal Dec 3 '12 at 19:56
Done, @cardinal. –  Stéphane Laurent Jul 1 '13 at 21:18

Below is the code. Not very pretty, but (as said in the update of my question) it works well.

# returns estimates
deming.estim <- function(x,y,lambda=1){  # lambda=sigmay²/sigmax²
n <- length(x)
my <- mean(y)
mx <- mean(x)
SSDy <- crossprod(y-my)[,]
SSDx <- crossprod(x-mx)[,]
SPDxy <- crossprod(x-mx,y-my)[,]
A <- sqrt((SSDy - lambda*SSDx)^2 + 4*lambda*SPDxy^2)
B <- SSDy - lambda*SSDx
beta <- (B + A) / (2*SPDxy)
alpha <- my - mx*beta
sigma.uu <- ( (SSDy + lambda*SSDx) - A ) /(2*lambda) / (n-1)
s.vv <- crossprod(y-my-beta*(x-mx))/(n-2)
# formula Gilard & Iles
sbeta2. <- (SSDx*SSDy-SPDxy^2)/n/(SPDxy^2/beta^2)
sbeta. <- sqrt(sbeta2.)
salpha2. <- s.vv/n + mx^2*sbeta2.
salpha. <- sqrt(salpha2.)
V <- rbind( c(salpha2., -mx*sbeta2.), c(-mx*sbeta2., sbeta2.) )
return(list(alpha=alpha,beta=beta,
V=V,
sigma=sqrt(sigma.uu*(n-1)/(n-2)))
)
}

# returns one-sided upper tolerance bound
deming.tolerance <- function(x,y,lambda=1,xnew,p=98/100,alph=5/100){
fit <- deming.estim(x,y,lambda)
V <- fit$V sigma.uu <- fit$sigma^2
sigma.ee <- lambda*sigma.uu
Vbeta.over.sigma <- V/(sigma.ee + sigma.uu)
Xnew <- as.matrix(c(1,xnew))
d <- sqrt( (t(Xnew)%*%Vbeta.over.sigma%*%Xnew) )
k1 <- d * qt(1-alph, length(x)-2, ncp=qnorm(p)/d)
S <- sqrt(sigma.ee + sigma.uu)
ynew <- fit$alpha+fit$beta*xnew
TolUp <- ynew - xnew + k1*S
TolUp
}


Some simulations to check the frequentist coverage:

#################
## simulations ##
#################
n <- 60
mu <- runif(n, 1, 100)
xnew <- 50
sigma.x <- 1
lambda0 <- 2
sigma.y <- sqrt(lambda0)*sigma.x
alpha <- 2
beta <- 3
p <- 0.98
# true value of the quantile:
qq <- qnorm(p, alpha+(beta-1)*xnew, sqrt(sigma.x^2+sigma.y^2))

nsims <- 5000
test <- rep(NA,nsims)
lambda <- lambda0
for(i in 1:nsims){
x <- rnorm(n, mu, sigma.x)
y <- rnorm(n, alpha+beta*mu, sigma.y)
test[i] <- ( qq < deming.tolerance(x,y,lambda=lambda,xnew=xnew,p=p) )
}

> mean(test)
[1] 0.9544

-