[R] deriv when one term is indexed

Ken Knoblauch knoblauch at lyon.inserm.fr
Sat Nov 18 09:42:56 CET 2006


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