Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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.

share|improve this question
1  
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
1  
@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
1  
Done, @cardinal. –  Stéphane Laurent Jul 1 '13 at 21:18
add comment

1 Answer 1

up vote 1 down vote accepted

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
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.