[R] calibration/inverse regression?

David Lucy dlucy at maths.ed.ac.uk
Fri Mar 15 18:00:35 CET 2002


Bill,

These should work OK - I did it a long time ago so the code is a little
crude. Watch out for some of the longer lines which may not wrap well in
ascii text.

Dave.


> I wonder if anyone out there has written a routine to solve the simple
> linear calibration problem?
> 
> - fit regression of y vs x
> - estimate the value x0 (with 95% CI) that gives y0


# R function to do classical calibration
# input the x and y vectors and the value of y to make an
# estimate for
# val is a value of dep
# returns the calibrated value of x and it's
# single confidence interval about x
# form of calibration particularly suited to jackknifing
calib <- function(indep, dep, val)
{

	# generate the model y on x and get out predicted values for x
	reg <- lm(dep~indep)
	xpre <- (val - coef(reg)[1])/coef(reg)[2]

	# generate a confidence
	yyat <- ((dep - predict(reg))^2)
	sigyyat <- ((sum(yyat)/(length(dep)-2))^0.5)
	term1 <- sigyyat/coef(reg)[2]
	sigxxbar <- sum((indep - mean(indep))^2)
	denom <- sigxxbar * ((coef(reg)[2])^2)
	conf <- abs(((1+(1/length(dep))+(((val -
mean(dep))^2)/denom))^2)*term1)

return(xpre, conf)
}


# R function to do classical calibration
# input the x and y vectors
# returns a data.frame defined in the main function with
# vectors xpre and conf which are the predicted x values and
# the single confidence interval about x
calib <- function(indep,dep)
{

	# generate the model y on x and get out predicted values for x
	reg <- lm(dep~indep)
	xpre <- (dep-coef(reg)[1])/coef(reg)[2]

	# generate a confidence
	yyat <- ((dep-predict(reg))^2)
	sigyyat <- ((sum(yyat)/(length(dep)-2))^0.5)
	term1 <- sigyyat/coef(reg)[2]
	sigxxbar <- sum((indep-mean(indep))^2)
	denom <- sigxxbar * ((coef(reg)[2])^2)
	conf <- 
abs(((1+(1/length(dep))+(((dep-mean(dep))^2)/denom))^2)*term1)

return(xpre, conf)
}






**********************************************************************
**  Dr. David Lucy                                                  **
**  Joseph Bell Centre for Forensic Statistics and Legal Reasoning  **
**  Department of Mathematics and Statistics                        **
**  The University of Edinburgh                                     **
**  James Clerk Maxwell Building                                    **
**  King's Buildings                                                **
**  Mayfield Road                                                   **
**  Edinburgh                                                       **
**  EH9 3JZ                                                         **
**                                                                  **
**  tel: 0131 650 5057                                              **
**  extension: 505057                                               **
**  e-mail: dlucy at maths.ed.ac.uk                                    **
**                                                                  **
**  Truncated address:                                              **
**                                                                  **
**  Dr. David Lucy                                                  **
**  Department of Mathematics                                       **
**  JCMB                                                            **
**  King's Buildings                                                **
**  Edinburgh University                                            **
**  Edinburgh                                                       **
**  EH9 3JZ                                                         **
**********************************************************************



-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list