[R] Adjust Richards model with nlme

[David Winsemius]

> On Jan 19, 2016, at 3:20 AM, raphael fernandes <raphaelfsa89 at yahoo.com.br> wrote:
> 
> Hello everyone,
> 
> I have tried to adjust the non linear Richards growth model with this script:
> 
> richards <- function(x,beta1,beta2,beta3,beta4)
>     beta1*(1-beta2*exp(-x*beta3))^beta4
> richards <- deriv(~beta1*(1-beta2*exp(-x*beta3))^beta4,c("beta1","beta2","beta3","beta4"),function(x,beta1,beta2,beta3,beta4){})
> FF.richards <- nls(Peso~richards(Idade,beta1,beta2,beta3,beta4),data=femeas,start=c(beta1=370,beta2=9.7,beta3=272,beta4=0.00003))
> 
> but I receive the follow error:
> 
> Erro em qr.default(.swts * attr(rhs, "gradient")) : 
> NA/NaN/Inf em chamada de função externa (argumento 1)
> Além disso: Mensagens de aviso perdidas:
> In log(.expr5) : NaNs produzidos

Yes. Any interim trials for that iterative procedure that had a negative value for expr5 in the `richards` object would trigger that error:

richards
function (x, beta1, beta2, beta3, beta4) 
{
    .expr3 <- exp(-x * beta3)
    .expr5 <- 1 - beta2 * .expr3
    .expr6 <- .expr5^beta4
    .expr9 <- .expr5^(beta4 - 1)
    .value <- beta1 * .expr6
    .grad <- array(0, c(length(.value), 4L), list(NULL, c("beta1", 
        "beta2", "beta3", "beta4")))
    .grad[, "beta1"] <- .expr6
    .grad[, "beta2"] <- -(beta1 * (.expr9 * (beta4 * .expr3)))
    .grad[, "beta3"] <- beta1 * (.expr9 * (beta4 * (beta2 * (.expr3 * 
        x))))
    .grad[, "beta4"] <- beta1 * (.expr6 * log(.expr5))
    attr(.value, "gradient") <- .grad
    .value
}

> 

> What can I do to fix it?

You might try posting some data and seeing if the optimization gurus have insights. I hesitate to offer the hack I constructed to sidestep the log(neg_num) problem since I have no way to testing it for sensibility , and I'm certainly not in the guru-category.

-- 
David
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.

David Winsemius
Alameda, CA, USA