Dave said my newdata data frame 'new' must have a column named 'area'.
It did. Nonetheless predict.lm throws an error with type = "terms" and
newdata = new. I see nothing in the predict.lm documentation that
bars this usage. Is there a bug?
To illustrate an OLS behavior, I had cited Ludbrook '12. Dave noted
the high reprint fee. So, from that paper,
#Ludbrook's data set S1 (except renaming
#his 'x' as 'concn' and his 'y' as 'area')
S1 <- data.frame(
area = c(2.4,2.6,6.0,6.5,8.9,),
concn = c(1.1,4,5,8.5,8.5))
#Regular and inverse OLS regressions
reg <- lm(area ~ concn, data = S1)
inv.reg <- lm(concn ~ area, data = S1)
#View the regular fit (blue line)
plot(area ~ concn, data = S1)
abline(reg, col = "blue") #OLS
#Add the inverse fit (red line)
lowHiAreas <- data.frame(area = c(0,8))
lowHiConcs <- predict(inv.reg,
newdata = lowHiAreas,
type = "response")
lines(x = lowHiConcs,
y = unlist(lowHiAreas), col = "red")
According to Ludbrook the intersection of these lines is at point
{mean(concn), mean(area)}. He also shows several Model II
regression lines (e.g., least product regression) that also intersect there.
Dave and Rui both note that the concentrations I wish to predict will
include error, whereas my calibration concentrations are fixed. This
issue was raised also by VitoshKa in 2010
http://stackoverflow.com/questions/4191857/inverse-of-predict-function-in-r
who suggested that a Bayesian approach may be called for.
Rui showed how to convert normalized/centered predicted concentrations
back to real ones.
Rui (and also James in the Stack Overflow thread) talked about using
the x' = (y' - b0)/b1 approach. James did it in a way that he said
preserved variable names.
The Stack Overflow discussion also mentioned the function
chemCal::inverse.predict(), which is based on Massart et al., '97. I
need to read about this.
I guess for now I'll use the x' = (y' - b0)/b1 approach.
As no one has explained why predict.lm fails with type = "terms"
and newdata, and since I'm pretty sure I successfully did this in
years past, I may document this and submit it as a possible bug.
Thanks all,
-John
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
