[R-sig-ME] nlmer model definition and starting values

[Sebastian P. Luque]

Sorry, I forgot to include the self-start function:


---<--------------------cut here---------------start------------------->---
SSGompertz <- selfStart(~ Linf * exp(-b * exp(-k * x)),
                        function(mCall, data, LHS) {
                            xy <- sortedXyData(mCall[["x"]], LHS, data)
                            if (nrow(xy) < 4) {
                                stop("too few distinct input values to fit Gompertz model")
                            }
                            xyL <- xy
                            xyL$y <- log(abs(xyL$y))
                            pars <- NLSstAsymptotic(xyL)
                            pars <- coef(nls(y ~ exp(-b2 * b3 ^ x), data=xy, algorithm="plinear",
                                             start=c(b2=pars[["b1"]], b3=exp(-exp(pars[["lrc"]])))))
                            val <- pars[c(3, 1, 2)]
                            val[3] <- -log(val[3])
                            names(val) <- mCall[c("Linf", "b", "k")]
                            val
                        }, c("Linf", "b", "k"))
---<--------------------cut here---------------end--------------------->---


Thanks.



On Wed, 29 Apr 2009 00:05:08 -0500,
"Sebastian P. Luque" <spluque at gmail.com> wrote:

> Hi, I can't understand how to specify starting values for the
> following nonlinear growth model.  The fake data produced below is
> designed to closely match a real case scenario.  The model is for the
> following Gompertz growth parameterization:

> y=Linf * exp(-b * exp(-k * x))

> where y is some size variable, and Linf, b, and k are parameters
> indicating asymptotic size, size at birth, and a growth constant,
> respectively.  I've hacked a self-start function for it.  Each row in
> the fake data represents a single individual, with each individual
> measured only once.  The data correspond to 4 species and 2 sexes,
> with each combination having its own set of parameters, and some
> normal heteroscedastic error added.

> I would like to start with the assumption that b is a random effect,
> while species and sex have fixed effects on Linf and k.


> "gompertz.fun" <- function(x, Linf, b, k) {Linf * exp(-b * exp(-k *
> x))}

> ## Set parameters for, say, 4 species and 2 sexes

> grps <- expand.grid(species=LETTERS[1:4], sex=c("M", "F")) parms <-
> c("Linf", "b", "k") parms.mat <- matrix(c(432.3, 411.5, 389.7, 448.2,
> 381.5, 360.5, 338.1, 385.7, 1.16, 0.59, 0.73, 0.61, 7.56, 0.49, 0.51,
> 0.51, 0.15, 0.08, 0.07, 0.09, 0.29, 0.14, 0.11, 0.14), ncol=3,
> dimnames=list(paste(grps[, 1], grps[, 2]), parms))

> ## Generate random x for each combination of species and sex, and y ##
> based on the parameters, with different y normal random noise ##
> (heteroscedastic) per combination set.seed(123) x <- rep(seq(1, 60,
> 3), 8) + runif(160, min=0, max=2) fake <-
> data.frame(grps[rep(seq(nrow(grps)), each=20), ], x) fake <-
> with(fake, fake[order(species, sex, x), ]) y <- apply(fake, 1,
> function(k) { grp <- paste(k[1], k[2]) parms <-
> parms.mat[rownames(parms.mat) %in% grp, ]
> gompertz.fun(as.numeric(k[3]), parms[1], parms[2], parms[3]) }) ## Add
> the increasing error fake$y <- rnorm(160, mean=y, sd=seq(5, 30,
> length.out=20)) ## Quick visualization xyplot(y ~ x | species,
> data=fake, groups=sex, aspect=0.7, scales=list(alternating=1,
> tck=c(0.75, 0), rot=c(0, 90)))


> How could I write the nlmer call for such a model, and how could I
> anticipate the order of starting values to supply.  Thanks in advance
> for any feedback on this.



-- 
Seb