The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Getting Started with rumenGP

Introduction

rumenGP provides a complete workflow for analyzing in vitro rumen gas production experiments.

The package supports:

ANKOM RF datasets
Manual gas-volume datasets
Pressure-based datasets
Twelve built-in kinetic models
User-defined kinetic models
Model comparison and ranking
Treatment-level model evaluation
Diagnostic and visualization tools

This vignette demonstrates a complete workflow using the packaged ANKOM example dataset.

library(rumenGP)

Load Example Data

The package includes a small example dataset.

files <- example_data()

files
#> $ankom
#> [1] "C:/Users/Arlan/AppData/Local/Temp/RtmpQ9phut/Rinst21142159614e/rumenGP/extdata/example_ankom.xlsx"
#> 
#> $metadata
#> [1] "C:/Users/Arlan/AppData/Local/Temp/RtmpQ9phut/Rinst21142159614e/rumenGP/extdata/example_metadata.xlsx"

Import ANKOM Data

Import the ANKOM RF output file and metadata table.

raw_data <- read_ankom(
  files$ankom
)

metadata <- read_metadata(
  files$metadata
)

Validate Metadata

Before processing data, validate the metadata table.

metadata <- validate_metadata(
  metadata
)
#> Metadata validation passed.
#> Heads: 24
#> Treatment: 5

Process ANKOM Data

Convert pressure measurements into cumulative gas production.

gp <- process_ankom(
  raw_data,
  metadata,
  headspace_ml = 210,
  temperature_c = 39,
  zero_negative_pressure = TRUE
)

Validate Processed Data

The resulting dataset is a standardized rumen_gp object.

gp <- validate_ankom(
  gp
)
#> rumenGP data validation passed.
#> Observations: 1752
#> Heads: 24
#> Treatments: 5

class(gp)
#> [1] "rumen_gp"   "tbl_df"     "tbl"        "data.frame"

Inspect the data:

head(gp)
#> # A tibble: 6 × 11
#>   time_raw Time_h Head  Gas_PSI Gas_kPa Gas_moles Gas_mL Bottle   Rep Treatment
#>   <chr>     <dbl> <chr>   <dbl>   <dbl>     <dbl>  <dbl>  <dbl> <dbl> <chr>    
#> 1 10:29:34      0 1           0       0         0      0      1     1 Plant_A  
#> 2 10:29:34      0 2           0       0         0      0      2     2 Plant_A  
#> 3 10:29:34      0 3           0       0         0      0      3     3 Plant_A  
#> 4 10:29:34      0 4           0       0         0      0      4     4 Plant_A  
#> 5 10:29:34      0 5           0       0         0      0      5     5 Plant_A  
#> 6 10:29:34      0 6           0       0         0      0      6     6 Plant_A  
#> # ℹ 1 more variable: pH <dbl>

Visualize Raw Gas Production

Individual bottle profiles can be visualized.

plot_gp(
  gp,
  head = "1"
)

Fit Kinetic Models

Several built-in models are available.

groot_fit <- fit_groot(gp)
#> rumenGP data validation passed.
#> Observations: 1752
#> Heads: 24
#> Treatments: 5

gompertz_fit <- fit_gompertz(gp)
#> rumenGP data validation passed.
#> Observations: 1752
#> Heads: 24
#> Treatments: 5

brody_fit <- fit_brody(gp)
#> rumenGP data validation passed.
#> Observations: 1752
#> Heads: 24
#> Treatments: 5

Summarize Model Fits

Each model provides parameter estimates and diagnostic statistics.

summary(groot_fit)
#> 
#> Groot model summary
#> -------------------
#> Total bottles: 24
#> Successful fits: 24
#> Failed fits: 0
#> Low R-squared (< 0.90): 3

Identify Potentially Problematic Bottles

flags <- flag_model(
  groot_fit
)

head(flags)
#>   Head Bottle Rep Treatment Converged Status       RSS        R2      RMSE
#> 1    1      1   1   Plant_A      TRUE     OK  34.55456 0.9990094 0.6927658
#> 2   10     10   4   Plant_B      TRUE     OK 111.43398 0.1535773 1.2440636
#> 3   11     11   5   Plant_B      TRUE     OK 158.20900 0.9956744 1.4823452
#> 4   12     12   6   Plant_B      TRUE     OK 201.79056 0.9948662 1.6741107
#> 5   13     13   1   Plant_C      TRUE     OK 240.73021 0.9946908 1.8285172
#> 6   14     14   2   Plant_C      TRUE     OK 415.12645 0.9908573 2.4011758
#>        AIC      BIC Lambda_Boundary Flag_FitFailed Flag_LowR2
#> 1 159.4700 168.5767           FALSE          FALSE      FALSE
#> 2 243.7743 252.8810           FALSE          FALSE       TRUE
#> 3 269.0092 278.1159           FALSE          FALSE      FALSE
#> 4 286.5278 295.6344           FALSE          FALSE      FALSE
#> 5 299.2319 308.3386           FALSE          FALSE      FALSE
#> 6 338.4652 347.5718           FALSE          FALSE      FALSE
#>   Flag_LambdaBoundary Overall_Flag
#> 1               FALSE           OK
#> 2               FALSE       LOW_R2
#> 3               FALSE           OK
#> 4               FALSE           OK
#> 5               FALSE           OK
#> 6               FALSE           OK

Plot Model Fits

Observed and predicted values can be visualized.

plot_fit(
  groot_fit,
  head = "1"
)

Plot Residuals

Residual plots help identify systematic deviations.

plot_residuals(
  groot_fit,
  head = "1"
)

Compare Models

Compare model performance using multiple metrics.

comparison <- compare_models(

  Groot = groot_fit,

  Gompertz = gompertz_fit,

  Brody = brody_fit

)

comparison
#>      Model Bottles Successful_Fits Failed_Fits   Mean_R2 Mean_RMSE  Mean_RSS
#> 1    Groot      24              24           0 0.9134994  2.154301  491.4637
#> 2 Gompertz      24              23           1 0.9545965  4.216094 1837.5449
#> 3    Brody      24              23           1 0.9453304  3.604892 1076.8697
#>   Mean_AIC Mean_BIC Lambda_Boundary
#> 1 297.0687 306.1754               0
#> 2 392.2245 401.3864               8
#> 3 393.6300 402.7918               0

The comparison table includes:

Mean R-squared
Mean RMSE
Mean RSS
Mean AIC
Mean BIC
Number of successful fits

Rank Models

rank_models(
  comparison
)
#>      Model Bottles Successful_Fits Failed_Fits   Mean_R2 Mean_RMSE  Mean_RSS
#> 1    Groot      24              24           0 0.9134994  2.154301  491.4637
#> 2 Gompertz      24              23           1 0.9545965  4.216094 1837.5449
#> 3    Brody      24              23           1 0.9453304  3.604892 1076.8697
#>   Mean_AIC Mean_BIC Lambda_Boundary Rank_R2 Rank_RMSE Rank_AIC Rank_BIC
#> 1 297.0687 306.1754               0       3         1        1        1
#> 2 392.2245 401.3864               8       1         3        2        2
#> 3 393.6300 402.7918               0       2         2        3        3

Compare Models by Treatment

Treatment-level comparisons are also available.

treatment_comparison <-
  compare_models_by_treatment(

    Groot = groot_fit,

    Gompertz = gompertz_fit,

    Brody = brody_fit

  )

treatment_comparison
#> # A tibble: 15 × 6
#>    Treatment Model    Mean_R2 Mean_RMSE Mean_AIC Mean_BIC
#>    <chr>     <chr>      <dbl>     <dbl>    <dbl>    <dbl>
#>  1 BLANK     Groot      0.714      1.75     257.     266.
#>  2 Plant_A   Groot      0.968      2.33     281.     290.
#>  3 Plant_B   Groot      0.854      1.73     281.     290.
#>  4 Plant_C   Groot      0.982      2.11     309.     319.
#>  5 TMR       Groot      0.988      3.14     377.     386.
#>  6 BLANK     Gompertz   0.943      2.01     293.     302.
#>  7 Plant_A   Gompertz   0.936      5.05     417.     426.
#>  8 Plant_B   Gompertz   0.969      4.16     393.     402.
#>  9 Plant_C   Gompertz   0.967      3.73     381.     390.
#> 10 TMR       Gompertz   0.957      5.83     464.     473.
#> 11 BLANK     Brody      0.907      1.77     299.     308.
#> 12 Plant_A   Brody      0.928      4.01     405.     415.
#> 13 Plant_B   Brody      0.976      3.86     411.     421.
#> 14 Plant_C   Brody      0.941      3.46     394.     403.
#> 15 TMR       Brody      0.977      4.50     435.     444.

Rank Models by Treatment

ranked_treatments <-
  rank_models_by_treatment(
    treatment_comparison
  )

ranked_treatments
#> # A tibble: 15 × 10
#>    Treatment Model    Mean_R2 Mean_RMSE Mean_AIC Mean_BIC Rank_R2 Rank_RMSE
#>    <chr>     <chr>      <dbl>     <dbl>    <dbl>    <dbl>   <int>     <int>
#>  1 BLANK     Groot      0.714      1.75     257.     266.       3         1
#>  2 Plant_A   Groot      0.968      2.33     281.     290.       1         1
#>  3 Plant_B   Groot      0.854      1.73     281.     290.       3         1
#>  4 Plant_C   Groot      0.982      2.11     309.     319.       1         1
#>  5 TMR       Groot      0.988      3.14     377.     386.       1         1
#>  6 BLANK     Gompertz   0.943      2.01     293.     302.       1         3
#>  7 Plant_A   Gompertz   0.936      5.05     417.     426.       2         3
#>  8 Plant_B   Gompertz   0.969      4.16     393.     402.       2         3
#>  9 Plant_C   Gompertz   0.967      3.73     381.     390.       2         3
#> 10 TMR       Gompertz   0.957      5.83     464.     473.       3         3
#> 11 BLANK     Brody      0.907      1.77     299.     308.       2         2
#> 12 Plant_A   Brody      0.928      4.01     405.     415.       3         2
#> 13 Plant_B   Brody      0.976      3.86     411.     421.       1         2
#> 14 Plant_C   Brody      0.941      3.46     394.     403.       3         2
#> 15 TMR       Brody      0.977      4.50     435.     444.       2         2
#> # ℹ 2 more variables: Rank_AIC <int>, Rank_BIC <int>

Determine the Best Model per Treatment

best_models <-
  best_model_by_treatment(
    ranked_treatments
  )

best_models
#> # A tibble: 5 × 12
#>   Treatment Model Mean_R2 Mean_RMSE Mean_AIC Mean_BIC Rank_R2 Rank_RMSE Rank_AIC
#>   <chr>     <chr>   <dbl>     <dbl>    <dbl>    <dbl>   <int>     <int>    <int>
#> 1 BLANK     Groot   0.714      1.75     257.     266.       3         1        1
#> 2 Plant_A   Groot   0.968      2.33     281.     290.       1         1        1
#> 3 Plant_B   Groot   0.854      1.73     281.     290.       3         1        1
#> 4 Plant_C   Groot   0.982      2.11     309.     319.       1         1        1
#> 5 TMR       Groot   0.988      3.14     377.     386.       1         1        1
#> # ℹ 3 more variables: Rank_BIC <int>, Total_Rank <int>, Overall_Rank <int>

Model Win Frequency

model_win_frequency(
  best_models
)
#> # A tibble: 1 × 2
#>   Model Treatments_Won
#>   <chr>          <int>
#> 1 Groot              5

Quality Control Workflow

A typical workflow is:

Import data
    ↓
Validate metadata
    ↓
Process ANKOM data
    ↓
Validate processed data
    ↓
Fit models
    ↓
Flag problematic bottles
    ↓
Inspect residuals
    ↓
Exclude problematic bottles
    ↓
Refit models
    ↓
Compare models

Example bottle exclusion:

gp_clean <- exclude_heads(
  gp,
  heads = c("10"),
  reason = "Sensor malfunction"
)

Available Models

Current built-in models:

Brody
Dual Logistic
EXP0
EXPL
Gompertz
Groot
LE0
LEL
Logistic
Mitscherlich
Michaelis-Menten
Orskov and McDonald

Note on Groot and Michaelis-Menten

The Groot and generalized Michaelis-Menten models are mathematically equivalent.

Parameter correspondence:

VF = A
b = K
k = c

Researchers may choose either formulation depending on the terminology commonly used in their field.

Next Steps

Additional package capabilities include:

Importing Manual Datasets

gp <- as_rumen_gp(
  data = my_data,
  head_col = "Bottle",
  time_col = "Time",
  gas_col = "Gas"
)

Importing Pressure Data

gp <- as_rumen_gp(
  data = my_data,
  head_col = "Bottle",
  time_col = "Time",
  pressure_col = "PSI",
  pressure_unit = "psi",
  headspace_volume = 60
)

User-Defined Models

custom_fit <- fit_custom(
  data = gp,
  formula =
    Gas_mL ~
      A *
      (
        Time_h /
        (
          Time_h + K
        )
      ),
  start = list(
    A = 150,
    K = 10
  ),
  lower = c(
    A = 0,
    K = 0
  ),
  model_name = "Hyperbolic"
)

See:

?as_rumen_gp
?fit_custom

for additional details.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.