4

In R, I would like a way to take symbolic derivatives of the right hand side of formulas which may include interaction terms, squared terms, etc.

For example, I would like to be able to take the derivative of the right hand side of each of the following two [edit:three] formulas with respect to x:

y~x+I(x^2)

y~x:z

EDIT: y~x*z

I would like a function which, when each of the above three formulas are input, returns 1+2x, z, and 1+z, respectively.

I've tried the following:

f1<-y~x+I(x^2)
deriv(f1,"x")
## Error in deriv.formula(f1, "x") : Function 'I' is not in the derivatives table

f2<-y~x:z
deriv(f2,"x")
## Error in deriv.formula(f2, "x") : Function '`:`' is not in the derivatives table

Is there any way to force R to recognize I(x^2) (or, similarly, I(x*z), etc.) as x^2 (respectively, x*z), x:z as x*z (in the mathematical sense), and x*z (in the formula sense) as x+z+x*z (in the mathematical sense) for purposes of calculating the derivative?

Second, is there a way to take the output from deriv() and reshape it to look like the right hand side of a formula? In particular, I know that D() will alleviate this issue and generate output in the form I desire (though D() can't handle a formula as input), but what if I want to take derivatives with respect to multiple variables? I can work around this by applying D() over and over for each variable I'd like to take the derivative with respect to, but it would be nice to simply input a character string of all such variables and receive output suitable to be placed on the right hand side of a formula.

Thank you!

BUML1290
  • 93
  • 6
  • 1
    Please explain in agonizing detail why `y~x:z` should return `1+x`. (In R, `x*z` is `1+x+z+x:z`, _not_ the other way around.) – IRTFM Oct 17 '13 at 23:33
  • @DWin : Sorry, that should say `x*z`, which actually raises another issue -- R interprets `x*z` in a formula as `1+x+z+x:z`, but if you pass `deriv(y~x*z,"x")`, this yields `z` as the derivative, when it should be `1+z`. Regarding `x:z`, it should actually return `z`, not `1+z` (edit made in original question). – BUML1290 Oct 17 '13 at 23:38
  • The `deriv` function is using mathematical definitions of "*" and there is no mathematical definition of ":". Please edit your question to clear up the errors and reframe it to make a legitimate request. It's not helpful to just say a function doesn't behave the way you imagined it. You need to offer reasonable and specific requests for a new function if that is what is needed. – IRTFM Oct 17 '13 at 23:42
  • Right. My question is whether I can coerce it to recognize `x:z`, in an input formula, as `x*z` mathematically (or, further, `x:z:a` as `x*z*a`) and `x*z`, in an input formula, as `1+x+z+x*z` mathematically. – BUML1290 Oct 17 '13 at 23:47
  • My apologies. I wasn't suggesting that there was any problem with `deriv`, simply asking if there's a way to coerce it to recognize the input I'm giving it in the way I want (which is, I think, reasonably intuitive). I've edited the question to reflect this, though I think it has had reasonable and specific requests all along. – BUML1290 Oct 17 '13 at 23:57
  • Having nearly the same question. I rely on `terms`. The `"term.labels"` attribute of the returned value is a character vector. On each element I use `gsub` to replace `:` by `*`, then use `parse`, `deriv` and `eval`. I also must get rid of `I( )` – Yves Jan 26 '22 at 17:01

2 Answers2

2

If you have a formula expression you can work with it using substitute():

substitute( x~x:z+x:y , list(`:`=as.name("*") ) )
x ~ x * z + x * y

And this will let you pass an expression object to substitute with it first being evaluated (which would otherwise not happen since substitute does not evaluate its first argument):

form1 <- expression(x ~ x : z + x : y)
rm(form2)
form2 <- do.call('substitute' , list(form , list(`:`=as.name("*") ) ))
form2
# expression(x ~ x * z + x * y)

This shows how to "reshape" the RHS so that y ~ x:z is handled like ~ x*z by extracting the RHS from its list structure where the tilde operator is being treated as a function and the LHS is the second element in (~ , <LHS>, <RHS>):

 f2<-y~x:z
 substar <- function(form) { 
            do.call('substitute' , list(form , list(`:`=as.name("*") ) )) }
 f3 <- substar(f2)
 deriv(f3[[3]],"x")
 #----------------------
expression({
    .value <- x * z
    .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
    .grad[, "x"] <- z
    attr(.value, "gradient") <- .grad
    .value
})

If you want to work with expressions it may help to understand that they are organized like lists and that the operators are really Lisp-like functions:

> Z <- y~x+I(x^2)
> Z
y ~ x + I(x^2)
> Z[[1]]
`~`
> Z[[2]]
y
> Z[[3]]
x + I(x^2)
> Z[[3]][[1]]
`+`
> Z[[3]][[2]]
x
> Z[[3]][[3]]
I(x^2)
> Z[[3]][[3]][[1]]
I
> Z[[3]][[3]][[2]]
x^2
> Z[[3]][[3]][[2]][[1]]
`^`

If you want to see a function that will traverse an expression tree, the inimitable Gabor Grothendieck constructed one a few years ago in Rhelp: http://markmail.org/message/25lapzv54jc4wfwd?q=list:org%2Er-project%2Er-help+eval+substitute+expression

IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • That's very helpful, thank you. Do you know of a website or other resource that lays out the hierarchy of such an expression (e.g., why is `Z[[1]]` the tilde, `Z[[2]]` y, etc.)? Further, is there some way to extract all terms of a certain type (e.g., all terms that are bracketed by `I()`)? – BUML1290 Oct 18 '13 at 00:04
  • You might take a look at the `findbars` function in https://github.com/lme4/lme4/blob/master/R/utilities.R -- not simple, but it does a similar task ... – Ben Bolker Oct 18 '13 at 00:22
1

the help file of deriv (?deriv)says that expr argument in deriv function is a " A expression or call or (except D) a formula with no lhs" . So you can't use left hand side of the equation in an expression.

On the second part of the question, if I correctly understood your question, you can do something like this: say your rhs is x^2+y^2 and you need to take partial derivative of this expression with x and y: 

myexp <- expression((x^2) + (y^2))
D.sc.x <- D(myexp, "x")
> D.sc.x
2 * x
D.sc.y <- D(myexp, "y")
> D.sc.y
2 * y

In one line: 

lapply(as.list(c("x","y")),function(a)D(myexp,a))
[[1]]
2 * x

[[2]]
2 * y
Metrics
  • 15,172
  • 7
  • 54
  • 83
  • This does not really answer my question. My issue is that I would like to be able to take derivatives of right hand sides of formulas [and by the way, deriv() is able to handle formulas -- try deriv(y~x^2,"x") -- it simply shears off the LHS and takes the derivative of the RHS]. But, if I pass a formula that includes something like I(x^2) or x:z, it doesn't know how to handle those things. Second, I don't want to take two partials. I want to take a second mixed partial, for example. – BUML1290 Oct 17 '13 at 23:25