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.08812  6.25917  6.60459
#> Km   76.42907    2.37059 71.78279 81.07534
#> 
#> 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.10464  6.27879  6.68897
#> Km   77.97383    2.36898 73.33072 82.61694
#> 
#> 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.259171  6.60459
#> Km   71.782792 81.07534
set.seed(1)
confint(fit, method = "bootstrap", B = 99)
#>          lower     upper
#> Vmax  6.224315  6.638768
#> Km   71.080391 80.980198
#> attr(,"bootstrap_draws")
#>            Vmax       Km
#>   [1,] 6.476013 79.83446
#>   [2,] 6.378488 74.52053
#>   [3,] 6.315635 73.40742
#>   [4,] 6.606759 81.49450
#>   [5,] 6.355271 72.46983
#>   [6,] 6.490772 78.61708
#>   [7,] 6.375798 76.89579
#>   [8,] 6.415435 77.91822
#>   [9,] 6.395333 75.79797
#>  [10,] 6.452104 75.77421
#>  [11,] 6.354902 73.64767
#>  [12,] 6.253890 72.46448
#>  [13,] 6.432755 76.48557
#>  [14,] 6.239816 72.73522
#>  [15,] 6.392421 72.35687
#>  [16,] 6.438676 75.82950
#>  [17,] 6.407854 76.57698
#>  [18,] 6.418602 74.41746
#>  [19,] 6.416150 78.11961
#>  [20,] 6.398709 76.94679
#>  [21,] 6.409231 78.06684
#>  [22,] 6.414848 73.83373
#>  [23,] 6.596570 80.87162
#>  [24,] 6.445099 78.07892
#>  [25,] 6.566444 78.84443
#>  [26,] 6.316286 74.66421
#>  [27,] 6.393557 76.88738
#>  [28,] 6.458692 75.04789
#>  [29,] 6.379901 76.13781
#>  [30,] 6.257576 71.96651
#>  [31,] 6.422571 75.72082
#>  [32,] 6.464795 76.89847
#>  [33,] 6.286778 73.34995
#>  [34,] 6.471574 77.19527
#>  [35,] 6.449847 76.42107
#>  [36,] 6.371706 74.37568
#>  [37,] 6.357428 74.21839
#>  [38,] 6.357334 74.64550
#>  [39,] 6.531733 77.41515
#>  [40,] 6.372838 75.80516
#>  [41,] 6.640540 81.91559
#>  [42,] 6.427200 77.41147
#>  [43,] 6.243263 73.29176
#>  [44,] 6.580359 78.40105
#>  [45,] 6.359399 73.67903
#>  [46,] 6.546413 78.07211
#>  [47,] 6.494002 77.34733
#>  [48,] 6.368514 74.86932
#>  [49,] 6.431841 75.52364
#>  [50,] 6.451744 76.90597
#>  [51,] 6.453025 77.24754
#>  [52,] 6.484998 78.84884
#>  [53,] 6.444247 76.11353
#>  [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.76627
#>  [61,] 6.449384 76.19105
#>  [62,] 6.350371 74.00295
#>  [63,] 6.389128 75.40213
#>  [64,] 6.509620 78.61135
#>  [65,] 6.498540 77.34006
#>  [66,] 6.297256 74.08349
#>  [67,] 6.304906 73.09760
#>  [68,] 6.650319 82.14755
#>  [69,] 6.259182 72.65698
#>  [70,] 6.620325 80.46960
#>  [71,] 6.444258 77.05716
#>  [72,] 6.618748 80.11027
#>  [73,] 6.430931 77.26469
#>  [74,] 6.651050 82.28823
#>  [75,] 6.398579 76.10855
#>  [76,] 6.531708 76.63966
#>  [77,] 6.470498 76.55241
#>  [78,] 6.501118 78.43638
#>  [79,] 6.485009 77.65041
#>  [80,] 6.485932 77.24606
#>  [81,] 6.415493 75.51380
#>  [82,] 6.443931 74.05106
#>  [83,] 6.325003 73.76093
#>  [84,] 6.323351 74.28235
#>  [85,] 6.537226 78.46156
#>  [86,] 6.425374 77.41009
#>  [87,] 6.631429 81.96381
#>  [88,] 6.380613 75.35240
#>  [89,] 6.483115 77.88983
#>  [90,] 6.338269 74.75711
#>  [91,] 6.560837 76.39026
#>  [92,] 6.553474 80.97582
#>  [93,] 6.393650 74.96486
#>  [94,] 6.500069 77.84720
#>  [95,] 6.469865 77.91694
#>  [96,] 6.341839 75.80376
#>  [97,] 6.459499 78.36328
#>  [98,] 6.358715 73.68683
#>  [99,] 6.322228 72.95036
#> 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.06989077 1.936539
#>   [2,] 0.10251590 2.722894
#>   [3,] 0.06738529 1.785029
#>   [4,] 0.08951382 2.473210
#>   [5,] 0.08394792 2.186288
#>   [6,] 0.08792374 2.399725
#>   [7,] 0.08220107 2.242207
#>   [8,] 0.08117004 2.224783
#>   [9,] 0.06968204 1.872313
#>  [10,] 0.06273594 1.670407
#>  [11,] 0.09727914 2.567992
#>  [12,] 0.07318858 1.936858
#>  [13,] 0.06390245 1.719946
#>  [14,] 0.06910804 1.838718
#>  [15,] 0.08946443 2.313397
#>  [16,] 0.06838104 1.825621
#>  [17,] 0.06265429 1.694593
#>  [18,] 0.09248357 2.438256
#>  [19,] 0.09934208 2.728423
#>  [20,] 0.07680999 2.088807
#>  [21,] 0.07702599 2.116604
#>  [22,] 0.09012642 2.361914
#>  [23,] 0.08551788 2.351381
#>  [24,] 0.07813264 2.135340
#>  [25,] 0.07840244 2.120297
#>  [26,] 0.07508338 2.017148
#>  [27,] 0.06360233 1.729910
#>  [28,] 0.11555513 3.049014
#>  [29,] 0.07792122 2.106610
#>  [30,] 0.09942731 2.614573
#>  [31,] 0.05740997 1.534723
#>  [32,] 0.08073298 2.171909
#>  [33,] 0.06927707 1.842359
#>  [34,] 0.12746778 3.436621
#>  [35,] 0.07652956 2.052902
#>  [36,] 0.06969284 1.850052
#>  [37,] 0.07618688 2.023402
#>  [38,] 0.06886232 1.837687
#>  [39,] 0.11088210 2.968962
#>  [40,] 0.06947000 1.873353
#>  [41,] 0.08290897 2.288869
#>  [42,] 0.06899230 1.877298
#>  [43,] 0.08401557 2.248397
#>  [44,] 0.11013911 2.958334
#>  [45,] 0.10514683 2.774708
#>  [46,] 0.11423646 3.073505
#>  [47,] 0.06591382 1.773855
#>  [48,] 0.07499842 2.002926
#>  [49,] 0.05972990 1.590973
#>  [50,] 0.10325916 2.783762
#>  [51,] 0.06371978 1.723841
#>  [52,] 0.06437887 1.762995
#>  [53,] 0.06790257 1.816941
#>  [54,] 0.10346651 2.754202
#>  [55,] 0.07594883 2.037174
#>  [56,] 0.10585135 2.775861
#>  [57,] 0.08406349 2.242129
#>  [58,] 0.07029539 1.885742
#>  [59,] 0.08049980 2.176977
#>  [60,] 0.09584361 2.512544
#>  [61,] 0.09475082 2.535482
#>  [62,] 0.10686321 2.834380
#>  [63,] 0.06498273 1.740118
#>  [64,] 0.08710949 2.370474
#>  [65,] 0.09181090 2.468871
#>  [66,] 0.08351503 2.235820
#>  [67,] 0.10290794 2.721024
#>  [68,] 0.09149037 2.528006
#>  [69,] 0.06794211 1.800483
#>  [70,] 0.08135356 2.219621
#>  [71,] 0.06828453 1.846042
#>  [72,] 0.11446427 3.112120
#>  [73,] 0.06076518 1.649863
#>  [74,] 0.07801390 2.158467
#>  [75,] 0.07511733 2.024229
#>  [76,] 0.10682903 2.836522
#>  [77,] 0.09609063 2.573085
#>  [78,] 0.08725895 2.373233
#>  [79,] 0.06574222 1.777480
#>  [80,] 0.05788599 1.558047
#>  [81,] 0.06424337 1.715368
#>  [82,] 0.09241696 2.416938
#>  [83,] 0.07136973 1.895368
#>  [84,] 0.09834156 2.627760
#>  [85,] 0.09618969 2.602374
#>  [86,] 0.06443494 1.753763
#>  [87,] 0.09362445 2.589511
#>  [88,] 0.07478399 2.004146
#>  [89,] 0.06271870 1.700589
#>  [90,] 0.10980937 2.942900
#>  [91,] 0.10050336 2.649498
#>  [92,] 0.07296224 2.021514
#>  [93,] 0.09537578 2.539818
#>  [94,] 0.05874964 1.588089
#>  [95,] 0.06071779 1.650185
#>  [96,] 0.07449414 2.018624
#>  [97,] 0.07668550 2.097467
#>  [98,] 0.09693830 2.558596
#>  [99,] 0.09750317 2.566723
#> 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.01710788     4.814198e-03  9.733315e-01 1.253459e+00
#> 3 32 1.898201 0.02295280     6.808304e-03  1.730339e+00 2.066063e+00
#> 4 48 2.481175 0.02456079     8.338436e-03  2.295840e+00 2.666510e+00
#> 5 64 2.931304 0.02479772     9.628397e-03  2.732938e+00 3.129671e+00
#> 6 80 3.289353 0.02489441     1.076487e-02  3.080227e+00 3.498479e+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.

User-defined fixed variance functions

For a single curve, a function can replace a built-in variance name. It defines the shape h(S) in Var(Y | S) = gamma h(S), not the standard deviation.

h <- variance_function(function(S) 1 + S / 20,
                       label = "Var = gamma (1 + S/20)")
custom_fit <- fit_mm(one_curve$s_uM, one_curve$v_uM_per_min, variance = h)
summary(custom_fit)
#> Call:
#> fit_mm(x = one_curve$s_uM, y = one_curve$v_uM_per_min, variance = h)
#> 
#> Method: profile-score 
#> Working variance: Var = gamma (1 + S/20) 
#> 
#> CI method: wald
#> 
#>      Estimate Std. Error Lower CI Upper CI
#> Vmax  6.44671    0.10224  6.24631  6.64710
#> Km   76.67404    2.41534 71.94007 81.40801
#> 
#> Scale:0.0019  RMSE:0.11999  Weighted RSS:28
#> Quasi-AIC:-45.57156  Quasi-BIC:-41.57495
confint(custom_fit)
#>          lower     upper
#> Vmax  6.246315  6.647104
#> Km   71.940071 81.408014
predict(custom_fit, newdata = c(0, 10, 40), interval = "prediction")
#>    x       fit     se.fit working_variance         lwr        upr
#> 1  0 0.0000000 0.00000000      0.001899016 -0.08541073 0.08541073
#> 2 10 0.7437878 0.01144077      0.002848524  0.63680508 0.85077051
#> 3 40 2.2101607 0.02121135      0.005697049  2.05649436 2.36382696

custom_shapes <- list(
  linear_floor = function(S) 1 + S / 20,
  saturation = function(S) 1 + 9 * S / (20 + S)
)
custom_screen <- screen_mm(
  one_curve$s_uM, one_curve$v_uM_per_min,
  variance_models = c("constant", "sqrt"),
  variance_functions = custom_shapes, quiet = TRUE
)
custom_screen$table[, c("model", "quasi_aic", "quasi_bic", "success")]
#>          model quasi_aic quasi_bic success
#> 1 linear_floor -45.57156 -41.57495       1
#> 2         sqrt -44.87655 -40.87994       1
#> 3   saturation -39.55276 -35.55615       1
#> 4     constant -33.44685 -29.45023       1

The function must be deterministic, vectorized, and pointwise, returning one finite strictly positive numeric value for each input concentration, including zero if used. Values must not depend on the other concentrations in the vector. For example, a constant shape should return rep(1, length(S)), not a scalar. Invalid values or lengths cause errors; screening records candidate-specific errors and continues with the remaining candidates. Names must be unique and must not duplicate built-in candidates. For custom-only screening, specify variance_models = character(0).

Only Km requires numerical search for a fixed shape. Vmax and gamma retain their closed-form updates, and the function is reused for prediction and bootstrap refits. Internal scaling preserves the supplied variance model. The reported gamma is relative to the original shape, so its value changes inversely when the shape is multiplied by a positive constant.

The constants 20 and 9 above are fixed, not estimated. Extra shape parameters, functions of unknown mean parameters, and data-adaptive nuisance estimation are not fitted by this interface. Intervals and information criteria condition on the supplied shape and omit external shape-estimation or selection uncertainty; prediction bands still use the working-normal approximation. User-defined variance shapes are not higher-moment efficient scores and do not imply smaller variance than identically weighted WNLS. They are currently supported for independent curves, not the clustered routine.

Optional higher-moment estimation for independent curves

The user-defined variance shape above does not by itself incorporate asymmetry into the estimating equations. Version 0.1.1 adds the separate fit_mm_moments() estimator, based on the conditional-moment construction of Kim and Ma (2012), Section 3 (doi:10.1007/s10463-011-0332-y). It does not require a fully specified error density.

set.seed(19)
skew_dat <- simulate_mm_data(n = 200, error = "skewed", skew_shape = 4)
assay_shape <- function(S) 1 + 9 * S / (20 + S)
fit_hm <- fit_mm_moments(skew_dat$x, skew_dat$y, variance = assay_shape)
summary(fit_hm)
#> Call:
#> fit_mm_moments(x = skew_dat$x, y = skew_dat$y, variance = assay_shape)
#> 
#> Method: higher-moment estimating equations 
#> Working variance: Var = gamma h(S) (user-defined) 
#> Working moments: estimate 
#> Covariance: model 
#> CI method: wald 
#> 
#>      Estimate Std. Error Lower CI  Upper CI
#> Vmax 99.32847    0.58349 98.18485 100.47210
#> Km   19.60062    0.36175 18.89160  20.30965
#> 
#> Variance scale: 0.9439  RMSE: 2.476 
#> Scaled score norm: 8.5468e-10
confint(fit_hm)
#>         lower     upper
#> Vmax 98.18485 100.47210
#> Km   18.89160  20.30965
fit_hm$vcov_joint
#>             Vmax           Km        gamma
#> Vmax  0.34046250  0.186485439  0.012018354
#> Km    0.18648544  0.130864749 -0.005169235
#> gamma 0.01201835 -0.005169235  0.020383201

fit_gamma <- fit_mm_moments(
  skew_dat$x, skew_dat$y, variance = assay_shape,
  moments = "gamma", shape = 4
)
coef(fit_gamma)
#>     Vmax       Km 
#> 99.40513 19.64990
predict(fit_gamma, c(10, 20, 80), interval = "confidence")
#>    x      fit    se.fit working_variance      lwr      upr
#> 1 10 33.52629 0.2499221         3.771718 33.03645 34.01613
#> 2 20 50.14142 0.2391705         5.186112 49.67266 50.61019
#> 3 80 79.80349 0.2502522         7.732022 79.31301 80.29398

The default estimates constant standardized skewness and kurtosis from the profile pilot’s residuals. The gamma option supplies fixed working moments of centered gamma errors; it does not fit a gamma density. Student-t working moments use moments = "student_t", df = 10. Custom standardized moment functions can be passed with moments = "custom", skewness, and kurtosis; the latter is raw kurtosis, not excess kurtosis.

For residual $r_i$ and variance $v_i=\gamma h(S_i)$, let $u_i=(r_i,r_i^2-v_i)^\top$. The conditional moment covariance is

\[ \Omega_i=\begin{pmatrix} v_i & s_i v_i^{3/2}\\ s_i v_i^{3/2} & (k_i-1)v_i^2 \end{pmatrix}, \]

where $s_i$ and $k_i$ are working skewness and kurtosis. With $g_i$ the MM mean gradient and $D_i=\begin{pmatrix}g_{i1}&g_{i2}&0\\0&0&h(S_i)\end{pmatrix}$, the joint equation is $\sum_i D_i^\top\Omega_i^{-1}u_i=0$. Positive definiteness requires $k_i-1-s_i^2>0$.

Unlike the default profile route, the higher-moment fit solves jointly for $V_{\max},K_m,\gamma$ numerically. Its default covariance = "model" uses joint inverse working information, with the mean covariance extracted after allowing for the estimated scale. This interpretation requires appropriate conditional third and fourth moments. Optional covariance = "sandwich" uses the complete three-parameter equations. Wald intervals use a normal cutoff for the estimator, not an assumption that the individual errors are normal.

Normal moments recover the original mean/variance equations. Constant symmetric Student-t moments recover the same point estimates; heavy tails alone do not give an efficiency improvement in this conditional-moment model. Under asymmetry, appropriate third/fourth moments can improve asymptotic efficiency, but pilot moment estimation can inflate finite-sample variance. Incorrect variance shapes can bias the new estimator, and a sandwich does not remove that bias. Neither exact 95% coverage nor uniform improvement is guaranteed.

set.seed(20)
confint(fit_hm, method = "bootstrap", B = 399)  # pairs only
report_mm(fit_hm, interval_type = "confidence")

Pairs bootstrap repeats the pilot and higher-moment estimation. Prediction bands are not supplied because a few moments do not identify predictive quantiles; mean-confidence bands are pointwise. Quasi-AIC/BIC are omitted for these non-likelihood roots. This opt-in estimator is not yet connected to grouped selection or clustered fitting; those workflows below retain their original estimators.

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,
  variance_functions = custom_shapes,
  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.12128901
#> 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.06026352
#> 5 enzyme=5555 auto(quasi_aic)       power(1)  -73.10378  -69.10717 0.09232694
#> 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

By default, cluster_mm() uses a rank-one 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 single-curve profile-score 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.

An optional method = "projected" uses an MM-specific adaptation of Ma and Genton (2010, Section 3.2), https://doi.org/10.1111/j.1467-9868.2010.00741.x. It leaves the random-effect density unspecified while retaining conditionally independent Gaussian residuals with a correct fixed variance shape. The identity-instrument projected score does not generally attain the semiparametric efficiency bound. Standard errors use a stacked independent-cluster sandwich, and only mean-confidence bands are returned. At least five clusters and two readings per cluster are required; this path is therefore not applied to the three-core example above. The package README provides a reproducible skewed-random-effect example. This is distinct from non-Gaussian independent residual estimation.

Projected failures report a numerical stage, reason and available diagnostics in inferMM_cluster_projected_error; no estimate or interval is returned and no automatic method switch occurs. An incomplete numerical search can retain a verified interior root but warns explicitly. Inspect fit$convergence, especially the scaled joint score norm and search diagnostics. A non-positive variance-component moment has a separate reporting warning. For the default profile method, an iteration-limit warning marks converged = FALSE; such results should not be interpreted as converged fits.