algal nonlinear mixed effects problem

Thanks, Ben, for the illustration of how to use derive(), deparse(),
and eval() for setting up such functions. The following is based on
lme4_0.999999-0.

I modified Ben's attempt at a fit in two ways to get convergence.
(1) The start argument threw an Error; I provide the values in a
vector rather than a list. This might be related to differences
between nlmer implementations.
(2) The start value for asymp.R2 had to be moved away from 0.

The result suggests that there is no reliable evidence for the
variance component associated with "asymp.R1"

A part of the protocol is listed below.
Happy New Year to everyone.

Reinhold Kliegl

Here is the protocol:

# 4. Attempt the fit
nlmerfit2 <- nlmer(

+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  list(nlpars=c(asymp.L=0.7,
+      asymp.R1=0.6,asymp.R2=0,asymp.R3=0,xmid=5,scale=1)),data=d)
Fehler: length(start$fixef) > 0 is not TRUE

nlmerfit2a <- nlmer(

+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  c(asymp.L=0.6,
asymp.R1=0,asymp.R2=0,asymp.R3=0,xmid=5,scale=1),data=d)
Fehler in mer_finalize(ans) : Downdated X'X is not positive definite, 5.

nlmerfit2b <- nlmer(

+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  c(asymp.L=0.6,
asymp.R1=.1,asymp.R2=0,asymp.R3=0,xmid=5,scale=1),data=d)

nlmerfit2b

Nonlinear mixed model fit by the Laplace approximation
Formula: X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid,
scale) ~      asymp.R1 | Individual
   Data: d
   AIC   BIC   logLik deviance
 16.11 32.02 -0.05555   0.1111
Random effects:
 Groups     Name     Variance   Std.Dev.
 Individual asymp.R1 8.9101e-22 2.985e-11
 Residual            2.0575e-03 4.536e-02
Number of obs: 54, groups: Individual, 9

Fixed effects:
           Estimate Std. Error t value
asymp.L  -2.173e+02  2.187e+05  -0.001
asymp.R1  6.140e+01  5.002e+04   0.001
asymp.R2 -1.917e-02  4.370e+00  -0.004
asymp.R3 -1.207e-01  2.752e+01  -0.004
xmid      6.152e+03  6.167e+06   0.001
scale     4.807e+03  3.610e+06   0.001

Correlation of Fixed Effects:
         asym.L asy.R1 asy.R2 asy.R3 xmid
asymp.R1 -0.357
asymp.R2 -0.677 -0.446
asymp.R3 -0.677 -0.446  1.000
xmid     -1.000  0.357  0.677  0.677
scale    -0.597  0.962 -0.185 -0.185  0.597


On Mon, Dec 31, 2012 at 10:09 PM, Ben Bolker <bbolker at gmail.com> wrote:

Ben Bolker <bbolker at ...> writes:

  Inspired by a recent question of Diego Pujoni's on the list
( http://article.gmane.org/gmane.comp.lang.r.lme4.devel/9363/match= ),
I started working on the problem, with two realizations:

 1. the problem really needs a nonlinear fit: see
http://rpubs.com/bbolker/3363

 2.  I have no idea how to do an nlmer fit with a non-trivial
fixed-effects model (e.g. in nlme we could do

nlmefit1 <- nlme(model = X ~ SSfpl(Day, asymp.L, asymp.R,
                                   xmid, scale),
   fixed = list(asymp.R ~ Group, xmid + scale + asymp.L ~ 1),
   random = asymp.R ~ 1 | Individual, ...)

  This was raised earlier on stack overflow:

http://stackoverflow.com/questions/11056625/how-to-add-fixed-effect-to-four-parameter-logistic-model-in-nlmer

  Reinhold Kliegl e-mailed off-list to say that he was able to do
it / had to do it by hacking up an objective function that accounts
for the grouping (i.e., more or less doing it by hand).  I have
implemented this and updated http://rpubs.com/bbolker/3363 , which
also shows a solution with AD Model Builder; it turns out that at
least the development version of nlmer is a little too fragile to
fit the model I wanted to fit (I haven't tested with the stable
version).

Thanks, Ben, for the illustration of how to use derive(), deparse(),
and eval() for setting up such functions. The following is based on
lme4_0.999999-0.

I modified Ben's attempt at a fit in two ways to get convergence.
(1) The start argument threw an Error; I provide the values in a
vector rather than a list. This might be related to differences
between nlmer implementations.
(2) The start value for asymp.R2 had to be moved away from 0.

The result suggests that there is no reliable evidence for the
variance component associated with "Individual asymp.R1"

A part of the protocol is listed below.
Happy New Year to everyone.

Reinhold Kliegl

Here is the protocol:

# 4. Attempt the fit
nlmerfit2 <- nlmer(

+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  list(nlpars=c(asymp.L=0.7,
+      asymp.R1=0.6,asymp.R2=0,asymp.R3=0,xmid=5,scale=1)),data=d)
Fehler: length(start$fixef) > 0 is not TRUE

# ERROR: length(start$fixef) > 0 is not TRUE

nlmerfit2a <- nlmer(

+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  c(asymp.L=0.6,
asymp.R1=0,asymp.R2=0,asymp.R3=0,xmid=5,scale=1),data=d)
Fehler in mer_finalize(ans) : Downdated X'X is not positive definite, 5.

# Fehler in mer_finalize(ans) : Downdated X'X is not positive definite, 5.

nlmerfit2b <- nlmer(

+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  c(asymp.L=0.6,
asymp.R1=.1,asymp.R2=0,asymp.R3=0,xmid=5,scale=1),data=d)

nlmerfit2b

Nonlinear mixed model fit by the Laplace approximation
Formula: X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid,
scale) ~      asymp.R1 | Individual
   Data: d
   AIC   BIC   logLik deviance
 16.11 32.02 -0.05555   0.1111
Random effects:
 Groups     Name     Variance   Std.Dev.
 Individual asymp.R1 8.9101e-22 2.985e-11
 Residual            2.0575e-03 4.536e-02
Number of obs: 54, groups: Individual, 9

Fixed effects:
           Estimate Std. Error t value
asymp.L  -2.173e+02  2.187e+05  -0.001
asymp.R1  6.140e+01  5.002e+04   0.001
asymp.R2 -1.917e-02  4.370e+00  -0.004
asymp.R3 -1.207e-01  2.752e+01  -0.004
xmid      6.152e+03  6.167e+06   0.001
scale     4.807e+03  3.610e+06   0.001

Correlation of Fixed Effects:
         asym.L asy.R1 asy.R2 asy.R3 xmid
asymp.R1 -0.357
asymp.R2 -0.677 -0.446
asymp.R3 -0.677 -0.446  1.000
xmid     -1.000  0.357  0.677  0.677
scale    -0.597  0.962 -0.185 -0.185  0.597


On Mon, Dec 31, 2012 at 10:09 PM, Ben Bolker <bbolker at gmail.com> wrote:

Ben Bolker <bbolker at ...> writes:

  Inspired by a recent question of Diego Pujoni's on the list
( http://article.gmane.org/gmane.comp.lang.r.lme4.devel/9363/match= ),
I started working on the problem, with two realizations:

 1. the problem really needs a nonlinear fit: see
http://rpubs.com/bbolker/3363

 2.  I have no idea how to do an nlmer fit with a non-trivial
fixed-effects model (e.g. in nlme we could do

nlmefit1 <- nlme(model = X ~ SSfpl(Day, asymp.L, asymp.R,
                                   xmid, scale),
   fixed = list(asymp.R ~ Group, xmid + scale + asymp.L ~ 1),
   random = asymp.R ~ 1 | Individual, ...)

  This was raised earlier on stack overflow:

http://stackoverflow.com/questions/11056625/how-to-add-fixed-effect-to-four-parameter-logistic-model-in-nlmer

  Reinhold Kliegl e-mailed off-list to say that he was able to do
it / had to do it by hacking up an objective function that accounts
for the grouping (i.e., more or less doing it by hand).  I have
implemented this and updated http://rpubs.com/bbolker/3363 , which
also shows a solution with AD Model Builder; it turns out that at
least the development version of nlmer is a little too fragile to
fit the model I wanted to fit (I haven't tested with the stable
version).

On 13-01-01 05:17 AM, Reinhold Kliegl wrote:

Thanks, Ben, for the illustration of how to use derive(), deparse(),
and eval() for setting up such functions. The following is based on
lme4_0.999999-0.

I modified Ben's attempt at a fit in two ways to get convergence.
(1) The start argument threw an Error; I provide the values in a
vector rather than a list. This might be related to differences
between nlmer implementations.
(2) The start value for asymp.R2 had to be moved away from 0.

The result suggests that there is no reliable evidence for the
variance component associated with "Individual asymp.R1"

  The results are crazy, though -- take a look at the fixed effect
coefficients, or the predicted values ...

A part of the protocol is listed below.
Happy New Year to everyone.

Reinhold Kliegl

Here is the protocol:

# 4. Attempt the fit
nlmerfit2 <- nlmer(

+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  list(nlpars=c(asymp.L=0.7,
+      asymp.R1=0.6,asymp.R2=0,asymp.R3=0,xmid=5,scale=1)),data=d)
Fehler: length(start$fixef) > 0 is not TRUE

# ERROR: length(start$fixef) > 0 is not TRUE

nlmerfit2a <- nlmer(

+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  c(asymp.L=0.6,
asymp.R1=0,asymp.R2=0,asymp.R3=0,xmid=5,scale=1),data=d)
Fehler in mer_finalize(ans) : Downdated X'X is not positive definite, 5.

# Fehler in mer_finalize(ans) : Downdated X'X is not positive definite, 5.

nlmerfit2b <- nlmer(

+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  c(asymp.L=0.6,
asymp.R1=.1,asymp.R2=0,asymp.R3=0,xmid=5,scale=1),data=d)

nlmerfit2b

Nonlinear mixed model fit by the Laplace approximation
Formula: X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid,
scale) ~      asymp.R1 | Individual
   Data: d
   AIC   BIC   logLik deviance
 16.11 32.02 -0.05555   0.1111
Random effects:
 Groups     Name     Variance   Std.Dev.
 Individual asymp.R1 8.9101e-22 2.985e-11
 Residual            2.0575e-03 4.536e-02
Number of obs: 54, groups: Individual, 9

Fixed effects:
           Estimate Std. Error t value
asymp.L  -2.173e+02  2.187e+05  -0.001
asymp.R1  6.140e+01  5.002e+04   0.001
asymp.R2 -1.917e-02  4.370e+00  -0.004
asymp.R3 -1.207e-01  2.752e+01  -0.004
xmid      6.152e+03  6.167e+06   0.001
scale     4.807e+03  3.610e+06   0.001

Correlation of Fixed Effects:
         asym.L asy.R1 asy.R2 asy.R3 xmid
asymp.R1 -0.357
asymp.R2 -0.677 -0.446
asymp.R3 -0.677 -0.446  1.000
xmid     -1.000  0.357  0.677  0.677
scale    -0.597  0.962 -0.185 -0.185  0.597


On Mon, Dec 31, 2012 at 10:09 PM, Ben Bolker <bbolker at gmail.com> wrote:

Ben Bolker <bbolker at ...> writes:

  Inspired by a recent question of Diego Pujoni's on the list
( http://article.gmane.org/gmane.comp.lang.r.lme4.devel/9363/match= ),
I started working on the problem, with two realizations:

 1. the problem really needs a nonlinear fit: see
http://rpubs.com/bbolker/3363

 2.  I have no idea how to do an nlmer fit with a non-trivial
fixed-effects model (e.g. in nlme we could do

nlmefit1 <- nlme(model = X ~ SSfpl(Day, asymp.L, asymp.R,
                                   xmid, scale),
   fixed = list(asymp.R ~ Group, xmid + scale + asymp.L ~ 1),
   random = asymp.R ~ 1 | Individual, ...)

  This was raised earlier on stack overflow:

http://stackoverflow.com/questions/11056625/how-to-add-fixed-effect-to-four-parameter-logistic-model-in-nlmer

  Reinhold Kliegl e-mailed off-list to say that he was able to do
it / had to do it by hacking up an objective function that accounts
for the grouping (i.e., more or less doing it by hand).  I have
implemented this and updated http://rpubs.com/bbolker/3363 , which
also shows a solution with AD Model Builder; it turns out that at
least the development version of nlmer is a little too fragile to
fit the model I wanted to fit (I haven't tested with the stable
version).