search

vincentarelbundock / marginaleffects

R package to compute and plot predictions, slopes, marginal means, and comparisons (contrasts, risk ratios, odds, etc.) for over 100 classes of statistical and ML models. Conduct linear and non-linear hypothesis tests, or equivalence tests. Calculate uncertainty estimates using the delta method, bootstrapping, or simulation-based inference
https://marginaleffects.com
Other
435 stars 45 forks source link

Feature request: compute and compare peak of the slope and the corresponding timing for brms model. #931

Closed Sandhu-SS closed 11 months ago

Sandhu-SS commented 11 months ago

In many disciplines (such as epidemiology), it is often of interest to find the peak velocity (PV) and the age at the peak velocity (APV).

It would be great if marginaleffects allows for estimating PV and APV for each posterior draw and then summaries across draws to get the estimate and corresponding uncertainty (i.e, CI).

I think it would be straightforward to implement it because already the slopes function available in the marginaleffects estimates the velocity curve.

A new function (which could be named as findpeak or getpeak etc.) with two arguments, one the velocity curve returned from the slopes for each draw and another time/x that is used in estimating the velocity curve, can used to get the required PV and APV estimates for each draw that are then summarized.

It would be great if users are allowed to pass their own function just like byfun argument in the marginaleffects::predictions(). However, unlike the byfun that is used to summarize the estimates obtained from the posterior draws, this findpeak will be used for each draw.

As an example, an additional argument findpeak can be added to the marginaleffects::slopes() that accepts a predefined function funx such as findpeak = funx where funx is defined as

funx <- function (x, y) {
   out <- getparms(x, y) 
   return(out)
}

Note that arguments x and y for funx will be set internally as time/x and the velocity curve returned from the marginaleffects::slopes()

Thank you.

vincentarelbundock commented 11 months ago

Thanks for the suggestion. I am not an epidemiologist and am not familiar with these concepts.

Could you give me a full replicable example with public data? For example, take one of the examples here: https://marginaleffects.com/articles/brms.html

Then, call posteriordraws() to extract the draws. Finally, show me how to manipulate them to obtain the quantities of interest.

Sandhu-SS commented 11 months ago

Thank you for looking into it.

Please see below an example of velocity curve along with the peak velocity (PV) and age at peak velocity (APV).

In this example, I have used the crude approach (index based) to find the PV and the APV. There are many packages that provide much better methods to find the maxima and the minima, such as findpeaksfunction from the pracma R package (findpeaks: Find Peaks in pracma: Practical Numerical Math Functions (rdrr.io)) and the getPeak function from the sitar package (getPeakTrough: Identify peak or trough on curve in sitar: Super Imposition by Translation and Rotation Growth Curve Analysis (rdrr.io)).

 # Load required packages
  suppressPackageStartupMessages({
    library(tidyverse, quietly = TRUE)
    library(marginaleffects, quietly = TRUE)
    library(brms, quietly = TRUE)
    })

  # Get publicly available growth data from the 'fda' R package
  # Data comprise of longitudinal height measurements made on 38 boys between
  # 6 and 18 years of age. Data analysed below are a subset of data.

  # data <- fda::growth$hgtm %>% data.frame() %>% 
  # tibble::rownames_to_column(., "age")

  # The same data is copied below to avoid installing the 'fda' package

  data <- read.table(text = "
    age boy01   boy02   boy03   boy04   boy05   boy06   boy07   boy08   boy09   boy10   boy11   boy12   boy13   boy14   boy15   boy16   boy17   boy18   boy19   boy20   boy21   boy22   boy23   boy24   boy25   boy26   boy27   boy28   boy29   boy30   boy31   boy32   boy33   boy34   boy35   boy36   boy37   boy38   boy39
1   1   81.3    76.2    76.8    74.1    74.2    76.8    72.4    73.8    75.4    78.8    76.9    81.6    78  76.4    76.4    76.2    75  79.7    70  73.8    75.3    76.6    73  75.1    77.4    73.2    78  76.5    78  74.9    74.2    78.7    78  76  76  68.8    82.2    77  75
2   1.25    84.2    80.4    79.8    78.4    76.3    79.1    76  78.7    81  83.3    79.9    83.7    81.8    79.4    81.2    79.2    78.4    81.3    74.5    76.7    79.1    79  76  78.6    81.6    76.2    82.2    80  80.6    78.7    78.3    83.3    82.6    81.2    80.1    73.5    88.9    82.5    79.2
3   1.5 86.4    83.2    82.6    82.6    78.3    81.1    79.4    83  84.9    87  84.1    86.3    85  83.4    86  82.3    82  83.3    78.7    79.3    82.6    82.6    78.7    81.9    84.8    79  86.1    82.8    85.4    81  82.3    87  85.8    84.6    85.9    77.5    90.9    85.6    82.2
4   1.75    88.9    85.4    84.7    85.4    80.3    84.4    82  85.8    87.9    89.6    88.5    88.8    86.4    87.6    89.2    85.4    84  86.5    82.3    82.8    85.1    86.6    82.2    84  88  81.8    88.9    85.4    90.2    84.9    86  90.4    88.4    87  91.4    80.8    93  87.9    84.6
5   2   91.4    87.6    86.7    88.1    82.2    87.4    84.2    88.4    90  91.4    90.6    92.2    87.1    91.4    92.2    88.4    85.9    88.9    85.1    86.4    87.2    89  84.6    85.6    89.6    84.8    91.4    88.4    93.5    87.7    89.1    93.1    90.9    89.1    95.2    84  95.3    90.3    87
6   3   101.1   97  94.2    98.6    89.4    94  93.2    97.3    97.3    100.4   96.6    99.3    96.2    101.2   101.3   101 95.6    99.4    93.6    95.3    95.6    95.3    94.7    92.8    96.5    95.3    101.3   95.2    103.6   93.2    97.1    103 99.1    94.6    103.5   95.2    103.8   98.3    95.8
7   4   109.5   104.6   100.4   104.4   96.9    102.4   101.6   106.7   102.8   111 104.8   105.8   104 106.4   110.3   107.4   102.1   104.1   101.9   102.2   102.6   102.1   104.9   100.6   103.3   101.6   108 102 111 104.2   104.9   109.7   105.8   102 111.6   103.1   111.8   107.2   103.1
8   5   115.8   112.3   107.1   111 104.1   109.2   109 112.8   108.2   118.5   112 112.6   111 113.3   117.3   114.8   109.2   112.5   108.8   109.5   109.7   108.5   109 106.7   108.7   109.3   114.9   107.5   117.9   110 112 118 112.2   107.4   118.1   111.1   117.1   114.6   110.2
9   6   121.9   118.9   112.3   116.3   110.7   116.1   117.3   119.1   113.5   125.8   118.6   119.6   117.3   120.3   122.5   120.7   117.5   119 115.3   116.1   116.8   111.5   115.3   112.5   113.8   116 121.5   114 126.7   116 118.6   124.6   117.6   113 126.3   118.4   123.7   121.9   116.1
10  7   130 125 118.6   123.2   115.8   121.9   122.9   125.7   120.2   131.1   124.8   126.9   124.1   128.8   129.5   128 124.8   128 122 121.7   121.7   119.4   121.2   118.6   120.1   122.7   128.3   119.4   133.6   123.8   124.7   132.3   123.4   119.4   132.8   124.9   130.4   129.2   122.2
11  8   138.2   130.1   124 129.9   121.7   127.8   128.3   132.3   125.7   138.4   130.3   133.8   130 135 134.8   134.5   131.5   133.6   127.9   128.3   127.1   126.1   126.5   123.7   125.8   129.9   134 124.2   140.3   130.6   131.8   138.4   129.6   124.2   139.1   130.7   136.4   135.9   127.7
12  8.5 141.1   133 126.5   133 125.2   129.8   130.9   135.4   127.9   141.5   133.1   136.4   133.6   137.8   137.4   137.4   134.9   136.6   130.9   131.4   129.8   128.8   130.4   126.4   128.9   133.1   137 126.6   143.9   133.3   135 141.6   132.3   127 143.4   134 139.3   139.1   130.7
13  9   144.3   135.4   128.9   136 128.5   132.4   133.7   138.3   130.3   144.5   135.8   139.9   136.8   140.6   140 139.9   137.8   139.5   134.1   134.3   132.5   131.4   134.3   129.1   132.1   136.2   140.1   129.1   147.4   136 138.2   144.8   135 129.8   147.5   137.5   142.2   142.2   133.7
14  9.5 147.5   137.5   131.2   138.7   131.2   135 136 141.3   133.2   147.5   138.7   143.2   139.6   143.4   142.5   142.5   140.7   141.9   136.8   137 134.7   134.3   136.6   131.8   134.3   139.4   142.5   131.6   150.4   138.9   141.2   147.4   137.3   132.9   150.6   140.3   145.2   145 136.8
15  10  150.5   139.7   133.4   141.4   134.2   137.3   138 144.7   135.9   150.5   141.7   146 142.2   146.5   145 145.7   143.8   144.8   139.3   140 136.7   136.8   139.1   134.5   136.9   142.6   145 134 153.4   141.9   143.9   150.5   139.7   136.1   154 143 148.3   147.8   139.8
16  10.5    153.4   142.2   135.8   144 137.4   139.7   139.9   147.9   137.9   153.6   144.2   149 144.8   149.6   147.8   148.7   146.8   148.5   141.6   143 138.9   139.7   141.6   136.8   139.2   144.7   147.4   136.5   156.1   144.1   146.4   153.4   142 139 157.6   145.5   151.2   150.5   142.6
17  11  156.2   144.2   138.4   146.4   140.5   142 141.6   151.2   139.6   156.8   146.6   151.8   147.4   152.2   150.5   151.3   149.5   152.3   144 145.6   141.5   142.2   143.8   139 141.6   147.9   149.5   138.8   159.1   146.2   149.5   156.1   144.4   141.7   160.9   148 154.2   153.5   145.1
18  11.5    159.7   146.2   141 148.8   143.3   144.1   143.4   154.5   141.5   160.5   148.7   154.8   149.8   155.3   153.4   154.3   152.4   156.3   146.1   148.7   143.6   144.1   146.3   141.8   144.2   150.8   151.4   141 162.8   148.1   153.8   159.3   146.5   144.4   163.8   150.8   157.6   157.5   147.3
19  12  163.8   148.1   143.6   151.4   145.9   146.9   145.5   158.2   143.7   164.6   150.5   158 151.9   159.2   157.4   158.1   155.3   161.2   148.6   152.3   145.8   146.6   149 145.3   146.7   153.5   153.2   143.3   167.1   150.1   158.4   163.7   148.6   147.4   166.3   153.6   161.8   161.5   150.2
20  12.5    168.8   149.8   146.8   154.5   148.7   150 147.8   162.3   146 169.2   152.4   161.5   154.1   163.9   162 162.8   158.9   165.9   151.7   157 148 150 151.6   149.7   149.6   156.1   155.2   145.7   171.9   152.5   163.2   168.8   150.7   150.9   168.6   157.3   166.5   165.5   153.3
21  13  174.9   151.6   150.8   157.7   152 153.4   150.3   166.5   148.6   174.3   154 165.3   158.2   168.6   166.6   167.5   163.6   169.5   155 162 150 154 154.2   154.4   153.3   158.7   157.1   148.1   177.1   155.1   168 174.2   153.1   155.9   171 162.1   171.1   169.7   156.7
22  13.5    181.2   153.8   155.1   162.2   155.8   157.3   152.8   171.1   151.5   178.8   158.9   169.6   163.3   172.8   170.7   171.5   168.6   172.1   158.9   166.4   152.5   158 157.2   158.8   157.7   161.1   159.2   150.7   181.8   158.3   172.4   179 156 161.5   173.8   167 174.9   173.1   160.5
23  14  186.3   156.3   159.5   167.4   160.4   161.6   155.8   175.2   154.3   181.9   161.7   174.3   168.1   176.6   174 175.1   173 174.2   163.3   169.8   155.4   162.1   160.6   162.2   162.3   163.7   161.4   154.4   185.7   162.4   175.7   183.1   159.6   166.2   177.3   171.2   177.9   175.6   164.1
24  14.5    189.6   159.2   163.3   172.5   165.1   165.2   159.2   178.4   157.8   184.2   165.1   178.8   173 179.4   176.1   178 176.3   175.6   167.4   172.1   158.9   166.5   164.3   164.8   166.1   166.8   163.6   158.9   188.8   167.3   177.6   186 163.5   169.6   180.7   174.4   180.1   177.2   167.7
25  15  191.3   163.3   166.8   176.3   168.5   168.8   163.1   180.8   161.8   185.6   169.6   182.1   175.2   181.4   177.3   179.8   178.5   176.4   170.4   173.7   163 169.2   168.6   167 168.6   170.3   166.6   162.7   190.8   171.6   178.7   188 167.8   171.9   183.4   176.8   181.3   178.3   171
26  15.5    192.1   167.7   167.8   178.5   170.4   170.1   166.5   182.2   165.6   186.6   174 184.6   177.4   182.6   177.9   181 180 177 171.9   174.7   166.8   170.7   173 168 170.5   174.3   170.9   165.2   192.2   174.9   179.8   189.6   170.3   172.9   185.3   178.4   181.8   179 173.1
27  16  192.8   171.5   168.8   179.8   171.6   171.5   168.9   183.2   168.7   187.1   176.6   186.1   179.2   183.7   178.2   181.9   181 177.3   173 175.5   169.8   171.9   176.2   168.8   172.2   178.3   175 166.9   193 177.4   180.4   190.7   172.8   173.4   186.2   179.1   182.2   179.2   174.7
28  16.5    193.2   174.3   169.8   180.7   172.4   172.1   170.2   184 171 187.4   178.6   187.5   180.6   184.4   178.3   182.7   181.8   177.6   173.7   176.1   172 172.8   178 169.4   173.1   181 178.1   168 193.6   179 180.8   191.3   174.5   173.8   186.4   179.7   182.7   179.3   175.6
29  17  193.8   176.1   170.9   181.4   172.5   172.6   171.1   184.7   171.6   187.8   179.6   188 181.4   185 178.5   183.3   182.2   177.3   174.2   176.4   173.5   173.4   179.1   169.7   173.8   183.1   180.5   168.7   194 180 180.9   191.7   175.5   174.2   186.7   180.1   183.1   179.8   176.1
30  17.5    194.3   177.4   171.2   181.6   172.5   173.2   171.8   185 172.3   188.2   180.2   188.6   181.9   185.7   178.5   183.4   182.6   177 174.6   176.5   174.9   173.8   180.2   169.8   174.4   184.6   182 169.1   194.1   180.7   180.8   192.2   176 174.2   187.2   180.5   183.4   180.3   176.3
31  18  195.1   178.7   171.5   181.8   172.5   173.8   172.6   185.2   172.9   188.4   181.6   189 182.4   185.8   178.7   183.4   183 177 175 176.8   176.4   174 181 170.1   175.2   185.1   182.7   169.4   194.3   181.1   181.3   192.8   176.4   174.1   188 180.8   183.7   180.7   176.4
", header = T, row.names = 1)

  # Convert data to long format and select a subset with age between 6 and 18 years
  data <- data %>% 
    gather(boy, height, -age) %>% 
    mutate(age = as.numeric(age)) %>% 
    mutate(boy = as.factor(boy)) %>% 
    relocate(boy, age, height) %>% 
    filter(age > 6)

  # Data structure
  glimpse(data)
#> Rows: 858
#> Columns: 3
#> $ boy    <fct> boy01, boy01, boy01, boy01, boy01, boy01, boy01, boy01, boy01, …
#> $ age    <dbl> 7.0, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13…
#> $ height <dbl> 130.0, 138.2, 141.1, 144.3, 147.5, 150.5, 153.4, 156.2, 159.7, …

  # Define model 
  # Preece and Baines (1978, Annals of Human Biology) to describe human growth.

  bf <- bf(height ~ h1 - 2 * (h1 - htheta) / (exp(s0 * (age - theta)) + (exp(s1 * (age - theta)))),
           h1     ~ 1,
           htheta ~ 1,
           s0     ~ 1,
           s1     ~ 1,
           theta  ~ 1,
           nl = TRUE)

  # Set priors
  bp = c(
    prior(normal(180, 10),   class = 'b', nlpar = "h1"),
    prior(normal(165, 10),   class = 'b', nlpar = "htheta"),
    prior(cauchy(0.1, 0.1),  class = 'b', nlpar = "s0"),
    prior(cauchy(1.5, 0.5),  class = 'b', nlpar = "s1"),
    prior(normal(14, 2),     class = 'b', nlpar = "theta"))

  # Set initials
  inits <- function(...) {
    list(
      b_h1     = array(rnorm(1, 180, 5),    dim = 1),
      b_htheta = array(rnorm(1, 160, 5),    dim = 1),
      b_s0     = array(rnorm(1, 0.2, 0.05), dim = 1),
      b_s1     = array(rnorm(1, 1.5, 0.25), dim = 1),
      b_theta  = array(rnorm(1, 14, 2),     dim = 1)
    )
  }

  # Fit model
  bfit <- brm(bf, prior = bp, data = data, init = inits, 
              chains = 4, cores = 4, iter = 2000, seed = 123)
#> Compiling Stan program...
#> Start sampling

  # Model fit summary
  summary(bfit)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: height ~ h1 - 2 * (h1 - htheta)/(exp(s0 * (age - theta)) + (exp(s1 * (age - theta)))) 
#>          h1 ~ 1
#>          htheta ~ 1
#>          s0 ~ 1
#>          s1 ~ 1
#>          theta ~ 1
#>    Data: data (Number of observations: 858) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Population-Level Effects: 
#>                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> h1_Intercept       180.43      0.87   178.97   182.32 1.00     1501     1626
#> htheta_Intercept   168.10      0.97   165.94   169.80 1.01     1545     1615
#> s0_Intercept         0.12      0.01     0.09     0.14 1.00     1009     1310
#> s1_Intercept         1.06      0.15     0.78     1.37 1.00     1058     1379
#> theta_Intercept     14.03      0.18    13.65    14.37 1.00     1713     1774
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     6.71      0.17     6.40     7.06 1.00     2422     1925
#> 
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

  # Get slope (i.e., dydx) using the marginaleffects::slopes
  get_slope <- marginaleffects::slopes(bfit)

  # Get velocity and time
  velocity <- get_slope$estimate
  time     <- get_slope$age

  # Crude method to find the peak velocity and the timing by using the index method
  get_index <- which.max(velocity)

  # Peak velocity (PV)
  pv <- velocity[get_index]

  # Age at peak velocity (APV)
  apv <- time[get_index]

  # Plot velocity curve and mark PV (horizontal dotted red line) and APV (vertical dotted red line)

  bfit$data %>% 
    ggplot(., aes(x = age)) +
    geom_line(aes(y = velocity)) +
    geom_hline(aes(yintercept = pv), lty = 2, col = 'red') +
    geom_vline(aes(xintercept = apv), lty = 2, col = 'red')

Created on 2023-10-16 with reprex v2.0.2

In practice, we need to get PV and APV for each posterior draw (i.e., finding PV and APV from the slope at draw level) and them summarizing them across all the draws.

Thanks

vincentarelbundock commented 11 months ago

Is this what you are looking for?

slo <- slopes(bfit, variables = "age")

k = slo |> 
    posterior_draws() |>
    filter(draw == max(draw), .by = "drawid") |>
    select(drawid, PV = draw, APV = age) |>
    unique() |>
    summarize(PV = median(PV), APV = median(APV))

plot_slopes(bfit, variables = "age", condition = "age") +
    geom_hline(aes(yintercept = k$PV), lty = 2, col = "red") +
    geom_vline(aes(xintercept = k$APV), lty = 2, col = "red")

Sandhu-SS commented 11 months ago

Thank you for a neat solution!

It would be great if we can integrate it into the general workflow of the marginaleffects package.

For example, consider an example below where we fit brms model to the data with categorical covariate study (with three levels, study1, study2 and study3) and want to compare these groups for differences in the PV and APV.

Is it possible to create functions similar to the marginaleffects::comparisons() and marginaleffects::avg_comparisons() to compare growth parameters, PV and APV (these can be named as marginaleffects::parms_comparisons () or marginaleffects::avg_parms_comparisons(). Or else, an additional argument parms = TRUE can be added to the already existing comparisons and avg_ comparisons** functions.

Apologies If already an easy solution is available within the marginaleffects framework and missed that.

  # Load required packages
  suppressPackageStartupMessages({
    library(tidyverse, quietly = TRUE)
    library(marginaleffects, quietly = TRUE)
    library(brms, quietly = TRUE)
  })

  # Get publicly available growth data from the 'fda' R package
  # Data comprise of longitudinal height measurements made on 38 boys between
  # 6 and 18 years of age. Data analysed below are a subset of data.

  # data <- fda::growth$hgtm %>% data.frame() %>% 
  # tibble::rownames_to_column(., "age")

  # The same data is copied below to avoid installing the 'fda' package

  data <- read.table(text = "
    age boy01   boy02   boy03   boy04   boy05   boy06   boy07   boy08   boy09   boy10   boy11   boy12   boy13   boy14   boy15   boy16   boy17   boy18   boy19   boy20   boy21   boy22   boy23   boy24   boy25   boy26   boy27   boy28   boy29   boy30   boy31   boy32   boy33   boy34   boy35   boy36   boy37   boy38   boy39
1   1   81.3    76.2    76.8    74.1    74.2    76.8    72.4    73.8    75.4    78.8    76.9    81.6    78  76.4    76.4    76.2    75  79.7    70  73.8    75.3    76.6    73  75.1    77.4    73.2    78  76.5    78  74.9    74.2    78.7    78  76  76  68.8    82.2    77  75
2   1.25    84.2    80.4    79.8    78.4    76.3    79.1    76  78.7    81  83.3    79.9    83.7    81.8    79.4    81.2    79.2    78.4    81.3    74.5    76.7    79.1    79  76  78.6    81.6    76.2    82.2    80  80.6    78.7    78.3    83.3    82.6    81.2    80.1    73.5    88.9    82.5    79.2
3   1.5 86.4    83.2    82.6    82.6    78.3    81.1    79.4    83  84.9    87  84.1    86.3    85  83.4    86  82.3    82  83.3    78.7    79.3    82.6    82.6    78.7    81.9    84.8    79  86.1    82.8    85.4    81  82.3    87  85.8    84.6    85.9    77.5    90.9    85.6    82.2
4   1.75    88.9    85.4    84.7    85.4    80.3    84.4    82  85.8    87.9    89.6    88.5    88.8    86.4    87.6    89.2    85.4    84  86.5    82.3    82.8    85.1    86.6    82.2    84  88  81.8    88.9    85.4    90.2    84.9    86  90.4    88.4    87  91.4    80.8    93  87.9    84.6
5   2   91.4    87.6    86.7    88.1    82.2    87.4    84.2    88.4    90  91.4    90.6    92.2    87.1    91.4    92.2    88.4    85.9    88.9    85.1    86.4    87.2    89  84.6    85.6    89.6    84.8    91.4    88.4    93.5    87.7    89.1    93.1    90.9    89.1    95.2    84  95.3    90.3    87
6   3   101.1   97  94.2    98.6    89.4    94  93.2    97.3    97.3    100.4   96.6    99.3    96.2    101.2   101.3   101 95.6    99.4    93.6    95.3    95.6    95.3    94.7    92.8    96.5    95.3    101.3   95.2    103.6   93.2    97.1    103 99.1    94.6    103.5   95.2    103.8   98.3    95.8
7   4   109.5   104.6   100.4   104.4   96.9    102.4   101.6   106.7   102.8   111 104.8   105.8   104 106.4   110.3   107.4   102.1   104.1   101.9   102.2   102.6   102.1   104.9   100.6   103.3   101.6   108 102 111 104.2   104.9   109.7   105.8   102 111.6   103.1   111.8   107.2   103.1
8   5   115.8   112.3   107.1   111 104.1   109.2   109 112.8   108.2   118.5   112 112.6   111 113.3   117.3   114.8   109.2   112.5   108.8   109.5   109.7   108.5   109 106.7   108.7   109.3   114.9   107.5   117.9   110 112 118 112.2   107.4   118.1   111.1   117.1   114.6   110.2
9   6   121.9   118.9   112.3   116.3   110.7   116.1   117.3   119.1   113.5   125.8   118.6   119.6   117.3   120.3   122.5   120.7   117.5   119 115.3   116.1   116.8   111.5   115.3   112.5   113.8   116 121.5   114 126.7   116 118.6   124.6   117.6   113 126.3   118.4   123.7   121.9   116.1
10  7   130 125 118.6   123.2   115.8   121.9   122.9   125.7   120.2   131.1   124.8   126.9   124.1   128.8   129.5   128 124.8   128 122 121.7   121.7   119.4   121.2   118.6   120.1   122.7   128.3   119.4   133.6   123.8   124.7   132.3   123.4   119.4   132.8   124.9   130.4   129.2   122.2
11  8   138.2   130.1   124 129.9   121.7   127.8   128.3   132.3   125.7   138.4   130.3   133.8   130 135 134.8   134.5   131.5   133.6   127.9   128.3   127.1   126.1   126.5   123.7   125.8   129.9   134 124.2   140.3   130.6   131.8   138.4   129.6   124.2   139.1   130.7   136.4   135.9   127.7
12  8.5 141.1   133 126.5   133 125.2   129.8   130.9   135.4   127.9   141.5   133.1   136.4   133.6   137.8   137.4   137.4   134.9   136.6   130.9   131.4   129.8   128.8   130.4   126.4   128.9   133.1   137 126.6   143.9   133.3   135 141.6   132.3   127 143.4   134 139.3   139.1   130.7
13  9   144.3   135.4   128.9   136 128.5   132.4   133.7   138.3   130.3   144.5   135.8   139.9   136.8   140.6   140 139.9   137.8   139.5   134.1   134.3   132.5   131.4   134.3   129.1   132.1   136.2   140.1   129.1   147.4   136 138.2   144.8   135 129.8   147.5   137.5   142.2   142.2   133.7
14  9.5 147.5   137.5   131.2   138.7   131.2   135 136 141.3   133.2   147.5   138.7   143.2   139.6   143.4   142.5   142.5   140.7   141.9   136.8   137 134.7   134.3   136.6   131.8   134.3   139.4   142.5   131.6   150.4   138.9   141.2   147.4   137.3   132.9   150.6   140.3   145.2   145 136.8
15  10  150.5   139.7   133.4   141.4   134.2   137.3   138 144.7   135.9   150.5   141.7   146 142.2   146.5   145 145.7   143.8   144.8   139.3   140 136.7   136.8   139.1   134.5   136.9   142.6   145 134 153.4   141.9   143.9   150.5   139.7   136.1   154 143 148.3   147.8   139.8
16  10.5    153.4   142.2   135.8   144 137.4   139.7   139.9   147.9   137.9   153.6   144.2   149 144.8   149.6   147.8   148.7   146.8   148.5   141.6   143 138.9   139.7   141.6   136.8   139.2   144.7   147.4   136.5   156.1   144.1   146.4   153.4   142 139 157.6   145.5   151.2   150.5   142.6
17  11  156.2   144.2   138.4   146.4   140.5   142 141.6   151.2   139.6   156.8   146.6   151.8   147.4   152.2   150.5   151.3   149.5   152.3   144 145.6   141.5   142.2   143.8   139 141.6   147.9   149.5   138.8   159.1   146.2   149.5   156.1   144.4   141.7   160.9   148 154.2   153.5   145.1
18  11.5    159.7   146.2   141 148.8   143.3   144.1   143.4   154.5   141.5   160.5   148.7   154.8   149.8   155.3   153.4   154.3   152.4   156.3   146.1   148.7   143.6   144.1   146.3   141.8   144.2   150.8   151.4   141 162.8   148.1   153.8   159.3   146.5   144.4   163.8   150.8   157.6   157.5   147.3
19  12  163.8   148.1   143.6   151.4   145.9   146.9   145.5   158.2   143.7   164.6   150.5   158 151.9   159.2   157.4   158.1   155.3   161.2   148.6   152.3   145.8   146.6   149 145.3   146.7   153.5   153.2   143.3   167.1   150.1   158.4   163.7   148.6   147.4   166.3   153.6   161.8   161.5   150.2
20  12.5    168.8   149.8   146.8   154.5   148.7   150 147.8   162.3   146 169.2   152.4   161.5   154.1   163.9   162 162.8   158.9   165.9   151.7   157 148 150 151.6   149.7   149.6   156.1   155.2   145.7   171.9   152.5   163.2   168.8   150.7   150.9   168.6   157.3   166.5   165.5   153.3
21  13  174.9   151.6   150.8   157.7   152 153.4   150.3   166.5   148.6   174.3   154 165.3   158.2   168.6   166.6   167.5   163.6   169.5   155 162 150 154 154.2   154.4   153.3   158.7   157.1   148.1   177.1   155.1   168 174.2   153.1   155.9   171 162.1   171.1   169.7   156.7
22  13.5    181.2   153.8   155.1   162.2   155.8   157.3   152.8   171.1   151.5   178.8   158.9   169.6   163.3   172.8   170.7   171.5   168.6   172.1   158.9   166.4   152.5   158 157.2   158.8   157.7   161.1   159.2   150.7   181.8   158.3   172.4   179 156 161.5   173.8   167 174.9   173.1   160.5
23  14  186.3   156.3   159.5   167.4   160.4   161.6   155.8   175.2   154.3   181.9   161.7   174.3   168.1   176.6   174 175.1   173 174.2   163.3   169.8   155.4   162.1   160.6   162.2   162.3   163.7   161.4   154.4   185.7   162.4   175.7   183.1   159.6   166.2   177.3   171.2   177.9   175.6   164.1
24  14.5    189.6   159.2   163.3   172.5   165.1   165.2   159.2   178.4   157.8   184.2   165.1   178.8   173 179.4   176.1   178 176.3   175.6   167.4   172.1   158.9   166.5   164.3   164.8   166.1   166.8   163.6   158.9   188.8   167.3   177.6   186 163.5   169.6   180.7   174.4   180.1   177.2   167.7
25  15  191.3   163.3   166.8   176.3   168.5   168.8   163.1   180.8   161.8   185.6   169.6   182.1   175.2   181.4   177.3   179.8   178.5   176.4   170.4   173.7   163 169.2   168.6   167 168.6   170.3   166.6   162.7   190.8   171.6   178.7   188 167.8   171.9   183.4   176.8   181.3   178.3   171
26  15.5    192.1   167.7   167.8   178.5   170.4   170.1   166.5   182.2   165.6   186.6   174 184.6   177.4   182.6   177.9   181 180 177 171.9   174.7   166.8   170.7   173 168 170.5   174.3   170.9   165.2   192.2   174.9   179.8   189.6   170.3   172.9   185.3   178.4   181.8   179 173.1
27  16  192.8   171.5   168.8   179.8   171.6   171.5   168.9   183.2   168.7   187.1   176.6   186.1   179.2   183.7   178.2   181.9   181 177.3   173 175.5   169.8   171.9   176.2   168.8   172.2   178.3   175 166.9   193 177.4   180.4   190.7   172.8   173.4   186.2   179.1   182.2   179.2   174.7
28  16.5    193.2   174.3   169.8   180.7   172.4   172.1   170.2   184 171 187.4   178.6   187.5   180.6   184.4   178.3   182.7   181.8   177.6   173.7   176.1   172 172.8   178 169.4   173.1   181 178.1   168 193.6   179 180.8   191.3   174.5   173.8   186.4   179.7   182.7   179.3   175.6
29  17  193.8   176.1   170.9   181.4   172.5   172.6   171.1   184.7   171.6   187.8   179.6   188 181.4   185 178.5   183.3   182.2   177.3   174.2   176.4   173.5   173.4   179.1   169.7   173.8   183.1   180.5   168.7   194 180 180.9   191.7   175.5   174.2   186.7   180.1   183.1   179.8   176.1
30  17.5    194.3   177.4   171.2   181.6   172.5   173.2   171.8   185 172.3   188.2   180.2   188.6   181.9   185.7   178.5   183.4   182.6   177 174.6   176.5   174.9   173.8   180.2   169.8   174.4   184.6   182 169.1   194.1   180.7   180.8   192.2   176 174.2   187.2   180.5   183.4   180.3   176.3
31  18  195.1   178.7   171.5   181.8   172.5   173.8   172.6   185.2   172.9   188.4   181.6   189 182.4   185.8   178.7   183.4   183 177 175 176.8   176.4   174 181 170.1   175.2   185.1   182.7   169.4   194.3   181.1   181.3   192.8   176.4   174.1   188 180.8   183.7   180.7   176.4
", header = T, row.names = 1)

  # Convert data to long format and select a subset with age between 6 and 18 years
  data <- data %>% 
    gather(boy, height, -age) %>% 
    mutate(age = as.numeric(age)) %>% 
    mutate(boy = as.factor(boy)) %>% 
    relocate(boy, age, height) %>% 
    filter(age > 6)

  # Create a categorical covariate, study 
  set.seed(123)
  data <- data %>% mutate(temp = as.numeric(boy)) %>% 
    mutate(study = cut(temp, breaks=3, labels=c('study1', 'study2', 'study3'))) %>% 
    select(-temp)

  # Data structure
  glimpse(data)
#> Rows: 858
#> Columns: 4
#> $ boy    <fct> boy01, boy01, boy01, boy01, boy01, boy01, boy01, boy01, boy01, …
#> $ age    <dbl> 7.0, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13…
#> $ height <dbl> 130.0, 138.2, 141.1, 144.3, 147.5, 150.5, 153.4, 156.2, 159.7, …
#> $ study  <fct> study1, study1, study1, study1, study1, study1, study1, study1,…

  # Define model with categorical covariate, study 
  # Preece and Baines (1978, Annals of Human Biology) to describe human growth.

  bf <- bf(height ~ h1 - 2 * (h1 - htheta) / (exp(s0 * (age - theta)) + (exp(s1 * (age - theta)))),
           h1     ~ 1 + study,
           htheta ~ 1 + study,
           s0     ~ 1 + study,
           s1     ~ 1 + study,
           theta  ~ 1 + study,
           nl = TRUE)

  # Set priors
  bp = c(
    prior(normal(180, 10),   class = 'b', coef = 'Intercept', nlpar = "h1"),
    prior(normal(165, 10),   class = 'b', coef = 'Intercept', nlpar = "htheta"),
    prior(cauchy(0.1, 0.1),  class = 'b', coef = 'Intercept', nlpar = "s0"),
    prior(cauchy(1.5, 0.5),  class = 'b', coef = 'Intercept', nlpar = "s1"),
    prior(normal(14, 2),     class = 'b', coef = 'Intercept', nlpar = "theta"),

    prior(normal(0, 5),   class = 'b', nlpar = "h1"),
    prior(normal(0, 5),   class = 'b', nlpar = "htheta"),
    prior(normal(0, 0.05),  class = 'b', nlpar = "s0"),
    prior(normal(0, 0.1),  class = 'b', nlpar = "s1"),
    prior(normal(0, 1),     class = 'b', nlpar = "theta")
  )

  # Set initials
  inits <- function(...) {
    list(
      b_h1     = array(c(rnorm(1, 180, 5), 0, 0),    dim = 3),
      b_htheta = array(c(rnorm(1, 160, 5), 0, 0),    dim = 3),
      b_s0     = array(c(rnorm(1, 0.2, 0.05), 0, 0), dim = 3),
      b_s1     = array(c(rnorm(1, 1.5, 0.25), 0, 0), dim = 3),
      b_theta  = array(c(rnorm(1, 14, 2), 0, 0),     dim = 3)
    )
  }

  # Fit model with categorical covariate, study 
  bfit <- brm(bf, prior = bp, data = data, init = inits, 
              threads = threading(8),
              chains = 2, cores = 16, iter = 2000, seed = 123)
#> Compiling Stan program...
#> Start sampling

  # Model fit summary
  summary(bfit)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: height ~ h1 - 2 * (h1 - htheta)/(exp(s0 * (age - theta)) + (exp(s1 * (age - theta)))) 
#>          h1 ~ 1 + study
#>          htheta ~ 1 + study
#>          s0 ~ 1 + study
#>          s1 ~ 1 + study
#>          theta ~ 1 + study
#>    Data: data (Number of observations: 858) 
#>   Draws: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 2000
#> 
#> Population-Level Effects: 
#>                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> h1_Intercept         180.78      1.10   178.86   183.26 1.00     1038     1187
#> h1_studystudy2        -1.89      1.23    -4.21     0.60 1.00     1129     1419
#> h1_studystudy3         1.09      1.22    -1.38     3.56 1.00     1161     1382
#> htheta_Intercept     167.92      1.33   165.02   170.32 1.00      791      921
#> htheta_studystudy2    -1.51      1.68    -4.79     1.83 1.00      895     1049
#> htheta_studystudy3     1.59      1.74    -1.81     4.92 1.00      848     1079
#> s0_Intercept           0.11      0.01     0.08     0.13 1.00      943      957
#> s0_studystudy2         0.01      0.01    -0.02     0.03 1.00     1144     1477
#> s0_studystudy3         0.01      0.01    -0.02     0.04 1.00     1269     1296
#> s1_Intercept           1.05      0.15     0.76     1.35 1.00     1018      932
#> s1_studystudy2        -0.02      0.09    -0.20     0.16 1.00     1546     1300
#> s1_studystudy3        -0.01      0.09    -0.18     0.17 1.00     1609     1369
#> theta_Intercept       14.21      0.24    13.71    14.68 1.00      836     1092
#> theta_studystudy2     -0.27      0.33    -0.92     0.38 1.00      928      984
#> theta_studystudy3     -0.32      0.34    -0.96     0.37 1.00      915     1300
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     6.57      0.16     6.27     6.88 1.00     1839     1510
#> 
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

  # Plot slopes for each level of Study
  plot_slopes(bfit, variables = "age", condition = c('age', 'study'))

Created on 2023-10-16 with reprex v2.0.2

In the above example, we are interested in summarizing the pairwise differences in the PV and APV (along with CIs) between study1, study2, and study3.

vincentarelbundock commented 11 months ago

Yes, that should be pretty easy. The beauty of bayesian models is that you can just extract the raw draws with posterior_draws(), and then manipulate those draws as you want. Just modify the code I gave above to summarize by group, and you should be fine.

I would discourage parms as a name. It's not obvious at all what this would do just by looking at that name.

Frankly, I'm not sure if I'll work on adding this feature to the package. I understand that people use this in epidemiology, but I don't in my own research, so the personal incentives to invest a lot of time developing this are limited.

My sense is that this would be pretty easy to do by "wrapping" functions from marginaleffects into a few separate functions, and post-processing the results. If you want to try your had at that and come up with a full-featured function, then I'd be happy to to consider a few avenues:

  1. Including a chapter on the website to illustrate computation of APV and PV: https://marginaleffects.com
  2. Merge the function into the package (I can't make any promises here, and I'm not yet convinced the use-case is common enough to justify a whole new function, as I try to keep the interface lean).
  3. Develop a separate new package which uses marginaleffects as a dependency to compute key quantities, but adds new functionality.
vincentarelbundock commented 11 months ago

Here’s another fun thing you can do: use custom functions in the comparison argument of the comparisons() function. This requires the ipw development branch of marginaleffects from Github:

remotes::install_github("vincentarelbundock/marginaleffects@ipw")
pv <- \(hi, lo) max((hi - lo) / 1e-6)

apv <- \(hi, lo, newdata) {
  dydx <- (hi - lo) / 1e-6
  idx <- which.max(dydx)
  out <- newdata$age[idx]
  return(out)
}

comparisons(bfit, variables = list(age = 1e-6), comparison = pv)
# 
#  Term Contrast Estimate 2.5 % 97.5 %
#   age   +1e-06     7.59  6.87   8.56
# 
# Columns: term, contrast, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx 
# Type:  response

comparisons(bfit, variables = list(age = 1e-6), comparison = apv)
# 
#  Term Contrast Estimate 2.5 % 97.5 %
#   age   +1e-06     13.5     7     14
# 
# Columns: term, contrast, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx 
# Type:  response
Sandhu-SS commented 11 months ago

Thank you for showing another way to get PV and APV estimates.

However, I feel this does not tap into the general (and great!) framework of the marginaleffects.

A great strength of your marginaleffects package is that it offers a great flexibility and ease of obtaining complex estimates such as pairwise comparisons across factor covariates.

I wish and hope you will add this feature to the marginaleffects whenever you get time.

Thank you very much for developing the marginaleffects.

vincentarelbundock commented 11 months ago

If I may, it fells like you might be missing the point of my example. The power of the comparisons() function is that we can combine the different operations.

Consider this simple GAM model:

library(marginaleffects)
library(itsadug)
library(mgcv)

simdat$Subject <- as.factor(simdat$Subject)
model <- bam(Y ~ Group + s(Time, by = Group) + s(Subject, bs = "re"),
             data = simdat)

plot_slopes(model, variables = "Time", condition = c("Time", "Group"))

Minimum velocity, overall:

comparisons(model,
  variables = list("Time" = 1e-6),
  vcov = FALSE,
  comparison = \(hi, lo) min((hi - lo) / 1e-6))
# 
#  Term Contrast Estimate
#  Time   +1e-06  -0.0267
# 
# Columns: term, contrast, estimate, predicted_lo, predicted_hi, predicted 
# Type:  response

Minimum velocity by group:

comparisons(model,
  variables = list("Time" = 1e-6),
  by = "Group",
  vcov = FALSE,
  comparison = \(hi, lo) min((hi - lo) / 1e-6))
# 
#  Term Contrast    Group Estimate
#  Time   +1e-06 Children  -0.0153
#  Time   +1e-06 Adults    -0.0267
# 
# Columns: term, contrast, Group, estimate 
# Type:  response

Difference between the minimum velocity of each group:

comparisons(model,
  variables = list("Time" = 1e-6),
  vcov = FALSE,
  by = "Group",
  hypothesis = "pairwise",
  comparison = \(hi, lo) min((hi - lo) / 1e-6))
# 
#               Term Estimate
#  Children - Adults   0.0114
# 
# Columns: term, estimate 
# Type:  response

I’m not 100% sure, but you may find that some of these operations are not available for Bayesian models. This is because (1) there are some very tricky programming challenges behind the scenes, and (2) it is usually trivial to manipulate the output of posterior_draws() to obtain the quantities of interest anyway, so the return to programming it internally is low.

Sandhu-SS commented 11 months ago

Thank you for your continued support.

It is very difficult to imagine the capabilities and the flexibility of the marginaleffects package. Amazing!

Earlier I failed to notice that we can combine hypothesis, comparison, and by arguments to achieve what we were looking for.

Your approach works flawlessly to the get the min/max of the velocity (dydx) for each group and to perform pairwise companions.

However, to fully extend the capabilities of this approach in getting the group specific estimates and to perform companions across groups, we need group specific time vector.

Below I explain my point by showing an example.

Fit a simple linear model:

set.seed(20)
x <- seq(from=1, to=20, by=0.1)
n <- length(x)
fac <- cut(x, breaks = 2, labels = c('g1', 'g2'))
y <- 10 + 1 * I(x^1) + 1 *  I(x^2) + 1 *  I(x^3) - 0.035 *  I(x^4)
e <- rnorm(n, mean=0, sd=2)
y <- y + e

data <- cbind.data.frame(y, x, fac)
model <- lm(y ~ fac * poly(x, degree = 4, raw = TRUE), data = data)

Plot fitted and velocity (dydx) curves:

# Plot fitted curve 
marginaleffects::plot_predictions(model, condition = c("x", "fac"))

# Plot velocity (dydx) curve 
marginaleffects::plot_slopes(model, variables = "x", condition = c("x", "fac"))

Define a function that will be called from within the marginaleffects::comparisons() via the comparison argument.

growthparameter_function <- function(hi, lo, newdata,...) {
  dydx  <- (hi - lo) / 1e-6
  out <- max(dydx)
  return(out)
}

Maximum velocity: overall, group specific, and differences between groups (as you have already shown earlier)

newdata <- data

# Overall
marginaleffects::comparisons(model, newdata = newdata,
            variables = list('x' = 1e-6),
            comparison = growthparameter_function)

# By group

marginaleffects::comparisons(model, newdata = newdata,
            variables = list('x' = 1e-6),
            by = "fac",
            comparison = growthparameter_function)

# Difference between groups

marginaleffects::comparisons(model, newdata = newdata,
            variables = list('x' = 1e-6),
            by = "fac",
            hypothesis = "pairwise",
            comparison = growthparameter_function)

Now consider a function (peak_fun(), see below) that takes two inputs, time and the dydx that are specified via arguments x and y, respectively. 

peak_fun <- function (x, y) {
  xy <- xy.coords(x, y)
  xy <- unique(as.data.frame(xy[1:2])[order(xy$x), ])
  x <- xy$x
  y <- xy$y
  ddy <- diff(diff(y) > 0)
  tp <- which(ddy == -1) + 1
  tp <- tp[which.max(y[tp])]
  peak = TRUE
  takeoff = FALSE
  if (length(tp) == 0 && (peak || takeoff)) 
    return(c(x = NA, y = NA))
  if (!peak) {
    if (!takeoff) 
      tp <- length(ddy)
    tp <- which(ddy[seq_len(tp)] == 1) + 1
    tp <- tp[which.min(y[tp])]
    if (length(tp) == 0) 
      return(c(x = NA, y = NA))
  }
  n <- 0
  repeat {
    n <- n + 1
    if (tp == n || tp + n > nrow(xy)) 
      break
    curve <- with(xy[(tp - n):(tp + n), ], lm(y ~ poly(x, 
                                                       2, raw = TRUE)))
    if (curve$rank == 3) 
      break
  }
  x <- -curve$coef[[2]]/curve$coef[[3]]/2
  return(x)
}

This peak_fun() function replaces the max() in growthparameter_function(). Let's call this growthparameter_function2().

growthparameter_function2 <- function(hi, lo, newdata,...) {
  dydx  <- (hi - lo) / 1e-6
 x <- newdata$x
  out <- peak_fun(x = x, y = dydx)
  return(out)
}

While this growthparameter_function2() works fine for the overall estimates, it fails when we try to get group specific estimates (i.e., by group) make comparisons (groups differences) because of the dimensions mismatch between the x and y that are passed to the 'peak_fun(). While the length of y (i.e., dydx) is same as the number of rows for the level of factor variable being evaluated, the x (i.e., age/time passed via the newdata$x) is the total number of rows.

# Overall - works fine
marginaleffects::comparisons(model, newdata = newdata,
            variables = list('x' = 1e-6),
            comparison = growthparameter_function2)

# By group -  fails

marginaleffects::comparisons(model, newdata = newdata,
            variables = list('x' = 1e-6),
            by = "fac",
            comparison = growthparameter_function2)

# Difference between groups -  fails

marginaleffects::comparisons(model, newdata = newdata,
            variables = list('x' = 1e-6),
            by = "fac",
            hypothesis = "pairwise",
            comparison = growthparameter_function2)
vincentarelbundock commented 11 months ago

This is tricker, but possible in this particular case. The custom function you feed to the comparison argument accepts an x argument. In your case, the name coincides with your variable of interest, but in comparisons(), x is just the generic name for the variable with respect to which we are taking the derivative: "x" in "dy/dx".

I use this internally for elasticities. See here: https://github.com/vincentarelbundock/marginaleffects/blob/main/R/sanitize_comparison.R#L12

Sandhu-SS commented 11 months ago

Hi @vincentarelbundock

Unfortunately, the x passed internally by the comparison and the newdata$x are exactly the same. Therefore, the same problem (dimension mismatch) exists even if we use x.

Ideally, I feel the x should correspond to the group for which dydx is computed.

Please see below the dimensions mismatch between dydx and ``x/newdata$x``` 

# Function to show  dimensions mismatch 

growthparameter_function2 <- function(hi, lo, newdata,..) {
  dydx  <- (hi - lo) / 1e-6
  message('length of dydx         = ', length(dydx))
  message('length of x            = ', length(x))
  message('length of newdata$x    = ', length(newdata$x))
  message('identical x and newdata$x: ',       identical(x, newdata$x))
  out <- peak_fun(x = newdata$x, y = dydx)
  return(out)
}

# Overall - works fine

marginaleffects::comparisons(model, #newdata = newdata,
                             variables = list('x' = 1e-6),
                             comparison = growthparameter_function2)
# By group -  fails

marginaleffects::comparisons(model, newdata = newdata,
                             variables = list('x' = 1e-6),
                             by = "fac",
                             comparison = growthparameter_function2)
vincentarelbundock commented 11 months ago

Thanks. Version 0.15.1.9017 from Github should fix this:

library(marginaleffects)
library(itsadug)
library(mgcv)

simdat$Subject <- as.factor(simdat$Subject)
model <- bam(Y ~ Group + s(Time, by = Group) + s(Subject, bs = "re"),
             data = simdat)

low = function(hi, lo, x) {
    dydx <- (hi - lo) / 1e-6
    dydx_min <- min(dydx)
    x[dydx == dydx_min][1]
}
comparisons(model,
  variables = list("Time" = 1e-6),
  vcov = FALSE,
  by = "Group",
  comparison = low
)
# 
#  Term Contrast    Group Estimate
#  Time   +1e-06 Children     1313
#  Time   +1e-06 Adults       1556
# 
# Columns: term, contrast, Group, estimate 
# Type:  response

BTW, you probably want to bootstrap (with inferences()) or something similar. Not sure the delta method is appropriate for max/min…

vincentarelbundock commented 11 months ago

I added these examples to the website: https://marginaleffects.com/articles/slopes.html#minimum-or-maximum-slope-velocity

Sandhu-SS commented 11 months ago

Works fine!

Will test it with brms models and will report back if find any issues/problems.

Thank you once again for your time and efforts.

Regards,

Sandhu-SS commented 11 months ago

Hi @vincentarelbundock

Sorry to bother you again but the same problem (dimensions mismatch between x and dydx) persists for the brms models when computing group specific estimates.

As you have shown earlier, it works fine with the mgcv example.

I am using the Version 0.15.1.9017 

utils::packageVersion("marginaleffects")
[1] ‘0.15.1.9017’

mgcv example - works fine for group specific estimates

suppressPackageStartupMessages({
  library(itsadug, quietly = TRUE)
  library(marginaleffects, quietly = TRUE)
  library(mgcv, quietly = TRUE)
})

simdat$Subject <- as.factor(simdat$Subject)
model <- bam(Y ~ Group + s(Time, by = Group) + s(Subject, bs = "re"),
             data = simdat)

# Plot slopes for each level of Study
# plot_slopes(model, variables = "Time", condition = c('Time', 'Group'))

# Function to show  dimensions mismatch 

growthparameter_function2 <- function(hi, lo, x,..) {
  dydx  <- (hi - lo) / 1e-6
  # message('length of dydx         = ', length(dydx))
  # message('length of x            = ', length(x))
  if(length(dydx) == length(x)) {
    out <- dydx
  } else {
    stop()
  }
  return(out)
}

# Overall - works fine

marginaleffects::comparisons(model, 
                             variables = list('Time' = 1e-6),
                             comparison = growthparameter_function2)
#> 
#>  Term Contrast Estimate Std. Error        z Pr(>|z|)   S   2.5 % 97.5 %
#>  Time   +1e-06  0.02612     0.0374  0.69889    0.485 1.0 -0.0471 0.0994
#>  Time   +1e-06  0.02611     0.0216  1.20880    0.227 2.1 -0.0162 0.0684
#>  Time   +1e-06  0.02608     0.1569  0.16622    0.868 0.2 -0.2814 0.3335
#>  Time   +1e-06  0.02601     0.1395  0.18645    0.852 0.2 -0.2474 0.2994
#>  Time   +1e-06  0.02589     0.1361  0.19025    0.849 0.2 -0.2408 0.2926
#> --- 75590 rows omitted. See ?avg_comparisons and ?print.marginaleffects --- 
#>  Time   +1e-06 -0.00613     0.2460 -0.02493    0.980 0.0 -0.4884 0.4761
#>  Time   +1e-06 -0.00597     2.0595 -0.00290    0.998 0.0 -4.0424 4.0305
#>  Time   +1e-06 -0.00587     2.1088 -0.00278    0.998 0.0 -4.1391 4.1274
#>  Time   +1e-06 -0.00582     2.1314 -0.00273    0.998 0.0 -4.1834 4.1717
#>  Time   +1e-06 -0.00581     0.0699 -0.08313    0.934 0.1 -0.1428 0.1312
#> Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted, Y, Group, Time, Subject 
#> Type:  response

# By group -  works fine

marginaleffects::comparisons(model, 
                             variables = list('Time' = 1e-6),
                             by = "Group",
                             comparison = growthparameter_function2)
#> 
#>  Term Contrast    Group Estimate Std. Error      z Pr(>|z|)   S  2.5 % 97.5 %
#>  Time   +1e-06 Children -0.00557     0.0140 -0.398     0.69 0.5 -0.033 0.0218
#>  Time   +1e-06 Adults   -0.00553     0.0314 -0.176     0.86 0.2 -0.067 0.0559
#> 
#> Columns: term, contrast, Group, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
#> Type:  response

utils::packageVersion("marginaleffects")
#> [1] '0.15.1.9017'

Created on 2023-10-20 with reprex v2.0.2

brms example - fails for group specific estimates

  # Load required packages
  suppressPackageStartupMessages({
    library(tidyverse, quietly = TRUE)
    library(marginaleffects, quietly = TRUE)
    library(brms, quietly = TRUE)
  })

  data <- read.table(text = "
    age boy01   boy02   boy03   boy04   boy05   boy06   boy07   boy08   boy09   boy10   boy11   boy12   boy13   boy14   boy15   boy16   boy17   boy18   boy19   boy20   boy21   boy22   boy23   boy24   boy25   boy26   boy27   boy28   boy29   boy30   boy31   boy32   boy33   boy34   boy35   boy36   boy37   boy38   boy39
1   1   81.3    76.2    76.8    74.1    74.2    76.8    72.4    73.8    75.4    78.8    76.9    81.6    78  76.4    76.4    76.2    75  79.7    70  73.8    75.3    76.6    73  75.1    77.4    73.2    78  76.5    78  74.9    74.2    78.7    78  76  76  68.8    82.2    77  75
2   1.25    84.2    80.4    79.8    78.4    76.3    79.1    76  78.7    81  83.3    79.9    83.7    81.8    79.4    81.2    79.2    78.4    81.3    74.5    76.7    79.1    79  76  78.6    81.6    76.2    82.2    80  80.6    78.7    78.3    83.3    82.6    81.2    80.1    73.5    88.9    82.5    79.2
3   1.5 86.4    83.2    82.6    82.6    78.3    81.1    79.4    83  84.9    87  84.1    86.3    85  83.4    86  82.3    82  83.3    78.7    79.3    82.6    82.6    78.7    81.9    84.8    79  86.1    82.8    85.4    81  82.3    87  85.8    84.6    85.9    77.5    90.9    85.6    82.2
4   1.75    88.9    85.4    84.7    85.4    80.3    84.4    82  85.8    87.9    89.6    88.5    88.8    86.4    87.6    89.2    85.4    84  86.5    82.3    82.8    85.1    86.6    82.2    84  88  81.8    88.9    85.4    90.2    84.9    86  90.4    88.4    87  91.4    80.8    93  87.9    84.6
5   2   91.4    87.6    86.7    88.1    82.2    87.4    84.2    88.4    90  91.4    90.6    92.2    87.1    91.4    92.2    88.4    85.9    88.9    85.1    86.4    87.2    89  84.6    85.6    89.6    84.8    91.4    88.4    93.5    87.7    89.1    93.1    90.9    89.1    95.2    84  95.3    90.3    87
6   3   101.1   97  94.2    98.6    89.4    94  93.2    97.3    97.3    100.4   96.6    99.3    96.2    101.2   101.3   101 95.6    99.4    93.6    95.3    95.6    95.3    94.7    92.8    96.5    95.3    101.3   95.2    103.6   93.2    97.1    103 99.1    94.6    103.5   95.2    103.8   98.3    95.8
7   4   109.5   104.6   100.4   104.4   96.9    102.4   101.6   106.7   102.8   111 104.8   105.8   104 106.4   110.3   107.4   102.1   104.1   101.9   102.2   102.6   102.1   104.9   100.6   103.3   101.6   108 102 111 104.2   104.9   109.7   105.8   102 111.6   103.1   111.8   107.2   103.1
8   5   115.8   112.3   107.1   111 104.1   109.2   109 112.8   108.2   118.5   112 112.6   111 113.3   117.3   114.8   109.2   112.5   108.8   109.5   109.7   108.5   109 106.7   108.7   109.3   114.9   107.5   117.9   110 112 118 112.2   107.4   118.1   111.1   117.1   114.6   110.2
9   6   121.9   118.9   112.3   116.3   110.7   116.1   117.3   119.1   113.5   125.8   118.6   119.6   117.3   120.3   122.5   120.7   117.5   119 115.3   116.1   116.8   111.5   115.3   112.5   113.8   116 121.5   114 126.7   116 118.6   124.6   117.6   113 126.3   118.4   123.7   121.9   116.1
10  7   130 125 118.6   123.2   115.8   121.9   122.9   125.7   120.2   131.1   124.8   126.9   124.1   128.8   129.5   128 124.8   128 122 121.7   121.7   119.4   121.2   118.6   120.1   122.7   128.3   119.4   133.6   123.8   124.7   132.3   123.4   119.4   132.8   124.9   130.4   129.2   122.2
11  8   138.2   130.1   124 129.9   121.7   127.8   128.3   132.3   125.7   138.4   130.3   133.8   130 135 134.8   134.5   131.5   133.6   127.9   128.3   127.1   126.1   126.5   123.7   125.8   129.9   134 124.2   140.3   130.6   131.8   138.4   129.6   124.2   139.1   130.7   136.4   135.9   127.7
12  8.5 141.1   133 126.5   133 125.2   129.8   130.9   135.4   127.9   141.5   133.1   136.4   133.6   137.8   137.4   137.4   134.9   136.6   130.9   131.4   129.8   128.8   130.4   126.4   128.9   133.1   137 126.6   143.9   133.3   135 141.6   132.3   127 143.4   134 139.3   139.1   130.7
13  9   144.3   135.4   128.9   136 128.5   132.4   133.7   138.3   130.3   144.5   135.8   139.9   136.8   140.6   140 139.9   137.8   139.5   134.1   134.3   132.5   131.4   134.3   129.1   132.1   136.2   140.1   129.1   147.4   136 138.2   144.8   135 129.8   147.5   137.5   142.2   142.2   133.7
14  9.5 147.5   137.5   131.2   138.7   131.2   135 136 141.3   133.2   147.5   138.7   143.2   139.6   143.4   142.5   142.5   140.7   141.9   136.8   137 134.7   134.3   136.6   131.8   134.3   139.4   142.5   131.6   150.4   138.9   141.2   147.4   137.3   132.9   150.6   140.3   145.2   145 136.8
15  10  150.5   139.7   133.4   141.4   134.2   137.3   138 144.7   135.9   150.5   141.7   146 142.2   146.5   145 145.7   143.8   144.8   139.3   140 136.7   136.8   139.1   134.5   136.9   142.6   145 134 153.4   141.9   143.9   150.5   139.7   136.1   154 143 148.3   147.8   139.8
16  10.5    153.4   142.2   135.8   144 137.4   139.7   139.9   147.9   137.9   153.6   144.2   149 144.8   149.6   147.8   148.7   146.8   148.5   141.6   143 138.9   139.7   141.6   136.8   139.2   144.7   147.4   136.5   156.1   144.1   146.4   153.4   142 139 157.6   145.5   151.2   150.5   142.6
17  11  156.2   144.2   138.4   146.4   140.5   142 141.6   151.2   139.6   156.8   146.6   151.8   147.4   152.2   150.5   151.3   149.5   152.3   144 145.6   141.5   142.2   143.8   139 141.6   147.9   149.5   138.8   159.1   146.2   149.5   156.1   144.4   141.7   160.9   148 154.2   153.5   145.1
18  11.5    159.7   146.2   141 148.8   143.3   144.1   143.4   154.5   141.5   160.5   148.7   154.8   149.8   155.3   153.4   154.3   152.4   156.3   146.1   148.7   143.6   144.1   146.3   141.8   144.2   150.8   151.4   141 162.8   148.1   153.8   159.3   146.5   144.4   163.8   150.8   157.6   157.5   147.3
19  12  163.8   148.1   143.6   151.4   145.9   146.9   145.5   158.2   143.7   164.6   150.5   158 151.9   159.2   157.4   158.1   155.3   161.2   148.6   152.3   145.8   146.6   149 145.3   146.7   153.5   153.2   143.3   167.1   150.1   158.4   163.7   148.6   147.4   166.3   153.6   161.8   161.5   150.2
20  12.5    168.8   149.8   146.8   154.5   148.7   150 147.8   162.3   146 169.2   152.4   161.5   154.1   163.9   162 162.8   158.9   165.9   151.7   157 148 150 151.6   149.7   149.6   156.1   155.2   145.7   171.9   152.5   163.2   168.8   150.7   150.9   168.6   157.3   166.5   165.5   153.3
21  13  174.9   151.6   150.8   157.7   152 153.4   150.3   166.5   148.6   174.3   154 165.3   158.2   168.6   166.6   167.5   163.6   169.5   155 162 150 154 154.2   154.4   153.3   158.7   157.1   148.1   177.1   155.1   168 174.2   153.1   155.9   171 162.1   171.1   169.7   156.7
22  13.5    181.2   153.8   155.1   162.2   155.8   157.3   152.8   171.1   151.5   178.8   158.9   169.6   163.3   172.8   170.7   171.5   168.6   172.1   158.9   166.4   152.5   158 157.2   158.8   157.7   161.1   159.2   150.7   181.8   158.3   172.4   179 156 161.5   173.8   167 174.9   173.1   160.5
23  14  186.3   156.3   159.5   167.4   160.4   161.6   155.8   175.2   154.3   181.9   161.7   174.3   168.1   176.6   174 175.1   173 174.2   163.3   169.8   155.4   162.1   160.6   162.2   162.3   163.7   161.4   154.4   185.7   162.4   175.7   183.1   159.6   166.2   177.3   171.2   177.9   175.6   164.1
24  14.5    189.6   159.2   163.3   172.5   165.1   165.2   159.2   178.4   157.8   184.2   165.1   178.8   173 179.4   176.1   178 176.3   175.6   167.4   172.1   158.9   166.5   164.3   164.8   166.1   166.8   163.6   158.9   188.8   167.3   177.6   186 163.5   169.6   180.7   174.4   180.1   177.2   167.7
25  15  191.3   163.3   166.8   176.3   168.5   168.8   163.1   180.8   161.8   185.6   169.6   182.1   175.2   181.4   177.3   179.8   178.5   176.4   170.4   173.7   163 169.2   168.6   167 168.6   170.3   166.6   162.7   190.8   171.6   178.7   188 167.8   171.9   183.4   176.8   181.3   178.3   171
26  15.5    192.1   167.7   167.8   178.5   170.4   170.1   166.5   182.2   165.6   186.6   174 184.6   177.4   182.6   177.9   181 180 177 171.9   174.7   166.8   170.7   173 168 170.5   174.3   170.9   165.2   192.2   174.9   179.8   189.6   170.3   172.9   185.3   178.4   181.8   179 173.1
27  16  192.8   171.5   168.8   179.8   171.6   171.5   168.9   183.2   168.7   187.1   176.6   186.1   179.2   183.7   178.2   181.9   181 177.3   173 175.5   169.8   171.9   176.2   168.8   172.2   178.3   175 166.9   193 177.4   180.4   190.7   172.8   173.4   186.2   179.1   182.2   179.2   174.7
28  16.5    193.2   174.3   169.8   180.7   172.4   172.1   170.2   184 171 187.4   178.6   187.5   180.6   184.4   178.3   182.7   181.8   177.6   173.7   176.1   172 172.8   178 169.4   173.1   181 178.1   168 193.6   179 180.8   191.3   174.5   173.8   186.4   179.7   182.7   179.3   175.6
29  17  193.8   176.1   170.9   181.4   172.5   172.6   171.1   184.7   171.6   187.8   179.6   188 181.4   185 178.5   183.3   182.2   177.3   174.2   176.4   173.5   173.4   179.1   169.7   173.8   183.1   180.5   168.7   194 180 180.9   191.7   175.5   174.2   186.7   180.1   183.1   179.8   176.1
30  17.5    194.3   177.4   171.2   181.6   172.5   173.2   171.8   185 172.3   188.2   180.2   188.6   181.9   185.7   178.5   183.4   182.6   177 174.6   176.5   174.9   173.8   180.2   169.8   174.4   184.6   182 169.1   194.1   180.7   180.8   192.2   176 174.2   187.2   180.5   183.4   180.3   176.3
31  18  195.1   178.7   171.5   181.8   172.5   173.8   172.6   185.2   172.9   188.4   181.6   189 182.4   185.8   178.7   183.4   183 177 175 176.8   176.4   174 181 170.1   175.2   185.1   182.7   169.4   194.3   181.1   181.3   192.8   176.4   174.1   188 180.8   183.7   180.7   176.4
", header = T, row.names = 1)

  # Convert data to long format and select a subset with age between 6 and 18 years
  data <- data %>% 
    gather(boy, height, -age) %>% 
    mutate(age = as.numeric(age)) %>% 
    mutate(boy = as.factor(boy)) %>% 
    relocate(boy, age, height) %>% 
    filter(age > 6)

  # Create a categorical covariate, study 
  set.seed(123)
  data <- data %>% mutate(temp = as.numeric(boy)) %>% 
    mutate(study = cut(temp, breaks=3, labels=c('study1', 'study2', 'study3'))) %>% 
    dplyr::select(-temp)

  bf <- bf(height ~ h1 - 2 * (h1 - htheta) / (exp(s0 * (age - theta)) + (exp(s1 * (age - theta)))),
           h1     ~ 1 + study,
           htheta ~ 1 + study,
           s0     ~ 1 + study,
           s1     ~ 1 + study,
           theta  ~ 1 + study,
           nl = TRUE)

  # Set priors
  bp = c(
    prior(normal(180, 10),   class = 'b', coef = 'Intercept', nlpar = "h1"),
    prior(normal(165, 10),   class = 'b', coef = 'Intercept', nlpar = "htheta"),
    prior(cauchy(0.1, 0.1),  class = 'b', coef = 'Intercept', nlpar = "s0"),
    prior(cauchy(1.5, 0.5),  class = 'b', coef = 'Intercept', nlpar = "s1"),
    prior(normal(14, 2),     class = 'b', coef = 'Intercept', nlpar = "theta"),

    prior(normal(0, 5),   class = 'b', nlpar = "h1"),
    prior(normal(0, 5),   class = 'b', nlpar = "htheta"),
    prior(normal(0, 0.05),  class = 'b', nlpar = "s0"),
    prior(normal(0, 0.1),  class = 'b', nlpar = "s1"),
    prior(normal(0, 1),     class = 'b', nlpar = "theta")
  )

  # Set initials
  inits <- function(...) {
    list(
      b_h1     = array(c(rnorm(1, 180, 5), 0, 0),    dim = 3),
      b_htheta = array(c(rnorm(1, 160, 5), 0, 0),    dim = 3),
      b_s0     = array(c(rnorm(1, 0.2, 0.05), 0, 0), dim = 3),
      b_s1     = array(c(rnorm(1, 1.5, 0.25), 0, 0), dim = 3),
      b_theta  = array(c(rnorm(1, 14, 2), 0, 0),     dim = 3)
    )
  }

  # Fit model with categorical covariate, study 
  bfit <- brm(bf, prior = bp, data = data, init = inits, 
              chains = 2, cores = 16, iter = 1000, seed = 123)
#> Compiling Stan program...
#> Start sampling

  # Model fit summary
  summary(bfit)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: height ~ h1 - 2 * (h1 - htheta)/(exp(s0 * (age - theta)) + (exp(s1 * (age - theta)))) 
#>          h1 ~ 1 + study
#>          htheta ~ 1 + study
#>          s0 ~ 1 + study
#>          s1 ~ 1 + study
#>          theta ~ 1 + study
#>    Data: data (Number of observations: 858) 
#>   Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
#>          total post-warmup draws = 1000
#> 
#> Population-Level Effects: 
#>                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> h1_Intercept         180.81      1.14   178.87   183.34 1.00      539      569
#> h1_studystudy2        -1.87      1.29    -4.48     0.53 1.00      630      614
#> h1_studystudy3         1.12      1.26    -1.38     3.60 1.00      718      690
#> htheta_Intercept     167.82      1.29   165.25   170.25 1.00      461      591
#> htheta_studystudy2    -1.44      1.69    -4.78     1.61 1.00      511      544
#> htheta_studystudy3     1.65      1.63    -1.60     4.79 1.00      430      615
#> s0_Intercept           0.11      0.01     0.08     0.13 1.00      498      371
#> s0_studystudy2         0.01      0.01    -0.02     0.03 1.00      624      612
#> s0_studystudy3         0.01      0.01    -0.02     0.04 1.00      507      437
#> s1_Intercept           1.04      0.15     0.76     1.34 1.00      521      385
#> s1_studystudy2        -0.02      0.09    -0.20     0.17 1.00      688      539
#> s1_studystudy3        -0.01      0.10    -0.20     0.18 1.00      839      647
#> theta_Intercept       14.19      0.23    13.76    14.65 1.00      456      612
#> theta_studystudy2     -0.26      0.33    -0.92     0.40 1.00      554      536
#> theta_studystudy3     -0.31      0.32    -0.95     0.25 1.00      441      640
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     6.56      0.16     6.26     6.90 1.00     1007      642
#> 
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

  # Plot slopes for each level of Study
  plot_slopes(bfit, variables = "age", condition = c('age', 'study'))


  # Function to show  dimensions mismatch 

  growthparameter_function2 <- function(hi, lo, x,..) {
    dydx  <- (hi - lo) / 1e-6
    # message('length of dydx         = ', length(dydx))
    # message('length of x            = ', length(x))
    if(length(dydx) == length(x)) {
      out <- dydx
    } else {
      stop()
    }
    return(out)
  }

  # Overall - works fine

  marginaleffects::comparisons(bfit, 
                               variables = list('age' = 1e-6),
                               comparison = growthparameter_function2)
#> 
#>  Term Contrast Estimate 2.5 % 97.5 %
#>   age   +1e-06    6.319 4.714  7.709
#>   age   +1e-06    5.748 4.584  6.793
#>   age   +1e-06    5.510 4.570  6.413
#>   age   +1e-06    5.336 4.574  6.064
#>   age   +1e-06    5.205 4.584  5.847
#> --- 848 rows omitted. See ?avg_comparisons and ?print.marginaleffects --- 
#>   age   +1e-06    2.242 1.299  3.158
#>   age   +1e-06    1.470 0.735  2.403
#>   age   +1e-06    0.929 0.401  1.819
#>   age   +1e-06    0.580 0.215  1.353
#>   age   +1e-06    0.353 0.114  0.962
#> Columns: rowid, term, contrast, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx, height, age, study 
#> Type:  response

  # By group -  fails

  marginaleffects::comparisons(bfit,
                               variables = list('age' = 1e-6),
                               by = "study",
                               comparison = growthparameter_function2)
#> Error: The function supplied to the `comparison` argument must accept two
#>   numeric vectors of predicted probabilities of length 286, and return a
#>   single numeric value or a numeric vector of length 286, with no missing
#>   value.

  utils::packageVersion("marginaleffects")
#> [1] '0.15.1.9017'

Created on 2023-10-20 with reprex v2.0.2

vincentarelbundock commented 11 months ago

Yep, I think that's a small bug. Can you try again with version 0.16.0.9000 from github?

Sandhu-SS commented 11 months ago

Hi,

I am not able to install the lates dev version 0.16.0.9000.

The following installs ‘0.16.0’

# Install the development version:

install.packages(
    c("marginaleffects", "insight"),
    repos = c("https://vincentarelbundock.r-universe.dev", "https://easystats.r-universe.dev"))

utils::packageVersion("marginaleffects")

[1] ‘0.16.0’
vincentarelbundock commented 11 months ago

Try 

remotes::install_github("vincentarelbundock/marginaleffects")
Sandhu-SS commented 11 months ago

Awesome!

Works perfectly!

I am really amazed by the speed and the hard work you are putting in resolving the issues quickly.

Thank you so much.

Kind regards

Version ‘0.16.0.9000’

utils::packageVersion("marginaleffects")
[1] ‘0.16.0.9000’

Overall - works fine


marginaleffects::comparisons(bfit, 
                               variables = list('age' = 1e-6),
                               comparison = growthparameter_function2)

Term Contrast Estimate  2.5 %   97.5 %
  age   +1e-06   1.6895 -22.63 687564.6
  age   +1e-06   1.6158 -16.54 199194.2
  age   +1e-06   1.7056 -13.28 107768.8
  age   +1e-06   1.7654 -12.30  57209.9
  age   +1e-06   1.9340 -14.50  30763.0
--- 848 rows omitted. See ?avg_comparisons and ?print.marginaleffects --- 
  age   +1e-06   0.2986  -9.47     35.7
  age   +1e-06   0.1638 -10.41     30.2
  age   +1e-06   0.0815  -9.86     23.1
  age   +1e-06   0.0450  -5.58     19.2
  age   +1e-06   0.0167  -4.05     14.0
Columns: rowid, term, contrast, estimate, conf.low, conf.high, predicted_lo, predicted_hi, predicted, tmp_idx, height, age, study 
Type:  response

By group - works fine

marginaleffects::comparisons(bfit,
                               variables = list('age' = 1e-6),
                               by = "study",
                               comparison = growthparameter_function2)

Term Contrast  study Estimate 2.5 % 97.5 %
  age   +1e-06 study1     2.87 -12.8  52163
  age   +1e-06 study2     2.87 -23.5 119044
  age   +1e-06 study3     2.49 -19.9  52807

Columns: term, contrast, study, estimate, conf.low, conf.high 
Type:  response
vincentarelbundock commented 11 months ago

I appreciate the thumbs up and the detailed reports. Helps for motivation and making the package better.