A Workflow for Variance-Aware Michaelis-Menten Analysis

inferMM targets a common enzyme-kinetic workflow: fit a Michaelis-Menten curve, compare plausible working variance functions, and summarize inference in a way that is less sensitive to heteroscedasticity than ordinary nonlinear least squares alone. The same interface also supports clustered assay layouts in which repeated measurements within a core, plate, or sample may be correlated.

Single-curve fitting

library(inferMM)

one_curve <- subset(sdl_demo, enzyme == "1111")

fit <- fit_mm(
  x = one_curve$s_uM,
  y = one_curve$v_uM_per_min,
  variance = "sqrt"
)

summary(fit)
#> Call:
#> fit_mm(x = one_curve$s_uM, y = one_curve$v_uM_per_min, variance = "sqrt")
#> 
#> Method: profile-score 
#> Working variance: Var = gamma S^(1/2) 
#> 
#> CI method: wald
#> 
#>      Estimate Std. Error Lower CI Upper CI
#> Vmax  6.43188    0.08808  6.25926   6.6045
#> Km   76.42907    2.36904 71.78583  81.0723
#> 
#> Scale:0.0012  RMSE:0.11978  Weighted RSS:28
#> Quasi-AIC:-44.87655  Quasi-BIC:-40.87994

fit_auto <- fit_mm(
  x = one_curve$s_uM,
  y = one_curve$v_uM_per_min,
  variance = "auto",
  power_selection = "quasi_aic"
)

summary(fit_auto)
#> Call:
#> fit_mm(x = one_curve$s_uM, y = one_curve$v_uM_per_min, variance = "auto", 
#>     power_selection = "quasi_aic")
#> 
#> Method: profile-score 
#> Working variance: Var = gamma S^0.85 (selected by quasi-AIC from log/power candidates) 
#> 
#> CI method: wald
#> 
#>      Estimate Std. Error Lower CI Upper CI
#> Vmax  6.48388    0.10439  6.27927  6.68849
#> Km   77.97381    2.36193 73.34452 82.60310
#> 
#> Scale:0.00026  RMSE:0.12129  Weighted RSS:28
#> Quasi-AIC:-46.89897  Quasi-BIC:-42.90235

The fitted object supports extraction of coefficients, Wald confidence intervals, bootstrap confidence intervals, and predictions.

coef(fit)
#>     Vmax       Km 
#>  6.43188 76.42907
confint(fit)
#>          lower     upper
#> Vmax  6.259256  6.604505
#> Km   71.785829 81.072301
set.seed(1)
confint(fit, method = "bootstrap", B = 99)
#>          lower    upper
#> Vmax  6.224301  6.63876
#> Km   71.079881 80.97992
#> attr(,"bootstrap_draws")
#>            Vmax       Km
#>   [1,] 6.476014 79.83447
#>   [2,] 6.378489 74.52054
#>   [3,] 6.315636 73.40745
#>   [4,] 6.606759 81.49450
#>   [5,] 6.355270 72.46983
#>   [6,] 6.490773 78.61711
#>   [7,] 6.375798 76.89579
#>   [8,] 6.415435 77.91821
#>   [9,] 6.395333 75.79797
#>  [10,] 6.452104 75.77422
#>  [11,] 6.354902 73.64769
#>  [12,] 6.253890 72.46448
#>  [13,] 6.432755 76.48557
#>  [14,] 6.239816 72.73521
#>  [15,] 6.392421 72.35687
#>  [16,] 6.438676 75.82950
#>  [17,] 6.407854 76.57698
#>  [18,] 6.418602 74.41747
#>  [19,] 6.416149 78.11960
#>  [20,] 6.398709 76.94679
#>  [21,] 6.409231 78.06683
#>  [22,] 6.414848 73.83375
#>  [23,] 6.596570 80.87162
#>  [24,] 6.445099 78.07890
#>  [25,] 6.566444 78.84445
#>  [26,] 6.316287 74.66422
#>  [27,] 6.393557 76.88738
#>  [28,] 6.458692 75.04789
#>  [29,] 6.379901 76.13779
#>  [30,] 6.257576 71.96651
#>  [31,] 6.422571 75.72082
#>  [32,] 6.464795 76.89847
#>  [33,] 6.286777 73.34994
#>  [34,] 6.471574 77.19526
#>  [35,] 6.449847 76.42107
#>  [36,] 6.371707 74.37569
#>  [37,] 6.357429 74.21841
#>  [38,] 6.357334 74.64551
#>  [39,] 6.531732 77.41515
#>  [40,] 6.372838 75.80516
#>  [41,] 6.640540 81.91559
#>  [42,] 6.427200 77.41146
#>  [43,] 6.243263 73.29175
#>  [44,] 6.580359 78.40103
#>  [45,] 6.359400 73.67905
#>  [46,] 6.546413 78.07210
#>  [47,] 6.494002 77.34733
#>  [48,] 6.368514 74.86933
#>  [49,] 6.431841 75.52364
#>  [50,] 6.451744 76.90597
#>  [51,] 6.453025 77.24754
#>  [52,] 6.484998 78.84887
#>  [53,] 6.444248 76.11355
#>  [54,] 6.416209 75.24676
#>  [55,] 6.458626 76.53895
#>  [56,] 6.304262 72.37091
#>  [57,] 6.390457 75.06245
#>  [58,] 6.369015 75.27856
#>  [59,] 6.411905 76.62328
#>  [60,] 6.335236 72.76626
#>  [61,] 6.449384 76.19105
#>  [62,] 6.350371 74.00297
#>  [63,] 6.389128 75.40213
#>  [64,] 6.509621 78.61138
#>  [65,] 6.498540 77.34006
#>  [66,] 6.297257 74.08351
#>  [67,] 6.304905 73.09759
#>  [68,] 6.650319 82.14755
#>  [69,] 6.259182 72.65697
#>  [70,] 6.620325 80.46960
#>  [71,] 6.444258 77.05716
#>  [72,] 6.618748 80.11028
#>  [73,] 6.430931 77.26469
#>  [74,] 6.651050 82.28822
#>  [75,] 6.398579 76.10857
#>  [76,] 6.531708 76.63966
#>  [77,] 6.470498 76.55241
#>  [78,] 6.501117 78.43636
#>  [79,] 6.485008 77.65040
#>  [80,] 6.485931 77.24606
#>  [81,] 6.415493 75.51380
#>  [82,] 6.443932 74.05108
#>  [83,] 6.325003 73.76095
#>  [84,] 6.323352 74.28236
#>  [85,] 6.537225 78.46154
#>  [86,] 6.425373 77.41009
#>  [87,] 6.631429 81.96381
#>  [88,] 6.380613 75.35240
#>  [89,] 6.483115 77.88982
#>  [90,] 6.338270 74.75712
#>  [91,] 6.560837 76.39026
#>  [92,] 6.553474 80.97582
#>  [93,] 6.393651 74.96486
#>  [94,] 6.500069 77.84719
#>  [95,] 6.469865 77.91693
#>  [96,] 6.341839 75.80376
#>  [97,] 6.459499 78.36327
#>  [98,] 6.358715 73.68685
#>  [99,] 6.322227 72.95035
#> attr(,"bootstrap_draws")attr(,"requested_B")
#> [1] 99
#> attr(,"bootstrap_draws")attr(,"successful_B")
#> [1] 99
#> attr(,"bootstrap_draws")attr(,"bootstrap_type")
#> [1] "wild"
#> attr(,"bootstrap_draws")attr(,"wild_weights")
#> [1] "mammen"
#> attr(,"bootstrap_draws")attr(,"se_draws")
#>              Vmax       Km
#>   [1,] 0.06985208 1.935124
#>   [2,] 0.10246839 2.721211
#>   [3,] 0.06735485 1.783951
#>   [4,] 0.08946312 2.471365
#>   [5,] 0.08391186 2.185031
#>   [6,] 0.08787750 2.398055
#>   [7,] 0.08215916 2.240689
#>   [8,] 0.08112749 2.223237
#>   [9,] 0.06964822 1.871104
#>  [10,] 0.06270605 1.669348
#>  [11,] 0.09723531 2.566446
#>  [12,] 0.07315613 1.935708
#>  [13,] 0.06387095 1.718819
#>  [14,] 0.06907691 1.837610
#>  [15,] 0.08942662 2.312088
#>  [16,] 0.06834826 1.824457
#>  [17,] 0.06262306 1.693470
#>  [18,] 0.09244142 2.436773
#>  [19,] 0.09928960 2.726513
#>  [20,] 0.07677103 2.087400
#>  [21,] 0.07698530 2.115122
#>  [22,] 0.09008626 2.360508
#>  [23,] 0.08547042 2.349661
#>  [24,] 0.07809181 2.133862
#>  [25,] 0.07836180 2.118843
#>  [26,] 0.07504769 2.015870
#>  [27,] 0.06357009 1.728746
#>  [28,] 0.11550177 3.047138
#>  [29,] 0.07788269 2.105225
#>  [30,] 0.09938417 2.613052
#>  [31,] 0.05738243 1.533743
#>  [32,] 0.08069294 2.170478
#>  [33,] 0.06924554 1.841238
#>  [34,] 0.12740395 3.434337
#>  [35,] 0.07649214 2.051568
#>  [36,] 0.06966066 1.848912
#>  [37,] 0.07615177 2.022158
#>  [38,] 0.06883002 1.836538
#>  [39,] 0.11082711 2.967009
#>  [40,] 0.06943604 1.872134
#>  [41,] 0.08286175 2.287154
#>  [42,] 0.06895697 1.876023
#>  [43,] 0.08397688 2.247013
#>  [44,] 0.11008318 2.956345
#>  [45,] 0.10509946 2.773038
#>  [46,] 0.11417859 3.071443
#>  [47,] 0.06588085 1.772678
#>  [48,] 0.07496306 2.001667
#>  [49,] 0.05970155 1.589968
#>  [50,] 0.10320772 2.781921
#>  [51,] 0.06368762 1.722687
#>  [52,] 0.06434465 1.761755
#>  [53,] 0.06786972 1.815773
#>  [54,] 0.10341771 2.752472
#>  [55,] 0.07591162 2.035847
#>  [56,] 0.10580533 2.774245
#>  [57,] 0.08402381 2.240719
#>  [58,] 0.07026170 1.884538
#>  [59,] 0.08045964 2.175535
#>  [60,] 0.09580167 2.511072
#>  [61,] 0.09470489 2.533849
#>  [62,] 0.10681429 2.832647
#>  [63,] 0.06495163 1.739009
#>  [64,] 0.08706395 2.368834
#>  [65,] 0.09176504 2.467235
#>  [66,] 0.08347602 2.234426
#>  [67,] 0.10286184 2.719394
#>  [68,] 0.09143797 2.526102
#>  [69,] 0.06791179 1.799408
#>  [70,] 0.08130941 2.218032
#>  [71,] 0.06825023 1.844812
#>  [72,] 0.11440296 3.109920
#>  [73,] 0.06073428 1.648749
#>  [74,] 0.07796900 2.156834
#>  [75,] 0.07508047 2.022909
#>  [76,] 0.10677764 2.834709
#>  [77,] 0.09604369 2.571415
#>  [78,] 0.08721350 2.371596
#>  [79,] 0.06570885 1.776283
#>  [80,] 0.05785708 1.557014
#>  [81,] 0.06421274 1.714279
#>  [82,] 0.09237578 2.415500
#>  [83,] 0.07133711 1.894212
#>  [84,] 0.09829563 2.626123
#>  [85,] 0.09614009 2.600597
#>  [86,] 0.06440193 1.752572
#>  [87,] 0.09357089 2.587562
#>  [88,] 0.07474818 2.002868
#>  [89,] 0.06268655 1.699433
#>  [90,] 0.10975734 2.941041
#>  [91,] 0.10045591 2.647836
#>  [92,] 0.07292105 2.020010
#>  [93,] 0.09533099 2.538228
#>  [94,] 0.05871973 1.587018
#>  [95,] 0.06068649 1.649058
#>  [96,] 0.07445737 2.017298
#>  [97,] 0.07664515 2.096006
#>  [98,] 0.09689460 2.557056
#>  [99,] 0.09746000 2.565203
#> attr(,"bootstrap_draws")attr(,"bootstrap_ci")
#> [1] "studentized"
head(predict(fit, newdata = seq(0, 80, length.out = 6), interval = "prediction"))
#>    x      fit     se.fit working_variance           lwr          upr
#> 1  0 0.000000 0.00000000     1.203550e-11 -6.799549e-06 6.799549e-06
#> 2 16 1.113395 0.01709891     4.814198e-03  9.733357e-01 1.253455e+00
#> 3 32 1.898201 0.02294258     6.808304e-03  1.730345e+00 2.066058e+00
#> 4 48 2.481175 0.02455237     8.338436e-03  2.295844e+00 2.666505e+00
#> 5 64 2.931304 0.02479206     9.628397e-03  2.732941e+00 3.129668e+00
#> 6 80 3.289353 0.02489137     1.076487e-02  3.080229e+00 3.498477e+00

For a one-step console summary plus a 95% confidence band, use:

report_mm(fit, interval_type = "confidence")
set.seed(1)
report_mm(fit, method = "bootstrap", B = 99,
          interval_type = "confidence")

Screening working variance models

screen <- screen_mm(
  x = one_curve$s_uM,
  y = one_curve$v_uM_per_min,
  power_values = c(0.4, 0.6),
  include_auto = TRUE,
  quiet = TRUE
)

screen$table[, c("model", "selected_model", "quasi_aic", "quasi_bic", "rmse")]
#>             model selected_model quasi_aic quasi_bic      rmse
#> 1 auto(quasi_aic)    power(0.85) -46.89897 -42.90235 0.1212890
#> 2      power(0.6)     power(0.6) -45.92753 -41.93092 0.1199642
#> 3            sqrt           sqrt -44.87655 -40.87994 0.1197806
#> 4      power(0.4)     power(0.4) -43.41216 -39.41554 0.1196889
#> 5        cuberoot       cuberoot -42.20181 -38.20520 0.1196584
#> 6             log            log -40.20657 -36.20995 0.1196680
#> 7        constant       constant -33.44685 -29.45023 0.1196306

The screening table is ordered by quasi-AIC, with quasi-BIC and RMSE reported as secondary summaries.

Grouped analyses

For real enzyme panels we often fit many curves at once and then compare the best model within each group.

grouped <- group_mm(
  data = sdl_demo,
  s = "s_uM",
  v = "v_uM_per_min",
  groups = "enzyme",
  variance_models = c("constant", "log", "sqrt", "cuberoot"),
  power_values = c(0.4, 0.6),
  include_auto = TRUE,
  quiet = TRUE
)

head(grouped$comparison$best_by_group[, c("group_label", "model", "selected_model", "quasi_aic", "quasi_bic", "rmse")])
#>   group_label           model selected_model  quasi_aic  quasi_bic       rmse
#> 1 enzyme=1111 auto(quasi_aic)    power(0.85)  -46.89897  -42.90235 0.12128898
#> 2 enzyme=2222 auto(quasi_aic)       power(1) -106.50725 -102.51063 0.06072664
#> 3 enzyme=3333 auto(quasi_aic)       power(1)  -96.64475  -92.64814 0.08058515
#> 4 enzyme=4444 auto(quasi_aic)       power(1) -100.28551  -96.28890 0.06026349
#> 5 enzyme=5555 auto(quasi_aic)       power(1)  -73.10378  -69.10717 0.09232693
#> 6 enzyme=6666 auto(quasi_aic)       power(1)  -41.49399  -37.49737 0.15304469

To visualize the current best model for each group, use the grouped plotting method.

plot(grouped, interval_type = "confidence")

Clustered repeated measurements

For repeated or clustered assays, cluster_mm() uses a Ma-Genton-style working covariance with a random effect on Vmax and a low-dimensional residual variance function. The bundled alves_demo data illustrate the interface on a small soil-enzyme panel.

cluster_fit <- cluster_mm(
  data = subset(alves_demo, enzyme == "BG"),
  s = "substrate_conc",
  v = "activity",
  cluster = "core",
  variance = "sqrt"
)

summary(cluster_fit)
#> Call:
#> cluster_mm(data = subset(alves_demo, enzyme == "BG"), s = "substrate_conc", 
#>     v = "activity", cluster = "core", variance = "sqrt")
#> 
#> Method: clustered-profile-score 
#> Clusters: 3 
#> Initialization: pooled variance-aware fit 
#> Converged: yes   Iterations: 4 
#> Residual working variance: Var = gamma S^(1/2) 
#> 
#> CI method: disabled by default for sparse clustered fits
#> 
#>        Estimate
#> Vmax 1435.43759
#> Km     34.29102
#> 
#> Interval inference is disabled by default for sparse clustered fits (3 clusters, 7 distinct substrate concentrations). Default interval reporting requires at least 4 clusters and at least 6 distinct concentrations.
#> 
#> Tau^2:242368  Gamma:1303.433  RMSE:362.7108
#> Quasi-AIC:1136.631  Quasi-BIC:1146.109
confint(cluster_fit)
#> Warning: Small-cluster warning: this fit uses 3 clusters across 7 distinct
#> substrate concentrations. Default interval reporting is disabled below 4
#> clusters or below 6 distinct concentrations because coverage may be unstable.
#> Use these intervals only as a sensitivity analysis.
#>          lower     upper
#> Vmax 874.21368 1996.6615
#> Km    29.34983   39.2322
head(predict(cluster_fit, newdata = seq(0, 700, length.out = 6),
             interval = "confidence"))
#>     x      fit   se.fit working_variance marginal_variance      lwr      upr
#> 1   0    0.000   0.0000     1.303433e-05      1.303433e-05   0.0000    0.000
#> 2 140 1153.021 228.8684     1.542243e+04      1.718025e+05 704.4475 1601.595
#> 3 280 1278.823 254.2365     2.181061e+04      2.141763e+05 780.5285 1777.117
#> 4 420 1327.087 264.0631     2.671244e+04      2.338723e+05 809.5330 1844.641
#> 5 560 1352.612 269.2834     3.084487e+04      2.460502e+05 824.8260 1880.398
#> 6 700 1368.403 272.5217     3.448561e+04      2.547453e+05 834.2708 1902.536

Because the bundled example contains only three cores, the printed summary reports point estimates by default and treats interval output as an explicit sensitivity analysis.

Simulation

set.seed(100)
sim_dat <- simulate_mm_data(variance_shape = "hill", error = "skewed")
head(sim_dat)
#>   x         y true_mean true_variance      error
#> 1 1  4.065991  4.761905      1.022444 -0.6882333
#> 2 2  7.609953  9.090909      1.089109 -1.4190785
#> 3 3 13.785376 13.043478      1.198044  0.6778096
#> 4 4 16.505541 16.666667      1.346154 -0.1388731
#> 5 5 20.075119 20.000000      1.529412  0.0607418
#> 6 6 22.084203 23.076923      1.743119 -0.7519053

The package is intentionally lightweight: it relies only on base and recommended R packages, while keeping the workflow close to the heteroscedastic and clustered Michaelis-Menten analyses used in the companion methodology paper.