The emaxnls package provides tools for nonlinear least squares estimation for Emax regression models.
Installation
You can install the development version of emaxnls from GitHub with:
# install.packages("pak")
pak::pak("djnavarro/emaxnls")Minimal example
library(tibble)
library(emaxnls)
set.seed(123)
emax_df
#> # A tibble: 400 × 10
#> dose exp_1 exp_2 rsp_1 rsp_2 cnt_a cnt_b cnt_c bin_d bin_e
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0 0 0 6.49 0 3.06 6.50 5.93 0 1
#> 2 0 0 0 7.81 0 5.72 3.84 5.60 0 1
#> 3 0 0 0 7.26 0 4.31 3.68 8.16 0 0
#> 4 0 0 0 7.45 0 4.03 2.86 9.38 1 1
#> 5 0 0 0 6.33 0 2.46 3.36 8.29 0 1
#> 6 0 0 0 7.13 0 4.87 8.90 7.06 0 0
#> 7 0 0 0 6.07 0 2.87 4.85 0.989 0 1
#> 8 0 0 0 5.47 1 1.07 5.34 3.34 1 1
#> 9 0 0 0 7.12 0 3.94 5.68 4.01 1 1
#> 10 0 0 0 8.21 1 7.46 8.16 6.92 1 1
#> # ℹ 390 more rows
emax_nls(
structural_model = rsp_1 ~ exp_1,
covariate_model = list(E0 ~ cnt_a, Emax ~ 1, logEC50 ~ 1),
data = emax_df
)
#> Structural model:
#>
#> Exposure: exp_1
#> Response: rsp_1
#> Emax type: hyperbolic
#>
#> Covariate model:
#>
#> E0: E0 ~ cnt_a
#> Emax: Emax ~ 1
#> logEC50: logEC50 ~ 1
#>
#> Coefficient table:
#>
#> label estimate std_error t_statistic p_value ci_lower ci_upper
#> 1 E0_Intercept 4.99 0.0740 67.5 3.21e-219 4.85 5.14
#> 2 E0_cnt_a 0.498 0.0113 44.2 4.30e-155 0.476 0.520
#> 3 Emax_Intercept 10.0 0.104 96.3 7.23e-277 9.80 10.2
#> 4 logEC50_Intercept 8.27 0.0366 226. 0 8.19 8.34Stepwise covariate modelling
base_model <- emax_nls(
structural_model = rsp_1 ~ exp_1,
covariate_model = list(E0 ~ 1, Emax ~ 1, logEC50 ~ 1),
data = emax_df
)
covariate_list <- list(
E0 = c("cnt_a", "cnt_b", "cnt_c", "bin_d", "bin_e"),
Emax = c("cnt_a", "cnt_b", "cnt_c", "bin_d", "bin_e")
)
final_mod <- base_model |>
emax_scm_forward(candidates = covariate_list, threshold = .01) |>
emax_scm_backward(candidates = covariate_list, threshold = .001)
emax_scm_history(final_mod)
#> # A tibble: 32 × 11
#> iteration attempt step action term_tested model_tested model_converged
#> <int> <int> <chr> <chr> <chr> <chr> <lgl>
#> 1 0 0 base model <NA> <NA> E0 ~ 1, Ema… TRUE
#> 2 1 1 forward add E0 ~ cnt_c E0 ~ 1 + cn… TRUE
#> 3 1 2 forward add Emax ~ bin_e E0 ~ 1, Ema… TRUE
#> 4 1 3 forward add E0 ~ cnt_b E0 ~ 1 + cn… TRUE
#> 5 1 4 forward add Emax ~ cnt_c E0 ~ 1, Ema… TRUE
#> 6 1 5 forward add Emax ~ cnt_a E0 ~ 1, Ema… TRUE
#> 7 1 6 forward add Emax ~ bin_d E0 ~ 1, Ema… TRUE
#> 8 1 7 forward add E0 ~ cnt_a E0 ~ 1 + cn… TRUE
#> 9 1 8 forward add Emax ~ cnt_b E0 ~ 1, Ema… TRUE
#> 10 1 9 forward add E0 ~ bin_e E0 ~ 1 + bi… TRUE
#> # ℹ 22 more rows
#> # ℹ 4 more variables: term_p_value <dbl>, model_aic <dbl>, model_bic <dbl>,
#> # model_updated <lgl>
final_mod
#> Structural model:
#>
#> Exposure: exp_1
#> Response: rsp_1
#> Emax type: hyperbolic
#>
#> Covariate model:
#>
#> E0: E0 ~ 1 + cnt_a
#> Emax: Emax ~ 1
#> logEC50: logEC50 ~ 1
#>
#> Coefficient table:
#>
#> label estimate std_error t_statistic p_value ci_lower ci_upper
#> 1 E0_Intercept 4.99 0.0740 67.5 3.21e-219 4.85 5.14
#> 2 E0_cnt_a 0.498 0.0113 44.2 4.30e-155 0.476 0.520
#> 3 Emax_Intercept 10.0 0.104 96.3 7.23e-277 9.80 10.2
#> 4 logEC50_Intercept 8.27 0.0366 226. 0 8.19 8.34Simulation
simulate(final_mod, nsim = 1)
#> # A tibble: 400 × 8
#> dat_id sim_id mu val E0_Intercept E0_cnt_a Emax_Intercept
#> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1 6.52 6.86 4.96 0.508 10.1
#> 2 2 1 7.84 8.12 4.96 0.508 10.1
#> 3 3 1 7.14 7.11 4.96 0.508 10.1
#> 4 4 1 7.00 6.85 4.96 0.508 10.1
#> 5 5 1 6.22 6.03 4.96 0.508 10.1
#> 6 6 1 7.42 7.07 4.96 0.508 10.1
#> 7 7 1 6.42 6.32 4.96 0.508 10.1
#> 8 8 1 5.53 4.89 4.96 0.508 10.1
#> 9 9 1 6.95 8.04 4.96 0.508 10.1
#> 10 10 1 8.71 9.31 4.96 0.508 10.1
#> # ℹ 390 more rows
#> # ℹ 1 more variable: logEC50_Intercept <dbl>