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

[ONKELINX, Thierry]

Dear Sebastian,

There is no need to define SSGompertz as SSgompertz is already included
in the stats package.

The helpfile for nlme() gives you a good example. Switching SSamymp into
SSgompertz is rather obvious. 

fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
            data = Loblolly,
            fixed = Asym + R0 + lrc ~ 1,
            random = Asym ~ 1,
            start = c(Asym = 103, R0 = -8.5, lrc = -3.3))
summary(fm1)
fm2 <- update(fm1, random = pdDiag(Asym + lrc ~ 1))
summary(fm2)
 
HTH,

Thierry

------------------------------------------------------------------------
----
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium 
tel. + 32 54/436 185
Thierry.Onkelinx at inbo.be 
www.inbo.be 

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to
say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of
data.
~ John Tukey

-----Oorspronkelijk bericht-----
Van: r-sig-mixed-models-bounces at r-project.org
[mailto:r-sig-mixed-models-bounces at r-project.org] Namens Sebastian P.
Luque
Verzonden: woensdag 29 april 2009 7:14
Aan: r-sig-mixed-models at r-project.org
Onderwerp: Re: [R-sig-ME] nlmer model definition and starting values

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

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer 
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in  this message 
and any annex are purely those of the writer and may not be regarded as stating 
an official position of INBO, as long as the message is not confirmed by a duly 
signed document.