Cardinal parameter models

library(predmicror)

Cardinal parameter models describe the effect of environmental factors on a microbial growth response.

In predmicror, the cardinal fitting wrapper expects the response to be the square root of the growth rate, usually in a column named sqrtGR. This vignette shows how to fit cardinal models for temperature, pH, water activity, and inhibitory compounds.

Available cardinal models

predmicror_models("cardinal")
#>       type model response_scale                 parameters
#> 1 cardinal  CMTI         sqrtGR    Tmax, Tmin, MUopt, Topt
#> 2 cardinal  CMAW         sqrtGR        AWmin, MUopt, AWopt
#> 3 cardinal  CMPH         sqrtGR pHmax, pHmin, MUopt, pHopt
#> 4 cardinal CMInh         sqrtGR          MIC, MUopt, alpha

Robust fitting helper

Nonlinear fitting can depend on starting values. The helper below keeps the vignette robust by reporting fitting problems without stopping the document build.

safe_fit <- function(expr) {
  withCallingHandlers(
    tryCatch(
      expr,
      error = function(e) {
        message("Model fit failed: ", conditionMessage(e))
        NULL
      }
    ),
    warning = function(w) {
      message("Model fit warning: ", conditionMessage(w))
      invokeRestart("muffleWarning")
    }
  )
}

Temperature model

The salmonella dataset contains temperature and growth-rate data for Salmonella typhimurium on cooked chicken.

data(salmonella)
head(salmonella)
#>   Temp    GR sqrtGR
#> 1    8 0.012  0.110
#> 2   10 0.025  0.158
#> 3   12 0.048  0.219
#> 4   14 0.078  0.279
#> 5   16 0.107  0.327
#> 6   18 0.155  0.394

temp_fit <- safe_fit(
  fit_cardinal(
    salmonella,
    model = "CMTI",
    x = "Temp",
    response = "sqrtGR",
    start = list(Tmax = 42, Tmin = 1, MUopt = 1, Topt = 37)
  )
)

temp_fit
#> predmicror fit
#>   type: cardinal
#>   model: CMTI
#>   response: sqrtGR (sqrtGR)
#>   formula: sqrtGR ~ CMTI(Temp, Tmax, Tmin, MUopt, Topt)
#> 
#> Nonlinear regression model
#>   model: sqrtGR ~ CMTI(Temp, Tmax, Tmin, MUopt, Topt)
#>    data: data
#>    Tmax    Tmin   MUopt    Topt 
#> 49.0326  5.0268  0.7521 40.3671 
#>  residual sum-of-squares: 0.005018
#> 
#> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr)
#> 
#> Number of iterations to convergence: 44 
#> Achieved convergence tolerance: 7.98e-17

if (!is.null(temp_fit)) {
  fit_metrics(temp_fit)
}
#>   model     type response response_scale  n p         SSE       RMSE        MAE
#> 1  CMTI cardinal   sqrtGR         sqrtGR 21 4 0.005018168 0.01545834 0.01253912
#>            bias        RSE        R2    adj_R2   logLik       AIC       BIC
#> 1 -4.168898e-05 0.01718099 0.9960034 0.9950043 57.76403 -105.5281 -100.3054
#>   converged
#> 1      TRUE

if (!is.null(temp_fit)) {
  grid <- data.frame(
    Temp = seq(min(salmonella$Temp), max(salmonella$Temp), length.out = 100)
  )
  pred <- predmicror_augment(temp_fit, newdata = grid)

  plot(
    sqrtGR ~ Temp,
    data = salmonella,
    xlab = "Temperature",
    ylab = "Square root growth rate",
    pch = 19
  )
  if (".fitted" %in% names(pred)) {
    lines(pred$Temp, pred[[".fitted"]], lwd = 2)
  }
}

pH model

data(ph)
head(ph)
#>    pH sqrtGR     GR
#> 1 4.4   0.14 0.0196
#> 2 4.5   0.20 0.0400
#> 3 4.6   0.39 0.1521
#> 4 4.8   0.46 0.2116
#> 5 5.0   0.58 0.3364
#> 6 5.1   0.59 0.3481

ph_fit <- safe_fit(
  fit_cardinal(
    ph,
    model = "CMPH",
    x = "pH",
    response = "sqrtGR",
    start = list(pHmax = 9.5, pHmin = 3.5, MUopt = 1, pHopt = 6.8)
  )
)

ph_fit
#> predmicror fit
#>   type: cardinal
#>   model: CMPH
#>   response: sqrtGR (sqrtGR)
#>   formula: sqrtGR ~ CMPH(pH, pHmax, pHmin, MUopt, pHopt)
#> 
#> Nonlinear regression model
#>   model: sqrtGR ~ CMPH(pH, pHmax, pHmin, MUopt, pHopt)
#>    data: data
#> pHmax pHmin MUopt pHopt 
#> 9.355 4.425 1.066 7.257 
#>  residual sum-of-squares: 0.03246
#> 
#> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr)
#> 
#> Number of iterations to convergence: 9 
#> Achieved convergence tolerance: 6.245e-17

if (!is.null(ph_fit)) {
  fit_metrics(ph_fit)
}
#>   model     type response response_scale  n p        SSE       RMSE        MAE
#> 1  CMPH cardinal   sqrtGR         sqrtGR 14 4 0.03245538 0.04814812 0.03263188
#>         bias        RSE        R2    adj_R2   logLik       AIC       BIC
#> 1 0.01020897 0.05696963 0.9664767 0.9515774 22.60348 -35.20697 -32.01168
#>   converged
#> 1      TRUE

if (!is.null(ph_fit)) {
  grid <- data.frame(pH = seq(min(ph$pH), max(ph$pH), length.out = 100))
  pred <- predmicror_augment(ph_fit, newdata = grid)

  plot(
    sqrtGR ~ pH,
    data = ph,
    xlab = "pH",
    ylab = "Square root growth rate",
    pch = 19
  )
  if (".fitted" %in% names(pred)) {
    lines(pred$pH, pred[[".fitted"]], lwd = 2)
  }
}

Water activity model

data(aw)
head(aw)
#>      aw sqrtGR       GR
#> 1 0.765  0.036 0.001296
#> 2 0.790  0.240 0.057600
#> 3 0.820  0.384 0.147456
#> 4 0.862  0.792 0.627264
#> 5 0.895  1.082 1.170724
#> 6 0.920  1.204 1.449616

aw_fit <- safe_fit(
  fit_cardinal(
    aw,
    model = "CMAW",
    x = "aw",
    response = "sqrtGR",
    start = list(AWmin = 0.90, MUopt = 1, AWopt = 0.995)
  )
)

aw_fit
#> predmicror fit
#>   type: cardinal
#>   model: CMAW
#>   response: sqrtGR (sqrtGR)
#>   formula: sqrtGR ~ CMAW(aw, AWmin, MUopt, AWopt)
#> 
#> Nonlinear regression model
#>   model: sqrtGR ~ CMAW(aw, AWmin, MUopt, AWopt)
#>    data: data
#>  AWmin  MUopt  AWopt 
#> 0.7717 1.5586 0.9460 
#>  residual sum-of-squares: 0.01396
#> 
#> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr)
#> 
#> Number of iterations to convergence: 12 
#> Achieved convergence tolerance: 6.072e-17

if (!is.null(aw_fit)) {
  fit_metrics(aw_fit)
}
#>   model     type response response_scale n p        SSE       RMSE        MAE
#> 1  CMAW cardinal   sqrtGR         sqrtGR 9 3 0.01396232 0.03938741 0.03342219
#>          bias        RSE       R2   adj_R2   logLik       AIC       BIC
#> 1 0.003427043 0.04823953 0.991345 0.986152 16.33833 -24.67667 -23.88777
#>   converged
#> 1      TRUE

if (!is.null(aw_fit)) {
  grid <- data.frame(aw = seq(min(aw$aw), max(aw$aw), length.out = 100))
  pred <- predmicror_augment(aw_fit, newdata = grid)

  plot(
    sqrtGR ~ aw,
    data = aw,
    xlab = "Water activity",
    ylab = "Square root growth rate",
    pch = 19
  )
  if (".fitted" %in% names(pred)) {
    lines(pred$aw, pred[[".fitted"]], lwd = 2)
  }
}

Inhibitor model

data(inh)
head(inh)
#>   Conce sqrtGR       GR
#> 1   0.0  0.864 0.746496
#> 2   0.6  0.704 0.495616
#> 3   1.0  0.683 0.466489
#> 4   1.9  0.523 0.273529
#> 5   2.7  0.384 0.147456
#> 6   3.8  0.338 0.114244

inh_fit <- safe_fit(
  fit_cardinal(
    inh,
    model = "CMInh",
    x = "Conce",
    response = "sqrtGR",
    start = list(MIC = max(inh$Conce) * 1.2, MUopt = 1, alpha = 1)
  )
)

inh_fit
#> predmicror fit
#>   type: cardinal
#>   model: CMInh
#>   response: sqrtGR (sqrtGR)
#>   formula: sqrtGR ~ CMInh(Conce, MIC, MUopt, alpha)
#> 
#> Nonlinear regression model
#>   model: sqrtGR ~ CMInh(Conce, MIC, MUopt, alpha)
#>    data: data
#>    MIC  MUopt  alpha 
#> 6.3231 0.7780 0.3569 
#>  residual sum-of-squares: 0.0166
#> 
#> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr)
#> 
#> Number of iterations to convergence: 14 
#> Achieved convergence tolerance: 2.082e-16

if (!is.null(inh_fit)) {
  fit_metrics(inh_fit)
}
#>   model     type response response_scale n p        SSE       RMSE        MAE
#> 1 CMInh cardinal   sqrtGR         sqrtGR 8 3 0.01659573 0.04554632 0.03882655
#>           bias        RSE       R2    adj_R2  logLik       AIC       BIC
#> 1 -0.003475956 0.05761204 0.961369 0.9323957 13.3607 -18.72139 -18.40363
#>   converged
#> 1      TRUE

if (!is.null(inh_fit)) {
  grid <- data.frame(Conce = seq(min(inh$Conce), max(inh$Conce), length.out = 100))
  pred <- predmicror_augment(inh_fit, newdata = grid)

  plot(
    sqrtGR ~ Conce,
    data = inh,
    xlab = "Inhibitor concentration",
    ylab = "Square root growth rate",
    pch = 19
  )
  if (".fitted" %in% names(pred)) {
    lines(pred$Conce, pred[[".fitted"]], lwd = 2)
  }
}

Collecting diagnostics

The fitted models above describe different environmental factors and should not be ranked against each other using information criteria. Still, collecting their metrics in one table is useful for reporting and quality control.

cardinal_fits <- Filter(Negate(is.null), list(
  temperature = temp_fit,
  pH = ph_fit,
  water_activity = aw_fit,
  inhibitor = inh_fit
))

if (length(cardinal_fits) > 0L) {
  diagnostics <- do.call(
    rbind,
    Map(function(name, fit) {
      out <- fit_metrics(fit)
      out$workflow <- name
      out
    }, names(cardinal_fits), cardinal_fits)
  )
  rownames(diagnostics) <- NULL
  diagnostics
}
#>   model     type response response_scale  n p         SSE       RMSE        MAE
#> 1  CMTI cardinal   sqrtGR         sqrtGR 21 4 0.005018168 0.01545834 0.01253912
#> 2  CMPH cardinal   sqrtGR         sqrtGR 14 4 0.032455382 0.04814812 0.03263188
#> 3  CMAW cardinal   sqrtGR         sqrtGR  9 3 0.013962315 0.03938741 0.03342219
#> 4 CMInh cardinal   sqrtGR         sqrtGR  8 3 0.016595735 0.04554632 0.03882655
#>            bias        RSE        R2    adj_R2   logLik        AIC        BIC
#> 1 -4.168898e-05 0.01718099 0.9960034 0.9950043 57.76403 -105.52805 -100.30544
#> 2  1.020897e-02 0.05696963 0.9664767 0.9515774 22.60348  -35.20697  -32.01168
#> 3  3.427043e-03 0.04823953 0.9913450 0.9861520 16.33833  -24.67667  -23.88777
#> 4 -3.475956e-03 0.05761204 0.9613690 0.9323957 13.36070  -18.72139  -18.40363
#>   converged       workflow
#> 1      TRUE    temperature
#> 2      TRUE             pH
#> 3      TRUE water_activity
#> 4      TRUE      inhibitor

Practical notes

Cardinal model parameters should be interpreted in the environmental units of the input data. For example, Tmin, Topt, and Tmax are in the same units as Temp, while AWmin and AWopt are in water-activity units.

Always check that the response scale is sqrtGR before fitting.