[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