[R] deriv when one term is indexed

Ken Knoblauch knoblauch at lyon.inserm.fr
Sat Nov 18 09:46:12 CET 2006


Sorry, a slight glitch got in because I subsetted the original data frame.
The last part, that works should read:

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]))


#####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



Ken Knoblauch 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/
>


-- 
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