[R] deriv when one term is indexed

Spencer Graves spencer.graves at pdf.com
Sun Nov 19 22:36:42 CET 2006


      You ask "whether there would be a way to adapt the deriv and 
deriv3 functions to deal with formulas that have an indexed term", 
effectively allowing differentiation with respect to a vector.  As Simon 
Blomberg famously said, "This is R. There is no if. Only how."(*)

      The implementation, however, might not be easy.  Let's start by 
looking at 'deriv': 

 > deriv
function (expr, ...)
UseMethod("deriv")
<environment: namespace:stats>

      To get beyond this, try the following: 

 > methods("deriv")
[1] deriv.default deriv.formula

      Next, look at 'deriv.default': 

 deriv.default
function (expr, namevec, function.arg = NULL, tag = ".expr",
    hessian = FALSE, ...)
.Internal(deriv.default(expr, namevec, function.arg, tag, hessian))
<environment: namespace:stats>

      The call to '.Internal' says that 'deriv.default' is written in 
compiled code of some kind.  If we similarly check 'deriv.formula', and 
'deriv3', we get the same results. 

      Someone else might offer a more optimistic perspective, but this 
seems like a project beyond my available time at the moment. 

      If you told us more about your motivating example, someone might 
be able to help with a special case. 

      Sorry I couldn't be more encouraging. 
      Spencer Graves
(*) Reference: 
library(fortunes)
Ftns <- read.fortunes()
sel <- regexpr("This is R", Ftns$quote)
which(sel>0)
[1] 109
 > fortunes(109)

Ken Knoblauch wrote:
> I don't know why I get a different result than you do, though I'll
> play around with it some and check the code to see if I can figure
> out what is going on, on my end.
>
> I will add that by transforming GL so that it spans the interval [0, 1]
> instead
> of [0, 255], the correlations among the parameters decreased, also
> evidenced in
> the plot and the pairs plots of the profile object, but this didn't change
> whether these functions worked for me on the model when i did not add
> derivative information.
>
> But, this is a bit of a diversion, since we haven't yet addressed the
> question
> to which my Subject heading refers and remains the question for which
> I'm most seeking an answer, direction or guidance, and that is
> whether there would be a way to adapt the deriv and deriv3 functions
> to deal with formulas that have an indexed term as in the one
> in my situation, i.e.,
>
> Lum ~ Blev + beta[Gun] * GL^gamm
>
> Any ideas on that?  Thank you.
>
> ken
>
>
>
>
> Gabor Grothendieck wrote:
>   
>> I mixed up examples.  Here it is again.  As with the last one
>> confint(dd.plin) gives an error (which I assume is a problem with
>> confint that needs to be fixed) but other than that it works without
>> issuing errors and I assume you don't need the confint(dd.plin)
>> in any case since dd.plin is just being used to get starting values.
>>
>>     
>>> gg <- model.matrix(~ Gun/GL - Gun, dd)
>>> dd.plin <- nls(Lum ~ gg^gamm, dd, start = list(gamm = 2.4),
>>>       
>> +                                alg = "plinear"
>> +                                )
>>     
>>> confint(dd.plin)
>>>       
>> Waiting for profiling to be done...
>> Error in eval(expr, envir, enclos) : (subscript) logical subscript too
>> long
>>     
>>> 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.nls)
>>>       
>> Waiting for profiling to be done...
>>                2.5%        97.5%
>> Blev  -1.612492e-01 2.988387e-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
>>     
>>> 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), data = dd,
>>>       
>> +                        start = list(Blev = B[2], beta = B[3:5], gamm =
>> B[1]))
>>     
>>> confint(dd.nls.d2)
>>>       
>> Waiting for profiling to be done...
>>                2.5%        97.5%
>> Blev  -1.612492e-01 2.988391e-02
>> beta1  1.269000e-05 1.792914e-05
>> beta2  3.844041e-05 5.388546e-05
>> beta3  6.108281e-06 8.762678e-06
>> gamm   2.481102e+00 2.542966e+00
>>     
>>> R.version.string # XP
>>>       
>> [1] "R version 2.4.0 Patched (2006-10-24 r39722)"
>>     
>> On 11/18/06, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:
>>     
>>> This works for me in terms of giving results without error messages
>>> except for the confint(dd.plin) which I assume you don't really need
>>> anyways.
>>>
>>>       
>>>> gg <- model.matrix(~ Gun/GL - Gun, dd)
>>>> dd.plin <- nls(Lum ~ gg^gamm, dd, start = list(gamm = 2.4),
>>>>         
>>> +                                alg = "plinear"
>>> +                                )
>>>       
>>>> confint(dd.plin)
>>>>         
>>> Waiting for profiling to be done...
>>> Error in eval(expr, envir, enclos) : (subscript) logical subscript too
>>> long
>>>       
>>>> B <- as.vector(coef(dd.nls0))
>>>> 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.nls)
>>>>         
>>> Waiting for profiling to be done...
>>>               2.5%        97.5%
>>> Blev  -1.612492e-01 2.988386e-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
>>>       
>>>> 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), data = dd,
>>>>         
>>> +                        start = list(Blev = B[2], beta = B[3:5], gamm =
>>> B[1]))
>>>       
>>>> confint(dd.nls.d2)
>>>>         
>>> Waiting for profiling to be done...
>>>               2.5%        97.5%
>>> Blev  -1.612492e-01 2.988391e-02
>>> beta1  1.269000e-05 1.792914e-05
>>> beta2  3.844041e-05 5.388546e-05
>>> beta3  6.108281e-06 8.762678e-06
>>> gamm   2.481102e+00 2.542966e+00
>>>       
>>>> R.version.string # XP
>>>>         
>>> [1] "R version 2.4.0 Patched (2006-10-24 r39722)"
>>>       
>>> On 11/18/06, Ken Knoblauch <knoblauch at lyon.inserm.fr> 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/
>>>>
>>>>
>>>>         
>
>
>



More information about the R-help mailing list