[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