Problem calling deriv function on fd objects

Question

I'm studying the book "Functional Data Analysis with R and MATLAB" using R and the R package fda and I have a small problem replicating Figure 1.12. It's very easy to fix it (see below) but I don't understand the underlying issue.

A reprex of the problem:

library(fda)
#> Loading required package: splines
#> Loading required package: fds
#> Loading required package: rainbow
#> Loading required package: MASS
#> Loading required package: pcaPP
#> Loading required package: RCurl
#> Loading required package: deSolve
#> 
#> Attaching package: 'fda'
#> The following object is masked from 'package:graphics':
#> 
#>     matplot

# Create an fd object represeting the weather in 4 canadian stations
fig1.11Stns = c('Montreal', 'Edmonton', 'Pr. Rupert', 'Resolute')
fig1.11Temp = CanadianWeather$dailyAv[
  , 
  fig1.11Stns, 
  'Temperature.C'
]

Temp.fourier   = create.fourier.basis(c(0, 365), 13)
fig1.11Temp.fd = Data2fd(day.5, fig1.11Temp, Temp.fourier)
class(fig1.11Temp.fd)
#> [1] "fd"

# Define a linear differential operator
yearRng      = c(0,365)
Lcoef        = matrix(c(0,(2*pi/365)^2,0),1,3)
harmaccelLfd = vec2Lfd(Lcoef, yearRng)

# Compute the derivative according to the linear differential operator
deriv(fig1.11Temp.fd, harmaccelLfd) # fails
#> Error in deriv.default(fig1.11Temp.fd, harmaccelLfd): invalid variable names

# However, when I explicitly call the S3 method, everything works fine
deriv.fd(fig1.11Temp.fd, harmaccelLfd) |> plot()

#> [1] "done"

^{Created on 2024-08-26 with reprex v2.0.2}

Why is there an error when I call the generic deriv but there is no error when I explicitly call the S3 method for fd objects?

A few notes:

1. the code (up to the last line) is taken from the companion scripts of the book, which are available inside the folder named scripts within the downloaded R package: system.file("scripts", package = "fda"). In particular, I'm using the script named fdarm-ch01.R.
2. I manually added the last line (i.e. deriv.fd instead of deriv).

Answer 1

It seems this is because the {fda} package does not register its deriv.fd() S3 method as one might expect.

If you look at the NAMESPACE, you will find this line:

S3method(deriv, fd)

The above NAMESPACE line is registering fda::deriv.fd() with the fda:::deriv() generic because fda:::deriv() is found before stats::deriv() "as seen from the package" (see this blog post). However, fda:::deriv() is not exported from the package.

The effect of this registration and lack of export is that the fda:::deriv() generic will not show up on the search path via library(fda). Therefore calls to deriv() as in your code use the stats::deriv() generic which does not know about the fda::deriv.fd() method.

If the fda:::deriv() generic was exported, when you library(fda) the fda::deriv() generic would take precedence over the pre-existing stats::deriv() generic, and R could then dispatch the call to fda::deriv.fd() based on the class of the first argument according to S3 dispatch.

I am unsure if this is intended behavior or not on the part of the package authors. To me the presence of a fda:::deriv() generic suggests the authors intended to export. As Ben suggests, you might consider contacting them for clarification.

There are a number of ways to change this behavior. Here are seven off the top of my head to illustrate the S3 registration and search path issues at hand:

1. Call fda::deriv.fd() S3 Method Directly
2. Call fda:::deriv() Generic Directly
3. Copies in Global Environment
4. Register fda::deriv.fd() with stats::deriv()
5. Fork and Export fda:::deriv() Generic
6. Fork and Remove fda:::derive() Generic
7. Fork and Register fda::deriv.fd() with stats::deriv()

1. Call fda::deriv.fd() S3 Method Directly

Bypass the S3 issues altogether and call fda::deriv.fd() directly as you did.

library(fda)

# Create an fd object represeting the weather in 4 canadian stations
fig1.11Stns = c('Montreal', 'Edmonton', 'Pr. Rupert', 'Resolute')
fig1.11Temp = CanadianWeather$dailyAv[
  , 
  fig1.11Stns, 
  'Temperature.C'
]

Temp.fourier   = create.fourier.basis(c(0, 365), 13)
fig1.11Temp.fd = Data2fd(day.5, fig1.11Temp, Temp.fourier)
class(fig1.11Temp.fd)
#> [1] "fd"

# Define a linear differential operator
yearRng      = c(0,365)
Lcoef        = matrix(c(0,(2*pi/365)^2,0),1,3)
harmaccelLfd = vec2Lfd(Lcoef, yearRng)

# Compute the derivative according to the linear differential operator
deriv(fig1.11Temp.fd, harmaccelLfd) # fails
#> Error in deriv.default(fig1.11Temp.fd, harmaccelLfd): invalid variable names

# call fda::deriv.fd() directly
deriv.fd(fig1.11Temp.fd, harmaccelLfd) |> plot()

#> [1] "done"

^{Created on 2024-08-28 with reprex v2.1.1}

2. Call fda:::deriv() Generic Directly

You could call fda:::deriv() directly:

library(fda)

# Create an fd object represeting the weather in 4 canadian stations
fig1.11Stns = c('Montreal', 'Edmonton', 'Pr. Rupert', 'Resolute')
fig1.11Temp = CanadianWeather$dailyAv[
  , 
  fig1.11Stns, 
  'Temperature.C'
]

Temp.fourier   = create.fourier.basis(c(0, 365), 13)
fig1.11Temp.fd = Data2fd(day.5, fig1.11Temp, Temp.fourier)
class(fig1.11Temp.fd)
#> [1] "fd"

# Define a linear differential operator
yearRng      = c(0,365)
Lcoef        = matrix(c(0,(2*pi/365)^2,0),1,3)
harmaccelLfd = vec2Lfd(Lcoef, yearRng)

# Compute the derivative according to the linear differential operator
deriv(fig1.11Temp.fd, harmaccelLfd) # fails
#> Error in deriv.default(fig1.11Temp.fd, harmaccelLfd): invalid variable names

# call fda:::deriv() generic
# for which fda::deriv.fd() has been registered via NAMESPACE
fda:::deriv(fig1.11Temp.fd, harmaccelLfd) |> plot()

#> [1] "done"

^{Created on 2024-08-28 with reprex v2.1.1}

3. Copies in Global Environment

You could make copies of fda:::deriv() and fda::deriv.fd() in global:

library(fda)

# Create an fd object represeting the weather in 4 canadian stations
fig1.11Stns = c('Montreal', 'Edmonton', 'Pr. Rupert', 'Resolute')
fig1.11Temp = CanadianWeather$dailyAv[
  , 
  fig1.11Stns, 
  'Temperature.C'
]

Temp.fourier   = create.fourier.basis(c(0, 365), 13)
fig1.11Temp.fd = Data2fd(day.5, fig1.11Temp, Temp.fourier)
class(fig1.11Temp.fd)
#> [1] "fd"

# Define a linear differential operator
yearRng      = c(0,365)
Lcoef        = matrix(c(0,(2*pi/365)^2,0),1,3)
harmaccelLfd = vec2Lfd(Lcoef, yearRng)

# Compute the derivative according to the linear differential operator
deriv(fig1.11Temp.fd, harmaccelLfd) # fails
#> Error in deriv.default(fig1.11Temp.fd, harmaccelLfd): invalid variable names

# place a copy of fda:::deriv() generic in global
deriv <- fda:::deriv
# and a copy of fda::deriv.fd() so that the new generic can find the fd method
deriv.fd <- fda::deriv.fd

# call your global deriv() generic
# - the global deriv() generic is found before stats::deriv() generic due to being defined in global environment
# - the global deriv() generic can find deriv.fd() due to it also being defined in the global environment
deriv(fig1.11Temp.fd, harmaccelLfd) |> plot()

#> [1] "done"

^{Created on 2024-08-28 with reprex v2.1.1}

4. Register fda::deriv.fd() with stats::deriv()

You could use registerS3method("deriv", "fd", deriv.fd) to register fda::deriv.fd() with stats::deriv():

library(fda)

# Create an fd object represeting the weather in 4 canadian stations
fig1.11Stns = c('Montreal', 'Edmonton', 'Pr. Rupert', 'Resolute')
fig1.11Temp = CanadianWeather$dailyAv[
  , 
  fig1.11Stns, 
  'Temperature.C'
]

Temp.fourier   = create.fourier.basis(c(0, 365), 13)
fig1.11Temp.fd = Data2fd(day.5, fig1.11Temp, Temp.fourier)
class(fig1.11Temp.fd)
#> [1] "fd"

# Define a linear differential operator
yearRng      = c(0,365)
Lcoef        = matrix(c(0,(2*pi/365)^2,0),1,3)
harmaccelLfd = vec2Lfd(Lcoef, yearRng)

# Compute the derivative according to the linear differential operator
deriv(fig1.11Temp.fd, harmaccelLfd) # fails
#> Error in deriv.default(fig1.11Temp.fd, harmaccelLfd): invalid variable names

# register fda::deriv.fd() S3 method with stats::deriv() generic
registerS3method("deriv", "fd", deriv.fd)

# call stats::deriv() generic
# - stats::deriv() generic is found before fda:::deriv() generic because fda:::deriv() is not exported and therefore not on search path
# - fda::deriv.fd() is registered to stats::deriv() via registerS3method()
deriv(fig1.11Temp.fd, harmaccelLfd) |> plot()

#> [1] "done"

^{Created on 2024-08-28 with reprex v2.1.1}

5. Fork and Export fda:::deriv() Generic

You could fork the source and add export(deriv) to the namespace file.

I have done so for you here.

# install patched version
remotes::install_github("the-mad-statter/fda@export-deriv-generic")

library(fda)

# Create an fd object represeting the weather in 4 canadian stations
fig1.11Stns = c('Montreal', 'Edmonton', 'Pr. Rupert', 'Resolute')
fig1.11Temp = CanadianWeather$dailyAv[
  , 
  fig1.11Stns, 
  'Temperature.C'
]

Temp.fourier   = create.fourier.basis(c(0, 365), 13)
fig1.11Temp.fd = Data2fd(day.5, fig1.11Temp, Temp.fourier)
class(fig1.11Temp.fd)
#> [1] "fd"

# Define a linear differential operator
yearRng      = c(0,365)
Lcoef        = matrix(c(0,(2*pi/365)^2,0),1,3)
harmaccelLfd = vec2Lfd(Lcoef, yearRng)

# call fda::deriv() generic
# - fda::deriv() generic is found before stats::deriv() generic because fda::deriv() generic is now exported via NAMESPACE and therefore on the search path before stats::deriv() generic due to library(fda)
# - fda::deriv.fd() is registered to fda::deriv() via NAMESPACE
deriv(fig1.11Temp.fd, harmaccelLfd) |> plot()

#> [1] "done"

^{Created on 2024-08-28 with reprex v2.1.1}

6. Fork and Remove fda:::derive() Generic

You could fork the source and remove the fda::deriv() generic from the package.

I have done so for you here.

# install patched version
remotes::install_github("the-mad-statter/fda@remove-generic")

library(fda)

# Create an fd object represeting the weather in 4 canadian stations
fig1.11Stns = c('Montreal', 'Edmonton', 'Pr. Rupert', 'Resolute')
fig1.11Temp = CanadianWeather$dailyAv[
  , 
  fig1.11Stns, 
  'Temperature.C'
]

Temp.fourier   = create.fourier.basis(c(0, 365), 13)
fig1.11Temp.fd = Data2fd(day.5, fig1.11Temp, Temp.fourier)
class(fig1.11Temp.fd)
#> [1] "fd"

# Define a linear differential operator
yearRng      = c(0,365)
Lcoef        = matrix(c(0,(2*pi/365)^2,0),1,3)
harmaccelLfd = vec2Lfd(Lcoef, yearRng)

# call stats::deriv() generic
# - fda:::deriv() has been removed from package
# - fda::deriv.fd() is registered to stats::deriv() via NAMESPACE
deriv(fig1.11Temp.fd, harmaccelLfd) |> plot()

#> [1] "done"

^{Created on 2024-08-28 with reprex v2.1.1}

7. Fork and Register fda::deriv.fd() with stats::deriv()

You could fork the source and change S3method(deriv, fd) to S3method(stats::deriv, fd) in the NAMESPACE file.

I have done so for you here.

# install patched version
remotes::install_github("the-mad-statter/fda@register-with-stats")

library(fda)

# Create an fd object represeting the weather in 4 canadian stations
fig1.11Stns = c('Montreal', 'Edmonton', 'Pr. Rupert', 'Resolute')
fig1.11Temp = CanadianWeather$dailyAv[
  , 
  fig1.11Stns, 
  'Temperature.C'
]

Temp.fourier   = create.fourier.basis(c(0, 365), 13)
fig1.11Temp.fd = Data2fd(day.5, fig1.11Temp, Temp.fourier)
class(fig1.11Temp.fd)
#> [1] "fd"

# Define a linear differential operator
yearRng      = c(0,365)
Lcoef        = matrix(c(0,(2*pi/365)^2,0),1,3)
harmaccelLfd = vec2Lfd(Lcoef, yearRng)

# call stats::deriv() generic
# - fda:::deriv() is not exported and therefore not on search path
# - fda::deriv.fd() is registered to stats::deriv() via NAMESPACE
deriv(fig1.11Temp.fd, harmaccelLfd) |> plot()

#> [1] "done"

^{Created on 2024-08-28 with reprex v2.1.1}