[R] Predicting x from y
(Ted Harding)
ted.harding at wlandres.net
Sat Nov 12 00:07:38 CET 2011
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.harding at wlandres.net>
Fax-to-email: +44 (0)870 094 0861
Date: 11-Nov-11 Time: 23:07:35
------------------------------ XFMail ------------------------------
More information about the R-help
mailing list