Thanks Ted! I really appreciate your time!
Thanks for the link about the 'problem of calibration', and your suggestion to reformulate my model. I had no idea about it before. I certainly learnt something today. I will try your suggestions later today and let you know how it works out. Stefan -----Original Message----- From: r-help-boun...@r-project.org on behalf of ted.hard...@wlandres.net Sent: Fri 11/11/2011 4:07 PM To: r-help@r-project.org Subject: Re: [R] Predicting x from y Follow-up: See at end. On 11-Nov-11 21:16:02, Ted Harding wrote: > On 11-Nov-11 14:51:19, Schreiber, Stefan wrote: >> Dear list members, >> >> I have just a quick question: >> >> I fitted a non-linear model y=a/x+b to describe my data >> (x=temperature and y=damage in %) and it works really nicely >> (see example below). I have 7 different species and 8 individuals >> per species. I measured damage for each individual per species >> at 4 different temperatures (e.g. -5, -10, -20, -40). >> Using the individuals per species, I fitted one model per species. >> Now I'd like to use the fitted model to go back and predict the >> temperature that causes 50% damage (and it's error). Basically, >> it pretty easy by just rearranging the formula to x=a/(y-b). >> But that way I don't get a measure of that temperature's error, >> do I? Can I use the residual standard error that R gave me for >> the non-linear model fit? Or do I have to fit 8 lines (each >> individual) per species, calculate x based on the 8 individuals >> and then take the mean? >> >> Unfortunately, dose.p from the MASS package doesn't work for >> non-linear models. When I take the log(abs(x)) the relationship >> becomes not satisfactory linear either. >> >> Any suggestions are highly appreciated! >> >> Thank you! >> Stefan >> >> EXAMPLE for species #1: >> >> y.damage<-c(5.7388985,1.7813519,3.7321461,2.9671031, >> 0.3223196,0.3207941,-1.4197658,-5.3472160, >> 41.1826677,29.3115243,31.3208841,35.3934115, >> 58.5848778,31.1541049,42.2983479,27.0615648, >> 64.1037728,54.7003353,66.7317044,65.4725881, >> 72.5755056,67.2683495,64.8717942,65.9603073, >> 75.0762273,56.7041960,60.0049429,70.0286506, >> 73.2801947,72.7015642,75.0944694,81.0361280) >> >> x.temp<-c(-5,-5,-5,-5,-5,-5,-5,-5,-10,-10,-10,-10,-10,-10,-10, >> -10,-20,-20,-20,-20,-20,-20,-20,-20,-40,-40,-40,-40,-40, >> -40,-40,-40) >> >> nls(y.damage~a/x.temp+b,start=list(a=400,b=80)) >> plot(y.damage~x.temp,xlab='Temperature',ylab='Damage [%]') >> curve(409.61/x+81.84,from=min(x.temp),to=max(x.temp),add=T) > > A couple of comments. > > First, in general it is not straightforward to estimate > the value of a covariate (here temperature) by inverting > the regression of a response (here damage) on that covariate. > This the "inverse regression" or "calibration" problem, > (and it is problematic)! For instance, in linear regression > the estimate obtained by inversion has (theoretically) > no expectation, and has infinite variance. For an outline, > and a few references, see the Wikipedia article: > > http://en.wikipedia.org/wiki/Calibration_(statistics) > > Second, I would be inclined to try nls() on a reformulated > version of the problem. Let T50 denote the temperature for > 50% damage, and introduce this as a parameter (displacing > your parameter "a"): > > y = 50*(b + T50)/(b + x) > > where T50 = a/50 - b in terms of your original parameters > "a" and "b". With this formula for the non-linear dependence > of damage on temperature, it is no longer necessAry to invert > the regression equation, since the parameter you want is > already there and will be estimated directly. > > Hoping this helps, > Ted. I think I have mis-read your model: I read it as y = a/(x+b) whereas you wrote "y/x+b" and your model formula in nls() is y.damage~a/x.temp+b, i.e. y.damage ~ (a/x.temp) + b which confirms it. In that case, you may be able to get a satisfactory result by using a linear regression with y.damage = a*z.temp + b where z.temp = 1/x.temp so the model formula would be y.damage ~ z.temp You then have the straightforward inverse regression problem (aka calibration problem). The solution to this takes a bit of explanation, which I do not have the time for right now. I will write further about it in the morning. Ted. -------------------------------------------------------------------- E-Mail: (Ted Harding) <ted.hard...@wlandres.net> Fax-to-email: +44 (0)870 094 0861 Date: 11-Nov-11 Time: 23:07:35 ------------------------------ XFMail ------------------------------ ______________________________________________ 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. [[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.
