[R] deriv when one term is indexed
Gabor Grothendieck
ggrothendieck at gmail.com
Sat Nov 18 02:46:37 CET 2006
Please provide reproducible code which shows the error.
> dd <- structure(list(Lum = c(0.15, 0.07, 0.1, 0.19, 0.4, 0.73, 1.2,
+ 1.85, 2.91, 3.74, 5.08, 6.43, 8.06, 9.84, 12, 14.2, 16.6, 0.1,
+ 0.1, 0.17, 0.46, 1.08, 2.22, 3.74, 5.79, 8.36, 11.6, 15.4, 19.9,
+ 24.6, 30.4, 36.1, 43, 49.9, 0.06, 0.06, 0.08, 0.13, 0.25, 0.43,
+ 0.66, 1.02, 1.46, 1.93, 2.49, 3.2, 3.96, 4.9, 5.68, 6.71, 7.93
+ ), GL = as.integer(c(0, 15, 31, 47, 63, 79, 95, 111, 127, 143,
+ 159, 175, 191, 207, 223, 239, 255, 0, 15, 31, 47, 63, 79, 95,
+ 111, 127, 143, 159, 175, 191, 207, 223, 239, 255, 0, 15, 31,
+ 47, 63, 79, 95, 111, 127, 143, 159, 175, 191, 207, 223, 239,
+ 255)), Gun = structure(as.integer(c(2, 2, 2, 2, 2, 2, 2, 2, 2,
+ 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
+ 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
+ )), .Label = c("bleu", "rouge", "vert"), class = "factor")), .Names
= c("Lum",
+ "GL", "Gun"), class = "data.frame", row.names = c("1", "2", "3",
+ "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15",
+ "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26",
+ "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37",
+ "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48",
+ "49", "50", "51"))
>
> 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).
>
> 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/
>
> ______________________________________________
> 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.
>
More information about the R-help
mailing list