5
$\begingroup$

I'm using gratia_0.7.3

Edits with changes from Gavin Simpson's comments.

Clarification: I am interested in knowing during which years CPUE changes significantly. I am running these models for 3 species in 3 rivers. I have 3 sets of models the first is changes in angler catch rate, the second is change in fisheries size relative abundance, and the third is changes in age-0 relative abundance. Data here is simulated based on angler catch rate (CPUE) of one species in one river. The second and third set of models use electrofishing survey data.

Based on Comments 1 and 2 I refit models with differing distributions and selected the Tweedie distribution with response CPUE in the model formula.

A few questions from comments and other posts I've read...

  • Assumption: Estimation of the first derivative is unbiased

    • I'm not sure how to check/ensure this?
  • To include uncertainty due to smoothing parameter being estimated set unconditional = TRUE.

    • Where should I use unconditional = TRUE in predict(), derivatives(), or both?

First, I need create new_data for the years 2001 to 2021. Below n = 200 rows are generated using seq_min_max() from gratia for years 2001 to 2021.

library(mgcv)
library(gratia)
library(dplyr)
library(magrittr)
library(ggplot2)

new_data <- with(sim_dat,
                 tibble( year = seq_min_max( c( min(year):max(year) ), 200) ) )

Second, I generate a prediction matrix from my current model for new_dat.

  • I am using type="link (as suggested comment 3).

  • I'm going to try including unconditional = TRUE to include uncertainty due to smooth parameter estimation.

    • Question: meant for predict(), derivatives(), or both?
pred_mat <- predict(tw_fit, 
                    newdata = new_data, 
                    type = "link", 
                    se.fit = TRUE,
                    unconditional = TRUE)

Third, I use derivatives() from gratia to get the first order derivative with simultaneous ($95\%$) confidence intervals from my gam() model fit (tw_fit) for the smooth term s(year) over the study period (2001 to 2021 in new_data).

Arguments

type = "central" type of finite difference used. using central derivative so I am not extrapolating too much beyond the start/end of the time series (as suggested comment 4).

n = 200 evaluating the derivative at 200 points suggested in comment 4.

eps = 1 to get derivative over an entire year (suggested comment 4)

unconditional = TRUE (I hope I'm not compounding this by using it in both pred() and derivatives()?)

fd_tw_fit <- derivatives(tw_fit,
                             term = "s(year)",
                             newdata = new_data,
                             order = 1L,
                             type = "central",
                             n = 200,
                             eps = 1,
                             interval = "simultaneous",
                             n_sim = 10000,
                             unconditional = TRUE)
signifD <- function(x, d, upper, lower, eval = 0) {
  miss <- upper > eval & lower < eval
  incr <- decr <- x
  want <- d > eval
  incr[!want | miss] <- NA
  want <- d < eval
  decr[!want | miss] <- NA
  list(incr = incr, decr = decr)
  }
m2.dsig <- signifD(x = pred_mat$fit, 
               d = fd_tw_fit$derivative,
                   upper = fd_tw_fit$upper, 
               lower = fd_tw_fit$lower,
                   eval = 0)
incr <- which( !(is.na(m2.dsig$incr)) ) %>%
  fd_tw_fit[.,]

decr <- which( !(is.na(m2.dsig$decr)) ) %>%
  fd_tw_fit[.,] 

Plotting the derivative and simultaneous confidence intervals.

draw(fd_tw_fit,
     alpha = 0.2) +
  theme(axis.text.x = element_text(angle = -45)) +
  scale_x_continuous(breaks = seq(2001, 2021, 1),
                     labels = seq(2001, 2021, 1)) +
  labs(title = "Simultaneous CI first derivative tw() GAM") +
  geom_hline(yintercept = 0)

First Derivative Plot

Plotting GAM fit on response scale Transforming predicted data to the response scale with exp() (comment 3).

# Calculate upper and lower confidence Intervals
upper_2se <- (pred_mat$fit) + (pred_mat$se.fit*2)  
lwr_2se <- (pred_mat$fit) - (pred_mat$se.fit*2)

# Transform to Response scale
upper_2se <- exp(upper_2se)
lwr_2se <- exp(lwr_2se)
pred_fit <- exp(pred_mat$fit)
new_data %<>% dplyr::mutate(pred = pred_fit,
                        lwr_2se = lwr_2se,
                        upr_2se = upper_2se,
                        change_point = ifelse(new_data$year %in% incr$data,
                                              "incr",
                                              ifelse(new_data$year %in% decr$data,
                                                         "decr",
                                                         "NS") ) %>%
                              factor())
raw_dat <- sim_dat %>%
  select(year,
         CPUE) %>%
  dplyr::mutate(change_point = "NS" %>% factor())

create ggplot Response scale with Change Points

ggplot(data = new_data, aes(x = year, y = pred, color = change_point))+
  geom_path(aes(group = 1),
            size = 1.25) +
  scale_color_manual(values = c("NS" = "black",
                                "incr" = "blue",
                                "decr" = "red")) +
  geom_ribbon(aes(ymin = lwr_2se, ymax = upr_2se),
              colour = "grey",
              alpha = 0.2,
              linetype = "blank") +
  geom_point(data = raw_dat,
             aes(x = year, y = CPUE, color = change_point),
             size = 1.5,
             alpha = 0.2,
             colour = "black")

Response scale Change Point Plot

New simulated data

sim_dat <- structure(list(year = c(2001, 2001, 2001, 2001, 2001, 2001, 2001, 
2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 
2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 2001, 
2001, 2001, 2001, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 
2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 
2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2002, 
2002, 2002, 2003, 2003, 2003, 2003, 2003, 2003, 2003, 2003, 2003, 
2003, 2003, 2003, 2003, 2003, 2003, 2003, 2003, 2003, 2003, 2003, 
2003, 2003, 2003, 2003, 2003, 2003, 2003, 2003, 2004, 2004, 2004, 
2004, 2004, 2004, 2004, 2004, 2004, 2004, 2004, 2004, 2004, 2004, 
2004, 2004, 2004, 2004, 2004, 2004, 2004, 2004, 2005, 2005, 2005, 
2005, 2005, 2005, 2005, 2005, 2005, 2005, 2005, 2005, 2005, 2006, 
2006, 2006, 2006, 2006, 2006, 2006, 2006, 2006, 2006, 2006, 2006, 
2006, 2006, 2006, 2006, 2006, 2007, 2007, 2007, 2007, 2007, 2007, 
2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 
2007, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 
2008, 2008, 2008, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 
2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2010, 2010, 2010, 
2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 
2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 
2011, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 
2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 
2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 
2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 
2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 
2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 
2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 
2014, 2014, 2014, 2014, 2014, 2014, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2017, 
2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 
2017, 2017, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 
2018, 2018, 2018, 2018, 2018, 2018, 2019, 2019, 2019, 2019, 2019, 
2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 
2019, 2019, 2019, 2019, 2019, 2019, 2020, 2020, 2020, 2020, 2020, 
2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 
2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 
2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021
), eff_hrs = c(232, 90, 696, 702, 192, 423, 280, 117, 120, 54, 
54, 1216, 240, 108, 1120, 437, 527.5, 95, 48, 702, 144, 45, 692, 
54, 466, 612, 48, 152, 810, 48, 5112, 108, 934, 162, 56, 1097, 
48, 64, 720, 405, 846, 2466, 289, 360, 567, 189, 702, 1440, 289, 
90, 928, 243, 1512, 98, 738, 48, 486, 64, 88, 232, 810, 128, 
304, 1216, 64, 136, 648, 99, 555, 1476, 315, 1392, 432, 555, 
1785, 540, 104, 656, 1288, 135, 516.5, 1785, 666, 112, 135, 72, 
378, 232, 144, 180, 50, 45, 126, 162, 112, 56, 1328, 252, 486, 
32, 585, 1008, 1926, 40, 648, 333, 119, 144, 32, 720, 117, 1926, 
648, 40, 559, 400, 352, 288, 1404, 468, 1404, 180, 468, 774, 
162, 162, 99, 216, 288, 936, 1156, 1071, 304, 1116, 1071, 1292, 
288, 324, 256, 992, 1216, 243, 630, 1296, 296, 702, 40, 99, 328, 
1072, 216, 48, 972, 1206, 80, 48, 64, 225, 70, 144, 765, 904, 
38.5, 14, 220, 270, 1098, 756, 756, 200, 680, 595, 266, 333, 
1782, 99, 32, 48, 252, 704, 351, 864, 972, 30, 48, 414, 2000, 
342, 112, 736, 144, 198, 96, 1984, 2232, 653, 570, 48, 342, 117, 
135, 88, 396, 135, 135, 108, 918, 2000, 48, 64, 117, 846, 639, 
234, 936, 144, 72, 468, 333, 72, 54, 459, 828, 684, 56, 752, 
42, 1022, 846, 944, 744, 200, 324, 900, 584, 80, 360, 540, 76.5, 
112, 2432, 99, 20, 288, 320, 408, 1008, 35, 112.5, 306, 792, 
378, 450, 1080, 24, 234, 288, 96, 68, 72, 243, 110, 75, 126, 
176, 49.5, 94, 64, 222, 648, 225, 49.5, 50, 90, 65, 90, 70, 81, 
264, 70, 558, 874, 128, 904, 371, 896, 96, 2340, 943.5, 1008, 
31.5, 32, 52, 504, 216, 28, 52, 32, 40, 560, 136, 180, 240, 720, 
64, 288, 800, 52, 736, 32, 32, 24, 288, 288, 352, 810, 384, 752, 
432, 810, 252, 720, 848, 594, 252, 136, 944, 368, 108, 80, 60, 
768, 405, 352, 1008, 70, 1009, 384, 85, 40, 112, 40, 1120, 72, 
4000, 324, 255, 55, 128, 32, 128, 176, 50, 928, 24, 240, 40, 
64, 128, 28, 184, 248, 32, 592, 20, 72, 800, 98, 48, 1170, 297, 
42, 30, 66, 44, 88, 1080, 96, 112, 336, 736, 192, 624, 369, 1088, 
176, 144, 104, 2196, 432, 128, 144, 585, 120, 216, 45, 2250, 
80, 816, 992, 256, 240, 540, 72, 45, 1170, 208, 200, 101, 224, 
387, 312, 18, 45), num = c(6, 0, 8, 4, 12, 16, 0, 18, 9, 86, 
8, 1, 6, 26, 4, 48, 21, 10, 0, 26, 3, 85, 9, 98, 5, 14, 21, 40, 
1, 15, 7, 53, 32, 8, 35, 20, 9, 36, 66, 39, 10, 37, 14, 14, 4, 
8, 11, 5, 4, 13, 3, 106, 10, 19, 125, 66, 56, 3, 2, 37, 3, 23, 
8, 23, 37, 191, 44, 14, 113, 4, 27, 6, 46, 11, 7, 51, 8, 3, 72, 
21, 11, 54, 24, 24, 47, 58, 3, 20, 5, 5, 23, 14, 2, 1, 2, 5, 
0, 0, 9, 10, 9, 1, 41, 16, 0, 4, 2, 1, 9, 2, 1, 40, 3, 0, 5, 
3, 10, 20, 9, 19, 14, 3, 2, 8, 9, 19, 7, 15, 40, 41, 100, 58, 
19, 42, 12, 71, 85, 14, 29, 2, 33, 34, 75, 72, 21, 59, 2, 17, 
10, 16, 7, 53, 6, 3, 11, 4, 3, 3, 80, 58, 3, 5, 51, 10, 3, 1, 
18, 0, 45, 14, 56, 0, 28, 53, 45, 29, 63, 79, 2, 13, 22, 34, 
9, 54, 122, 41, 105, 105, 0, 14, 8, 35, 38, 77, 89, 37, 202, 
87, 10, 12, 4, 62, 10, 20, 12, 225, 18, 74, 0, 177, 85, 16, 16, 
12, 30, 34, 22, 118, 195, 154, 87, 21, 0, 180, 61, 70, 58, 124, 
76, 58, 21, 123, 19, 85, 60, 67, 72, 175, 17, 13, 48, 5, 9, 20, 
28, 9, 12, 190, 33, 65, 6, 106, 24, 13, 163, 25, 47, 18, 41, 
103, 35, 16, 0, 12, 26, 188, 12, 19, 15, 21, 127, 39, 58, 49, 
73, 183, 181, 208, 110, 131, 8, 10, 28, 60, 4, 20, 16, 63, 2, 
40, 8, 4, 18, 133, 8, 1, 26, 9, 6, 1, 54, 7, 8, 18, 15, 84, 3, 
20, 51, 24, 18, 46, 39, 75, 62, 90, 85, 90, 64, 49, 65, 99, 97, 
81, 65, 31, 91, 6, 34, 44, 75, 34, 10, 20, 27, 15, 15, 26, 25, 
25, 57, 5, 3, 4, 34, 25, 138, 1, 12, 11, 17, 19, 3, 0, 35, 5, 
3, 2, 4, 70, 5, 12, 73, 4, 8, 10, 6, 44, 4, 31, 18, 31, 8, 7, 
27, 9, 20, 24, 77, 64, 1, 18, 29, 74, 27, 4, 53, 60, 40, 37, 
33, 45, 133, 101, 6, 87, 4, 83, 135, 18, 50, 10, 92, 23, 43, 
135, 61, 86, 6, 75, 83, 177, 57, 5), CPUE = c(0.0638297872340425, 
0.0690789473684211, 0.180555555555556, 0.0833333333333333, 0.021505376344086, 
0.0211640211640212, 0.0065359477124183, 0.0305555555555556, 0.0375, 
0.075, 0.00997150997150997, 0.0446428571428571, 0.0493827160493827, 
0.100473933649289, 0.125, 0.0555555555555556, 0.102564102564103, 
0.0229885057471264, 0.0749536178107607, 0.109375, 0.037037037037037, 
0.0718954248366013, 0.09375, 0.0683760683760684, 0.0555555555555556, 
0.0555555555555556, 0.175, 0.0592592592592593, 0.0648148148148148, 
0.0276854928017719, 0.0245901639344262, 0.037037037037037, 0.0972222222222222, 
0.0803443328550933, 0.075, 0.0129310344827586, 0.0608108108108108, 
0.0963597430406852, 0.16, 0.0816326530612245, 0.0375, 0.0584795321637427, 
0.0501355013550136, 0.2, 0.0617283950617284, 0.111111111111111, 
0.274509803921569, 0.127941176470588, 0.159090909090909, 0.037037037037037, 
0.0522875816993464, 0.0790598290598291, 0.0582010582010582, 0.1875, 
0.0625, 0.0375, 0.0493827160493827, 0.0596107055961071, 0.0346020761245675, 
0.0474020054694622, 0.14375, 0.1125, 0.0652777777777778, 0.046875, 
0.0862068965517241, 0.018018018018018, 0.0147058823529412, 0.0285714285714286, 
0.0471380471380471, 0.0444444444444444, 0.107142857142857, 0.0384615384615385, 
0.0611111111111111, 0.0990853658536585, 0.062888198757764, 0.0315315315315315, 
0.137212643678161, 0.0925925925925926, 0.111111111111111, 0.0851851851851852, 
0.0318428184281843, 0.0517241379310345, 0.109375, 0.0476190476190476, 
0.162162162162162, 0.111111111111111, 0.0101010101010101, 0.036036036036036, 
0.0285714285714286, 0.0126984126984127, 0.0493827160493827, 0.0166666666666667, 
0.0138888888888889, 0.0756302521008403, 0.0476190476190476, 0.0138888888888889, 
0.0154320987654321, 0.0265625, 0.0238095238095238, 0.00467289719626168, 
0.0488095238095238, 0.0246913580246914, 0.00793650793650794, 
0.03125, 0.0169082125603865, 0.0113636363636364, 0.0123456790123457, 
0.0188253012048193, 0.040983606557377, 0.00555555555555556, 0, 
0.0175438596491228, 0.0125, 0.00467289719626168, 0.0142450142450142, 
0.0111111111111111, 0.0250447227191413, 0.03125, 0.0142045454545455, 
0.0303030303030303, 0.0170940170940171, 0.0432098765432099, 0.0555555555555556, 
0.00641025641025641, 0.025, 0.0135327635327635, 0.0245478036175711, 
0.0672268907563025, 0.030982905982906, 0.0616776315789474, 0.0634920634920635, 
0.0659722222222222, 0.0432098765432099, 0.00657894736842105, 
0.0694444444444444, 0.0416666666666667, 0.0263157894736842, 0.0662931839402428, 
0.0856854838709677, 0.172839506172839, 0.0447530864197531, 0.16015625, 
0.0285467128027682, 0.0896057347670251, 0.0625, 0.00925925925925926, 
0.333333333333333, 0.0444444444444444, 0.0196078431372549, 0.075, 
0.109375, 0.0740740740740741, 0.09375, 0.0833333333333333, 0.00712250712250712, 
0.0299145299145299, 0.101307189542484, 0.0416666666666667, 0.0127659574468085, 
0.0640243902439024, 0.0142450142450142, 0.0689300411522634, 0.138888888888889, 
0.0693069306930693, 0, 0.0555555555555556, 0.0977443609022556, 
0.183060109289617, 0, 0.141203703703704, 0.0324074074074074, 
0.1, 0.0535714285714286, 0.0277777777777778, 0.0130813953488372, 
0.0994152046783626, 0.154761904761905, 0.0674603174603175, 0.126436781609195, 
0.166666666666667, 0.0729166666666667, 0.0990338164251208, 0.142663043478261, 
0.090625, 0.111111111111111, 0.0803571428571429, 0.171052631578947, 
0.0484330484330484, 0.139204545454545, 0.0936213991769547, 0.118055555555556, 
0.0416666666666667, 0.101814516129032, 0.0555555555555556, 0.0458937198067633, 
0.0310077519379845, 0.0772569444444444, 0.0727969348659004, 0.0751953125, 
0.133333333333333, 0.111111111111111, 0.0436681222707424, 0.0551301684532925, 
0.0833333333333333, 0.0625, 0.13302034428795, 0.244444444444444, 
0.261363636363636, 0.180555555555556, 0.472222222222222, 0.0681818181818182, 
0.276353276353276, 0.375, 0.28125, 0.316239316239316, 0.115, 
0, 0.0847953216374269, 0.339285714285714, 0.395833333333333, 
0.132716049382716, 0.479166666666667, 0.201058201058201, 0.222222222222222, 
0.2, 0.137096774193548, 0.235042735042735, 0.201388888888889, 
0.211352657004831, 0.388888888888889, 0.35, 0.25, 0.182795698924731, 
0.18018018018018, 0.25, 0.158666666666667, 0.2, 0.245, 0.6, 0.299145299145299, 
0.244444444444444, 0.37037037037037, 0.142857142857143, 0, 0.137345679012346, 
0.119047619047619, 0.277777777777778, 0.221590909090909, 0.424242424242424, 
0.09375, 0.260416666666667, 0.305555555555556, 0.228535353535354, 
0.17989417989418, 0.0990990990990991, 0.2, 0.4, 0.282608695652174, 
0.106666666666667, 0.166666666666667, 0.144230769230769, 0.242990654205607, 
0.0696022727272727, 0.171717171717172, 0.457516339869281, 0.189964157706093, 
0.311274509803922, 0.129496402877698, 0.253472222222222, 0.177295918367347, 
0.107954545454545, 0.346666666666667, 0.202020202020202, 0.363636363636364, 
0.56, 0.282738095238095, 0.148148148148148, 0.2, 0.270833333333333, 
0.189542483660131, 0.186507936507937, 0.14537037037037, 0.0987654320987654, 
0.101709401709402, 0.0208333333333333, 0.0667726550079491, 0.126984126984127, 
0.204545454545455, 0.0892857142857143, 0.0539083557951483, 0.03125, 
0.0491071428571429, 0.0568376068376068, 0.0961538461538462, 0.346153846153846, 
0.0714285714285714, 0.222222222222222, 0.12037037037037, 0.0515734265734266, 
0.0833333333333333, 0.0833333333333333, 0.114130434782609, 0.125, 
0.0821428571428571, 0.0714285714285714, 0.0944444444444444, 0.375, 
0.103333333333333, 0.134408602150538, 0.125, 0.109375, 0.0833333333333333, 
0.0511363636363636, 0.101190476190476, 0.0524691358024691, 0.134722222222222, 
0.0754716981132075, 0.127604166666667, 0.143097643097643, 0.208333333333333, 
0.11968085106383, 0.257936507936508, 0.1, 0.257936507936508, 
0.122222222222222, 0.12037037037037, 0.190775681341719, 0.119791666666667, 
0.0677083333333333, 0.25, 0.0294117647058824, 0.0880681818181818, 
0.0756578947368421, 0.0572916666666667, 0.227941176470588, 0.126984126984127, 
0.4, 0.06, 0.0679347826086956, 0.0727272727272727, 0.140625, 
0.1, 0.0196078431372549, 0.0797413793103448, 0.0475, 0.0625, 
0.0222222222222222, 0.0617283950617284, 0.058641975308642, 0.223214285714286, 
0.0661764705882353, 0.2, 0.125, 0.0277777777777778, 0.0357142857142857, 
0.032258064516129, 0.0816326530612245, 0.104166666666667, 0.214285714285714, 
0.075, 0.0422297297297297, 0.2421875, 0.166666666666667, 0.1125, 
0.277777777777778, 0.35, 0.0833333333333333, 0.183333333333333, 
0.0134680134680135, 0.0978260869565217, 0.0608888888888889, 0.0925925925925926, 
0.125, 0.241071428571429, 0.0497685185185185, 0.149590163934426, 
0.0902777777777778, 0.0677966101694915, 0.25, 0.178197064989518, 
0.204545454545455, 0.09, 0.339285714285714, 0.015625, 0.0542005420054201, 
0.104166666666667, 0.0303308823529412, 0.00925925925925926, 0.0605646630236794, 
0.110294117647059, 0.0815217391304348, 0.119047619047619, 0.0727040816326531, 
0.165441176470588, 0.0666666666666667, 0.248161764705882, 0.129166666666667, 
0.0769230769230769, 0.17037037037037, 0.174358974358974, 0.306930693069307, 
0.375, 0.271875, 0.0866935483870968, 0.166666666666667, 0.133333333333333, 
0.240384615384615, 0.161835748792271, 0.166666666666667, 0.0786666666666667, 
0.140939597315436, 0.152573529411765, 0.231060606060606, 0.133333333333333
)), class = "data.frame", row.names = c(NA, -414L))
$\endgroup$
3
  • $\begingroup$ What do you (or your committee) mean by a "test statistic". The derivative itself is a (scaled) difference of means (the respective means being the the two data points over which we computed the derivative. You could use that as the test statistic, and if you want a p-value, you could invert the confidence/simultaneous interval to get one. You could also divide the derivative by it's standard error to get a t-like statistic. $\endgroup$ Jan 6, 2023 at 16:11
  • $\begingroup$ Or are you interested in the timing of the changepoints? I.e. not the uncertainty in the response or the derivative but the uncertainty in year as to when the changes took place? We could do that; its a bit tedious but quite easy once you know how; but no point answering that if it isn't the actual question you/they want an answer to. $\endgroup$ Jan 6, 2023 at 16:12
  • $\begingroup$ @GavinSimpson Thank you for your explanation. I am interested in when the changes take place. I was using gamm() because I had a ARMA() correlation structure in there but wasn't necessary with the data I simulated for the example and left it out to not add more complexity. I'll try your suggestions and post back. Thanks again $\endgroup$
    – THATguy
    Jan 9, 2023 at 21:24

1 Answer 1

1
$\begingroup$

There are a few things wrong here:

  1. I'd fit this with the tw() family using gam() not gamm(). You have some 0s, and these aren't allowed in the Gamma family so you are adding 1 to your response, but I don't see any reason why you need to use gamm() here.

  2. Better, assuming you have the effort, variable would be to model these as actual counts with poisson() or nb() families and include the effort via an offset term - add + offset(log(effort)) to your model formula, where effort is the variable that normalises the counts to catch per unit effort. Absent that variable, the tw() is the best option.

  3. You are computing the confidence interval (upper_2se etc) using predictions from the response scale and that is incorrect: as your plot shows, the interval contains negative values which just doesn't make sense. instead, use type = "link" when predicting, and then form the confidence interval as you do it now, but then backtransform the upper and lower limits and the fitted values to the response scale using the inverse of the link function (in all the cases mentioned — Gamma, tw, nb poisson — the link function is log(), so you'd want exp() for the inverse. That way you'll get a confidence interval that respects the properties of the data.

  4. I'd compute the derivative over a much finer grid than you do — you are estimating the derivative over a very small interval around the integer year values.

    This is actually causing you some problems as shown in the last plot; each of the blue spikes should be associated with a red decrease, but because of the crude grid of points you are evaluating the derivative at, you simply aren't seeing the significant decreases. That last red section is also likely spurious.

    If you want an instantaneous measure of the derivative you should do like you do when you plot the smooth and do it for a finer grid, say 100-200 points over the range min(year) to max(year) with seq(min(year), max(year), length = 200) say.

    Or, if you want the derivative over the entire year (?) you would have to increase eps to be some much larger value, like 1, and I would suggest you do a central derivative so you aren't extrapolating too much beyond the start/end of the time series.

Remember these are derivatives on the link (log) scale; they aren't response level changes. Response derivatives is a bit more involved but should be in {gratia} in a few weeks.

I'll try to work up a full example from your data, but I also have some questions which I'll post as comments first. Consider this a WIP.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.