[R] deriv when one term is indexed

Gabor Grothendieck ggrothendieck at gmail.com
Sat Nov 18 15:36:05 CET 2006


I mixed up examples.  Here it is again.  As with the last one
confint(dd.plin) gives an error (which I assume is a problem with
confint that needs to be fixed) but other than that it works without
issuing errors and I assume you don't need the confint(dd.plin)
in any case since dd.plin is just being used to get starting values.

> gg <- model.matrix(~ Gun/GL - Gun, dd)
> dd.plin <- nls(Lum ~ gg^gamm, dd, start = list(gamm = 2.4),
+                                alg = "plinear"
+                                )
> confint(dd.plin)
Waiting for profiling to be done...
Error in eval(expr, envir, enclos) : (subscript) logical subscript too long
>
>
> B <- as.vector(coef(dd.plin))
> st <-list(Blev = B[2], beta = B[3:5], gamm = B[1])
>
>
>
> dd.nls <- nls(Lum ~ Blev + beta[Gun] * GL^gamm,
+                        data = dd, start = st
+                        )
>
> confint(dd.nls)
Waiting for profiling to be done...
               2.5%        97.5%
Blev  -1.612492e-01 2.988387e-02
beta1  6.108282e-06 8.762679e-06
beta2  1.269000e-05 1.792914e-05
beta3  3.844042e-05 5.388546e-05
gamm   2.481102e+00 2.542966e+00
>
> dd.deriv2 <- function (Blev, beta, gamm, GL)
+ {
+    .expr1 <- GL^gamm
+    .value <- Blev + rep(beta, each = 17) * .expr1
+    .grad <- array(0, c(length(.value), 5), list(NULL, c("Blev",
+        "beta.rouge", "beta.vert", "beta.bleu", "gamm")))
+    .grad[, "Blev"] <- 1
+    .grad[1:17, "beta.rouge"] <- .expr1[1:17]
+    .grad[18:34, "beta.vert"] <- .expr1[1:17]
+    .grad[35:51, "beta.bleu"] <- .expr1[1:17]
+    .grad[, "gamm"] <- ifelse(GL == 0, 0, rep(beta, each = 17) * (.expr1 *
+ log(GL)))
+    attr(.value, "gradient") <- .grad
+    .value
+ }
>
>
> dd.nls.d2 <- nls(Lum ~ dd.deriv2(Blev, beta, gamm, GL), data = dd,
+                        start = list(Blev = B[2], beta = B[3:5], gamm = B[1]))
>
> confint(dd.nls.d2)
Waiting for profiling to be done...
               2.5%        97.5%
Blev  -1.612492e-01 2.988391e-02
beta1  1.269000e-05 1.792914e-05
beta2  3.844041e-05 5.388546e-05
beta3  6.108281e-06 8.762678e-06
gamm   2.481102e+00 2.542966e+00
>
> R.version.string # XP
[1] "R version 2.4.0 Patched (2006-10-24 r39722)"
>


On 11/18/06, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:
> This works for me in terms of giving results without error messages
> except for the confint(dd.plin) which I assume you don't really need anyways.
>
> > gg <- model.matrix(~ Gun/GL - Gun, dd)
> > dd.plin <- nls(Lum ~ gg^gamm, dd, start = list(gamm = 2.4),
> +                                alg = "plinear"
> +                                )
> > confint(dd.plin)
> Waiting for profiling to be done...
> Error in eval(expr, envir, enclos) : (subscript) logical subscript too long
> >
> >
> > B <- as.vector(coef(dd.nls0))
> > st <-list(Blev = B[2], beta = B[3:5], gamm = B[1])
> >
> >
> >
> > dd.nls <- nls(Lum ~ Blev + beta[Gun] * GL^gamm,
> +                        data = dd, start = st
> +                        )
> >
> > confint(dd.nls)
> Waiting for profiling to be done...
>               2.5%        97.5%
> Blev  -1.612492e-01 2.988386e-02
> beta1  6.108282e-06 8.762679e-06
> beta2  1.269000e-05 1.792914e-05
> beta3  3.844042e-05 5.388546e-05
> gamm   2.481102e+00 2.542966e+00
> >
> > dd.deriv2 <- function (Blev, beta, gamm, GL)
> + {
> +    .expr1 <- GL^gamm
> +    .value <- Blev + rep(beta, each = 17) * .expr1
> +    .grad <- array(0, c(length(.value), 5), list(NULL, c("Blev",
> +        "beta.rouge", "beta.vert", "beta.bleu", "gamm")))
> +    .grad[, "Blev"] <- 1
> +    .grad[1:17, "beta.rouge"] <- .expr1[1:17]
> +    .grad[18:34, "beta.vert"] <- .expr1[1:17]
> +    .grad[35:51, "beta.bleu"] <- .expr1[1:17]
> +    .grad[, "gamm"] <- ifelse(GL == 0, 0, rep(beta, each = 17) * (.expr1 *
> + log(GL)))
> +    attr(.value, "gradient") <- .grad
> +    .value
> + }
> >
> >
> > dd.nls.d2 <- nls(Lum ~ dd.deriv2(Blev, beta, gamm, GL), data = dd,
> +                        start = list(Blev = B[2], beta = B[3:5], gamm = B[1]))
> >
> > confint(dd.nls.d2)
> Waiting for profiling to be done...
>               2.5%        97.5%
> Blev  -1.612492e-01 2.988391e-02
> beta1  1.269000e-05 1.792914e-05
> beta2  3.844041e-05 5.388546e-05
> beta3  6.108281e-06 8.762678e-06
> gamm   2.481102e+00 2.542966e+00
> >
> > R.version.string # XP
> [1] "R version 2.4.0 Patched (2006-10-24 r39722)"
> >
>
>
> On 11/18/06, Ken Knoblauch <knoblauch at lyon.inserm.fr> wrote:
> > Thank you for your rapid response.
> >
> > This is reproducible on my system.  Here it is again, with, I hope,
> > sufficient detail to properly document what does not work and what
> > does on my system,
> >
> > But my original question, properly motivated or not, concerns whether
> > there is a way to use or adapt deriv() to work if a term is indexed,
> > as here my term beta is indexed by Gun?
> >
> > In any case, it is puzzling that the error is not reproducible and so
> > I would be curious to track that down, if it is specific to my system.
> >
> > Thank you. Before retrying, I upgraded to the latest patched version
> > (details at the end).
> >
> > ####The data, again
> > dd
> > Lum     GL      Gun     mBl
> > 0.15    0       rouge   0.09
> > 0.07    15      rouge   0.01
> > 0.1     31      rouge   0.04
> > 0.19    47      rouge   0.13
> > 0.4     63      rouge   0.34
> > 0.73    79      rouge   0.67
> > 1.2     95      rouge   1.14
> > 1.85    111     rouge   1.79
> > 2.91    127     rouge   2.85
> > 3.74    143     rouge   3.68
> > 5.08    159     rouge   5.02
> > 6.43    175     rouge   6.37
> > 8.06    191     rouge   8
> > 9.84    207     rouge   9.78
> > 12      223     rouge   11.94
> > 14.2    239     rouge   14.14
> > 16.6    255     rouge   16.54
> > 0.1     0       vert    0.04
> > 0.1     15      vert    0.04
> > 0.17    31      vert    0.11
> > 0.46    47      vert    0.4
> > 1.08    63      vert    1.02
> > 2.22    79      vert    2.16
> > 3.74    95      vert    3.68
> > 5.79    111     vert    5.73
> > 8.36    127     vert    8.3
> > 11.6    143     vert    11.54
> > 15.4    159     vert    15.34
> > 19.9    175     vert    19.84
> > 24.6    191     vert    24.54
> > 30.4    207     vert    30.34
> > 36.1    223     vert    36.04
> > 43      239     vert    42.94
> > 49.9    255     vert    49.84
> > 0.06    0       bleu    0
> > 0.06    15      bleu    0
> > 0.08    31      bleu    0.02
> > 0.13    47      bleu    0.07
> > 0.25    63      bleu    0.19
> > 0.43    79      bleu    0.37
> > 0.66    95      bleu    0.6
> > 1.02    111     bleu    0.96
> > 1.46    127     bleu    1.4
> > 1.93    143     bleu    1.87
> > 2.49    159     bleu    2.43
> > 3.2     175     bleu    3.14
> > 3.96    191     bleu    3.9
> > 4.9     207     bleu    4.84
> > 5.68    223     bleu    5.62
> > 6.71    239     bleu    6.65
> > 7.93    255     bleu    7.87
> >
> > ###For initial values - this time using plinear algorithm insted of optim
> > gg <- model.matrix(~-1 + Gun/GL, dd)[ , c(4:6)]
> > dd.plin <- nls(Lum ~ cbind(rep(1, 51), gg^gamm), data = dd,
> >                                start = list(gamm = 2.4),
> >                                alg = "plinear"
> >                                )
> > B <- as.vector(coef(dd.plin))
> > st <-list(Blev = B[2], beta = B[3:5], gamm = B[1])
> >
> >
> > dd.nls <- nls(Lum ~ Blev + beta[Gun] * GL^gamm,
> >                        data = dd, start = st
> >                        )
> > confint(dd.plin)
> >
> > Waiting for profiling to be done...
> > Error in eval(expr, envir, enclos) : (subscript) logical subscript too long
> > ###  Here is the error that I observe
> > confint(dd.nls)
> > Waiting for profiling to be done...
> > Error in prof$getProfile() : step factor 0.000488281 reduced below
> > 'minFactor' of 0.000976562
> >
> > dd.deriv2 <- function (Blev, beta, gamm, GL)
> > {
> >    .expr1 <- GL^gamm
> >    .value <- Blev + rep(beta, each = 17) * .expr1
> >    .grad <- array(0, c(length(.value), 5), list(NULL, c("Blev",
> >        "beta.rouge", "beta.vert", "beta.bleu", "gamm")))
> >    .grad[, "Blev"] <- 1
> >    .grad[1:17, "beta.rouge"] <- .expr1[1:17]
> >    .grad[18:34, "beta.vert"] <- .expr1[1:17]
> >    .grad[35:51, "beta.bleu"] <- .expr1[1:17]
> >    .grad[, "gamm"] <- ifelse(GL == 0, 0, rep(beta, each = 17) * (.expr1 *
> > log(GL)))
> >    attr(.value, "gradient") <- .grad
> >    .value
> > }
> >
> > dd.nls.d2 <- nls(Lum ~ dd.deriv2(Blev, beta, gamm, GL), subset(dd, Gun !=
> > "gris"),
> >                        start = list(Blev = B[2], beta = B[3:5], gamm = B[1]))
> >
> >
> > #####But not here
> > confint(dd.nls.d2)
> > Waiting for profiling to be done...
> >               2.5%        97.5%
> > Blev  -1.612492e-01 2.988387e-02
> > beta1  1.269000e-05 1.792914e-05
> > beta2  3.844042e-05 5.388546e-05
> > beta3  6.108282e-06 8.762679e-06
> > gamm   2.481102e+00 2.542966e+00
> >
> > R version 2.4.0 Patched (2006-11-16 r39921)
> > i386-apple-darwin8.8.1
> >
> > locale:
> > C
> >
> > attached base packages:
> > [1] "stats"     "graphics"  "grDevices" "utils"     "datasets"
> > [6] "methods"   "base"
> >
> > other attached packages:
> >     MASS   lattice
> >  "7.2-29" "0.14-13"
> >
> > best,
> >
> > ken
> >
> >
> > Gabor Grothendieck wrote:
> > > Please provide reproducible code which shows the error.
> > >
> > >>
> > >> st <- list(Blev = -0.06551802, beta = c(1.509686e-05, 4.55525e-05,
> > > + 7.32272e-06), gamm = 2.51187)
> > >>
> > >> dd.nls <- nls(Lum ~ Blev + beta[Gun] * GL^gamm,
> > > +                        data = dd, start = st)
> > >>
> > >> confint(dd.nls)
> > > Waiting for profiling to be done...
> > >                2.5%        97.5%
> > > Blev  -1.612492e-01 2.988388e-02
> > > beta1  6.108282e-06 8.762679e-06
> > > beta2  1.269000e-05 1.792914e-05
> > > beta3  3.844042e-05 5.388546e-05
> > > gamm   2.481102e+00 2.542966e+00
> > >> R.version.string
> > > [1] "R version 2.4.0 Patched (2006-10-24 r39722)"
> > >
> > >
> > >
> > >
> > > On 11/17/06, Ken Knoblauch <knoblauch at lyon.inserm.fr> wrote:
> > >> Hi,
> > >>
> > >> I'm fitting a standard nonlinear model to the luminances measured
> > >> from the red, green and blue guns of a TV display, using nls.
> > >>
> > >> The call is:
> > >>
> > >> dd.nls <- nls(Lum ~ Blev + beta[Gun] * GL^gamm,
> > >>                        data = dd, start = st)
> > >> where st was initally estimated using optim()
> > >>
> > >> st
> > >> $Blev
> > >> [1] -0.06551802
> > >>
> > >> $beta
> > >> [1] 1.509686e-05 4.555250e-05 7.322720e-06
> > >>
> > >> $gamm
> > >> [1] 2.511870
> > >>
> > >> This works fine but I received an error message when I tried to
> > >> use confint().  I thought that getting derivatives with deriv might
> > >> help but found that deriv does not automatically handle the
> > >> indexing of the beta parameter.  I modified the output of deriv
> > >> from the case when the term beta is not indexed to give:
> > >>
> > >> dd.deriv2 <- function (Blev, beta, gamm, GL)
> > >> {
> > >>    .expr1 <- GL^gamm
> > >>    .value <- Blev + rep(beta, each = 17) * .expr1
> > >>    .grad <- array(0, c(length(.value), 5), list(NULL, c("Blev",
> > >>        "beta.rouge", "beta.vert", "beta.bleu", "gamm"  )))
> > >>    .grad[, "Blev"] <- 1
> > >>    .grad[1:17, "beta.rouge"] <- .expr1[1:17]
> > >>    .grad[18:34, "beta.vert"] <- .expr1[1:17]
> > >>    .grad[35:51, "beta.bleu"] <- .expr1[1:17]
> > >>    .grad[, "gamm"] <- ifelse(GL == 0, 0, rep(beta, each = 17) * (.expr1
> > >> *
> > >> log(GL)))
> > >>    attr(.value, "gradient") <- .grad
> > >>    .value
> > >> }
> > >>
> > >> which is not general but
> > >> now I can get a result from confint. My question:  Can deriv()
> > >> be made to handle an indexed term more automatically (elegantly)?
> > >> I think that this would become more urgent (or require more work
> > >> on my part) if I were to want the Hessian, too, for example, in
> > >> anticipation of using rms.curv as described in the on-line
> > >> complements for MASS.
> > >>
> > >> The plinear algorithm can be used for this model (it is similar in some
> > >> respects to the example given on p 218 of MASS, but the intercept
> > >> terms are indexed instead, there).
> > >>
> > >> Thanks for any suggestions.
> > >>
> > >> best,
> > >>
> > >> Ken
> > >>
> > >>
> > >> --
> > >> Ken Knoblauch
> > >> Inserm U371
> > >> Institut Cellule Souche et Cerveau
> > >> Département Neurosciences Intégratives
> > >> 18 avenue du Doyen Lépine
> > >> 69500 Bron
> > >> France
> > >> tel: +33 (0)4 72 91 34 77
> > >> fax: +33 (0)4 72 91 34 61
> > >> portable: +33 (0)6 84 10 64 10
> > >> http://www.lyon.inserm.fr/371/
> > >>
> > >> ______________________________________________
> > >> R-help at stat.math.ethz.ch mailing list
> > >> 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.
> > >>
> > >
> >
> >
> > --
> > Ken Knoblauch
> > Inserm U371
> > Institut Cellule Souche et Cerveau
> > Département Neurosciences Intégratives
> > 18 avenue du Doyen Lépine
> > 69500 Bron
> > France
> > tel: +33 (0)4 72 91 34 77
> > fax: +33 (0)4 72 91 34 61
> > portable: +33 (0)6 84 10 64 10
> > http://www.lyon.inserm.fr/371/
> >
> >
>



More information about the R-help mailing list