Phil,
This is great! I had no idea nls would accept functions in the formula
position. My apologies for not including data to reproduce my issue.
dat<-data.frame(X=c(-13.0,-11.0,-10.0,-9.5,-9.0,-8.5,-8.0,-7.5,-7.0,-6.5,-6.0,-5.0,
-13.0,-11.0,-10.0,-9.5,-9.0,-8.5,-8.0,-7.5,-7.0,-6.5,-6.0,-5.0),
Y=c(3146,3321,2773,2415,2183,1091,514,191,109,65,54,50,
3288,3243,2826,2532,2060,896,517,275,164,106,202,53))
With your suggestion, I've changed the formula in nls to the following
function:
myfunc<-function(NS,LogKi,BMax)with(dat,{
KdCPM = KdnM*SpAct*Vol*1000
R<-NS+1
S<-(1+10^(X-LogKi))*KdCPM+Hot
a<-(-1*R)
b<-R*S+NS*Hot+BMax
c<--1*Hot*(S*NS+BMax)
(-1*b+sqrt(b*b-4*a*c))/(2*a)
})
But to get it to compute without errors, I also had to increase the
tolerance level: the step factor keeps being reduced below the min
factor. Looking at the trace of the nls though, I don't see any changes
after the 10th iteration or so; would increasing the tolerance cause any
issue that I'm not thinking of?
KdnM <- .8687
SpAct <- 4884
Vol <- .125
Hot <- 10191.0
nls(Y~myfunc(NS,LogKi,BMax),data=dat,start=list(NS=.01,LogKi=-7,BMax=10*max(dat['Y'])),control=nls.control(warnOnly=TRUE,minFactor=1e-5),trace=TRUE)
Also, I've found that if the start value I provide for BMax is too
inaccurate (ex. max(dat['Y']), nls generates the 'singular gradient' error
message, which isn't something I'm used to. Usually nls is kind enough to
inform me that the initial parameter estimates are what caused the problem.
Has the error message changed in a recent update, or is this a different
error message than what I'm thinking about?
Thanks again for all your help,
Jared
On Mon, Dec 13, 2010 at 1:23 PM, Phil Spector <spec...@stat.berkeley.edu>wrote:
> Jared -
> nls will happily accept a function on the right hand side
> of the ~ -- you don't have to write out the formula in such
> detail.
> What you provided isn't reproducible because you didn't provide data, and
> it's not clear what "Y" in the formula
> represents. Let me provide you with an admittedly simpler
> reproducible example.
>
> Suppose we want to estimate the model
>
> response = v * dose / (k + dose)
>
> where response and dose are variables in a data frame called "dat",
> and v and k are the parameters to be estimated.
>
> Here's the data:
>
> dat = data.frame(dose=c(0.027,0.044,0.073,0.102,0.175,0.257,0.483,0.670),
>>
> + response=c(12.7,16.0,20.4,22.3,26.0,28.2,29.6,31.4))
>
> Normally we would fit such a simple model as
>
> nls(response ~ v*dose / (k + dose),data=dat,start=list(v=30,k=.05))
>>
> Nonlinear regression model
> model: response ~ v * dose/(k + dose)
> data: dat
> v k 32.94965 0.04568
> residual sum-of-squares: 1.091
>
> Number of iterations to convergence: 4 Achieved convergence tolerance:
> 8.839e-07
>
> Alternatively, I can write a function like this:
>
> myfunc = function(v,k)with(dat,v * dose / (k + dose))
>>
>
> and use the following call to nls:
>
> nls(response ~ myfunc(v,k),data=dat,start=list(v=30,k=.05))
>>
> Nonlinear regression model
> model: response ~ myfunc(v, k)
> data: dat
> v k 32.94965 0.04568
> residual sum-of-squares: 1.091
>
> Number of iterations to convergence: 4 Achieved convergence tolerance:
> 8.839e-07
>
> which gets the identical results.
>
> Admittedly this function is trivial, but perhaps in your case
> you could use the segments from prism to construct a function
> for the right-hand side of your nls call.
>
> Hope this helps.
> - Phil Spector
> Statistical Computing Facility
> Department of Statistics
> UC Berkeley
> spec...@stat.berkeley.edu
>
>
>
>
>
>
> On Mon, 13 Dec 2010, Jared Blashka wrote:
>
> I'm attempting to calculate a regression in R that I normally use Prism
>> for,
>> because the formula isn't pretty by any means.
>>
>> Prism presents the formula (which is in the Prism equation library as
>> Heterologous competition with depletion, if anyone is curious) in these
>> segments:
>>
>> KdCPM = KdnM*SpAct*Vol*1000
>> R=NS+1
>> S=(1+10^(X-LogKi))*KdCPM+Hot
>> a=-1*R
>> b=R*S+NS*Hot+BMax
>> c = -1*Hot*(S*MS+BMax)
>> Y = (-1*b+sqrt(b*b-4*a*c))/(2*a)
>>
>> I'm only trying to solve for NS, LogKi, and BMax. I have everything else
>> (KdnM, SpAct, Vol, Hot).
>>
>> I would use the simple formula at the bottom and then backsolve for the
>> terms I'm looking for, but the simple formula at the bottom takes out the
>> X
>> term, which is contained within S, which it itself contained in both b and
>> c.
>> So I tried to substitute all the terms back into Y and got the following
>>
>> formula<-as.formula("Y ~
>>
>> (-1*(((NS+1)*((1+10^(X-LogKi))*(KdnM*SpAct*Vol*1000)+Hot))+NS*Hot+BMax)+sqrt((((NS+1)*((1+10^(X-LogKi))*(KdnM*SpAct*Vol*1000)+Hot))+NS*Hot+BMax)*(((NS+1)*((1+10^(X-LogKi))*(KdnM*SpAct*Vol*1000)+Hot))+NS*Hot+BMax)-4*(-1*(NS+1))*(-1*Hot*(((1+10^(X-LogKi))*(KdnM*SpAct*Vol*1000)+Hot)*NS+BMax))))/(2*-1*(NS+1))")
>>
>> But trying to use that formula gives me the single gradient message, which
>> I
>> wasn't entirely surprised by.
>> fit<-nls(formula=formula,data=data,start=list(NS=.01,LogKi=-7,BMax=33000))
>> Error in nls(formula = formula, data = all_no_outliers, start = list(NS =
>> 0.01, :
>> singular gradient
>>
>> I've never worked with a formula this complicated in R. Is it even
>> possible
>> to do something like this? Any ideas or points in the right direction
>> would
>> be greatly appreciated.
>>
>> Thanks,
>> Jared
>>
>> [[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.
>>
>>
[[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.
