[R] Writing successful self-starting models for nlme
I R Cleasby
bop06irc at sheffield.ac.uk
Fri Sep 28 15:49:45 CEST 2007
Hi fellow R users,
My problem is the following:
I wish to use the nlme package to analyse the growth curves of house sparrow
chicks,
especially in relation to factors such as sex, position in the weight hierarchy
etc. as fixed effects and chick ID nested within BroodID as random effects.
The equation I wish to use to describe the chick growth is the logistic
growth function:
W = Asym/ (1+exp(-K*(x-xmid))
Where Asym = Asymptotic size, K is the rate constant of the equation and xmid
is the inflection point
This is very similar to the forumla used by the self-start function SSlogis,
the Asym and xmid
parameters are the same, however parameter K in the above equation is not the
same as parameter scal in
the simple logistic equation used by R.
My initial solution to this was to create my own self-start function using
the logistInit
function in Pinheiro & Bates (2000) as a guide. Thus, my function went as
follows:
logisticSS<-
function(mCall, LHS, data)
{
xy <- sortedXyData(mCall[["x"]], LHS, data)
if (nrow(xy) <3){
stop ("Too few distinct input values to fit a logistic")
}
Asym <- max(abs(xy[,"y"]))
xmid<-NLSstClosestX(xy, 0.5*Asym)
K <- (0.781/(NLSstClosestX(xy, 0.8*Asym)-NLSstClosestX(xy, 0.15*Asym)))*4
value<-c(Asym, xmid, K)
names(value)<- mCall[c("Asym", "xmid", "K")]
value
}
I then constructed a selfStart model by:
logistic<-
selfStart ( ~ Asym/ (1 + exp ( - K * ( x - xmid))), initial = logisticSS,
parameters = c ("Asym", "xmid", "K"))
for which the getInitial function is able to provide initial values.
Problems arise when I try and use this function in an nlsList call, where I get
a warning saying: number of iterations exceeds 50, and when I use it with the
nlme, which can spend hours computing. only to result in failure to converge.
My question is how would I go about improving my self start function to
allieviate these problems. I am sure there is a way to do it but I have been
unsuccessful in finding it so far. Any help would be much appreciated.
Thank you,
Ian Roberts
More information about the R-help
mailing list