# Get unbiased derivative from GAM

I have time-dependent data and I'm interested in the gradient at time zero. From theory I expect the gradient of the true relationship to be highest at time zero.

I would like to estimate the derivative at time zero using a GAM. However, the estimates are biased (too low). How can I reduce or preferably eliminate that bias?

#HM function
HM <- function(t, f0, phi, kappa) phi + f0 * exp(-kappa * t)/(-kappa)

t <- 0:180
phi <- 530
kappa <- 0.025
f0 <- 5
Ctrue <- HM(t, f0, phi, kappa)

sdC <- 5

library(mgcv)
library(gratia)

set.seed(1)
res <- replicate(100, {

C <- Ctrue + rnorm(length(t), sd = sdC)

DF <- data.frame(C, t)

fit <- gam(C ~ s(t, bs = "bs", k = 5), data = DF)
#gam.check(fit)

#ggplot(DF, aes(t, C)) +
#  geom_point() +
#  stat_function(fun = \(t) predict(fit, newdata = data.frame(t)))

derivatives(fit, data = data.frame(t = 0))

}, simplify = FALSE)


res <- do.call(rbind, res)

summary(res$derivative) # Min. 1st Qu. Median Mean 3rd Qu. Max. #4.266 4.432 4.508 4.517 4.589 4.773 #True value is 5  • Not good enough for an answer: When you use shape-constrained GAM (package scam), you can force the smooth to be monotonously positive. The derivative function does not work (but there is a stand-in), and while it does not fully remove the bias in my simulation, but at least the replicates now overlap with 5. Here the essential changes in code: library(scam); fit <- scam(C ~ s(t, bs="mpi", k = 5), data = DF) and derivative.scam(fit)$d[1,]. Commented Dec 14, 2023 at 17:10
• @Carsten Thanks. If I also enforce concave smoothers, there is a small further improvement. But there is still bias. (Also, something isn't right with the derivative.scam function. It doesn't apear to give correct derivatives with bs = "micv". I needed to do it manually: diff(predict(fits[[i]], newdata = data.frame(t = c(0, 1e-7)))) / 1e-7) Commented Dec 15, 2023 at 8:58

## 1 Answer

I suggest you may need to increase the potential complexity of the basis (i.e. increase k) and perhaps even extend the range of the penalty beyond zero to get more accurate and correct estimates of the slope at t = 0. I used your same data-generating process and, with these modifications, got a better answer (note that I'm using {marginaleffects} for these calculations, but I'm certain you could do this just as well with {gratia}:

library(mgcv)
library(marginaleffects)
library(ggplot2)
theme_set(theme_bw())

#HM function
HM <- function(t, f0, phi, kappa) phi + f0 * exp(-kappa * t)/(-kappa)

set.seed(1)
t <- 0:180
phi <- 530
kappa <- 0.025
f0 <- 5
Ctrue <- HM(t, f0, phi, kappa)
sdC <- 5
C <- Ctrue + rnorm(length(t), sd = sdC)
DF <- data.frame(C, t)

# The model
fit <- gam(C ~
# Use larger 'k' for more flexibility,
# particularly near the endpoints
s(t, bs = "bs", k = 10),

# Extend penalty beyond zero for more effective
# estimates of derivatives at the bound
knots = list(t = c(-5, 0, 180, 185)),
data = DF)
#> Warning in smooth.construct.bs.smooth.spec(object, dk$$data, dk$$knots): there is
#> *no* information about some basis coefficients

# The estimated function, on the response scale
plot_predictions(fit, condition = 't',
points = 0.5)



# Estimated slopes (using a finer grid of points for a smoother image)
plot_slopes(fit,
newdata = datagrid(t = seq(0, 180, by = 0.2)),
variables = 't', by = 't')



# Estimated slope at t = 0, with appropriate uncertainty envelope
slopes(fit,
newdata = datagrid(t = seq(0, 180, by = 0.2)),
variables = 't', by = 't')[1,]
#>
#>  Term    Contrast t Estimate Std. Error  z Pr(>|z|)     S 2.5 % 97.5 %
#>     t mean(dY/dX) 0     4.72      0.336 14   <0.001 146.0  4.06   5.38
#>
#> Columns: rowid, term, contrast, t, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
#> Type:  response

# Does the estimate at t = 0 differ 'significantly' from 5?
hypotheses(slopes(fit,
newdata = datagrid(t = seq(0, 180, by = 0.2)),
variables = 't', by = 't'),
hypothesis = 5)[1,]
#>
#>  Term    Contrast t Estimate Std. Error      z Pr(>|z|)   S 2.5 % 97.5 %
#>     t mean(dY/dX) 0     4.72      0.336 -0.841      0.4 1.3  4.06   5.38
#>
#> Columns: rowid, term, contrast, t, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
#> Type:  response


Created on 2024-03-25 with reprex v2.1.0

You can read more about how the penalty can be extended in this way using b-splines in Gavin Simpson's helpful post on spline extrapolation