[R] deriv when one term is indexed

[Ken Knoblauch]

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

Here are the data, if of interest:
dd
     Lum  GL   Gun
1   0.15   0 rouge
2   0.07  15 rouge
3   0.10  31 rouge
4   0.19  47 rouge
5   0.40  63 rouge
6   0.73  79 rouge
7   1.20  95 rouge
8   1.85 111 rouge
9   2.91 127 rouge
10  3.74 143 rouge
11  5.08 159 rouge
12  6.43 175 rouge
13  8.06 191 rouge
14  9.84 207 rouge
15 12.00 223 rouge
16 14.20 239 rouge
17 16.60 255 rouge
18  0.10   0  vert
19  0.10  15  vert
20  0.17  31  vert
21  0.46  47  vert
22  1.08  63  vert
23  2.22  79  vert
24  3.74  95  vert
25  5.79 111  vert
26  8.36 127  vert
27 11.60 143  vert
28 15.40 159  vert
29 19.90 175  vert
30 24.60 191  vert
31 30.40 207  vert
32 36.10 223  vert
33 43.00 239  vert
34 49.90 255  vert
35  0.06   0  bleu
36  0.06  15  bleu
37  0.08  31  bleu
38  0.13  47  bleu
39  0.25  63  bleu
40  0.43  79  bleu
41  0.66  95  bleu
42  1.02 111  bleu
43  1.46 127  bleu
44  1.93 143  bleu
45  2.49 159  bleu
46  3.20 175  bleu
47  3.96 191  bleu
48  4.90 207  bleu
49  5.68 223  bleu
50  6.71 239  bleu
51  7.93 255  bleu

and here is the sessionInfo
R version 2.4.0 Patched (2006-11-10 r39843)
i386-apple-darwin8.8.1

locale:
C

attached base packages:
[1] "stats"     "graphics"  "grDevices" "utils"     "datasets"
[6] "methods"   "base"

other attached packages:
     boot      MASS   lattice
 "1.2-26"  "7.2-29" "0.14-13"

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/