# Date: 2019-03-04 # internal version number 35 # clear existing variables rm(list = ls(all=TRUE)) library(readxl) library(tidyverse) # User data FileName <- c("inplant input data.xlsx") #FileName <- c("inplant input data000a.xlsx", # "inplant input dataALT1.xlsx", # "inplant input dataALT2a.xlsx", # "inplant input dataALT2b.xlsx", # "inplant input dataALT3.xlsx") ############################################################# ### Function: ### Purpose: ### Inputs: ### Returns: ### Modifies (global): ### Calls: ### Comment: ############################################################# ############################################################# ### Function: egr5c ### Purpose: generate exponential growth rate at 5C ### Comment: ############################################################# egr5c <- function(dm.id, frac_GIP) { # Generate n values of EGR at 5C using the parameters from the # FDA 2003 model. # Warning: since it is not an inactivation # model: values < 0 are set to 0 dmt <- as.vector(table(dm.id)) if(frac_GIP==1) dmtg <- c(0, dmt[1], 0, dmt[2], 0, dmt[3]) if(frac_GIP==0) dmtg <- c(dmt[1], 0, dmt[2], 0, dmt[3], 0) turkey.wo <- rGenerate(dmtg[1], egr.turkey.without.Data, empiric) turkey.w <- rGenerate(dmtg[2], egr.turkey.with.Data, empiric) ham.wo <- rGenerate(dmtg[3], egr.ham.without.Data, empiric) ham.w <- rGenerate(dmtg[4], egr.ham.with.Data, empiric) beef.wo <- rGenerate(dmtg[5], egr.beef.without.Data, empiric) beef.w <- rGenerate(dmtg[6], egr.beef.with.Data, empiric) res <- c(turkey.wo, turkey.w, ham.wo, ham.w, beef.wo, beef.w) res[res < 0] <- 0 return(res) } ############################################################# ### Function: egrT ### Purpose: adjust growth rate from 5C to input temperature ### Comment: ############################################################# egrT <- function(egr5, T){ egrT <- egr5/(6.18/(T+1.18))^2 return(egrT) } ############################################################# ### Function: F2C ### Purpose: Convert degrees farenheit to centigrade ### Comment: ############################################################# F2C <- function(F) (F-32) * 5/9 ############################################################# ### Function: fixtable ### Purpose: Match the row names of the table ### Comment: next version, use dplyr ############################################################# fixtable <- function(alt.ID, tbl) { # Forms the union of two tables s # The output has the same rows as alt.ID and in the same order # Create a matrix with 0 counts,the tmp <- matrix(rep.int(0,length(alt.ID)*2), nrow=length(alt.ID), ncol=2, dimnames = list(alt.ID, colnames(tbl))) rtbl <- tmp # union the two tables, with numbers coming from tbl first, # then filling in with the 0 counts if needed nt <- rbind(tbl, tmp[setdiff(rownames(tmp), rownames(tbl)),]) # reorder the table to match alt.ID row order for ( i in 1:length(alt.ID)) { idex <- match(alt.ID[i], row.names(nt)) rtbl[i,] <- nt[idex,] } return(rtbl) } ############################################################# ### Function: input.control ### Purpose: read in the control variables ### Comment: variables are returned as globals ############################################################# input.control <- function(inputFile) { controlVariables <- read_excel(inputFile, sheet="control") for(iVar in 1:length(controlVariables$LogicalVariables)) { Name <- controlVariables$LogicalVariables[iVar] if(is.na(Name) || Name == "") next assign(Name, as.logical(as.numeric(controlVariables$LogicValue[iVar])), envir = .GlobalEnv) } for(iVar in 1:length(controlVariables$NumericVariables)) { Name <- controlVariables$NumericVariables[iVar] if(is.na(Name) || Name == "") next assign(Name, controlVariables$NumValue[iVar], envir = .GlobalEnv) } for(iVar in 1:length(controlVariables$CharacterVariables)){ Name <- controlVariables$CharacterVariables[iVar] if(is.na(Name) || Name == "") next assign(Name, controlVariables$CharValue[iVar], envir = .GlobalEnv) } # following need to be in global environment LmOutputPerc <<- seq(0, 1, by=deltaProb) } ############################################################# ### Function: input.distributions ### Purpose: read in the distribution data for growth, lag time etc ### Comment: variables are returned as globals ############################################################# input.distributions <- function(inputFile) { dfDistributions <- read_excel(inputFile, sheet="distributions") dfDistr <- dfDistributions[1:25, 1:6] egr.turkey.with.Data <<- dfDistr[tolower(dfDistr$Parameter)=='egr 5c turkey with gi',] egr.turkey.without.Data <<- dfDistr[tolower(dfDistr$Parameter)=='egr 5c turkey without gi',] egr.ham.with.Data <<- dfDistr[tolower(dfDistr$Parameter)=='egr 5c ham with gi',] egr.ham.without.Data <<- dfDistr[tolower(dfDistr$Parameter)=='egr 5c ham without gi',] egr.beef.with.Data <<- dfDistr[tolower(dfDistr$Parameter)=='egr 5c beef with gi',] egr.beef.without.Data <<- dfDistr[tolower(dfDistr$Parameter)=='egr 5c beef without gi',] lagt.turkey.with.Data <<- dfDistr[tolower(dfDistr$Parameter)=='lag time 5c turkey with gi',] lagt.turkey.without.Data <<- dfDistr[tolower(dfDistr$Parameter)=='lag time 5c turkey without gi',] lagt.ham.with.Data <<- dfDistr[tolower(dfDistr$Parameter)=='lag time 5c ham with gi',] lagt.ham.without.Data <<- dfDistr[tolower(dfDistr$Parameter)=='lag time 5c ham without gi',] lagt.beef.with.Data <<- dfDistr[tolower(dfDistr$Parameter)=='lag time 5c beef with gi',] lagt.beef.without.Data <<- dfDistr[tolower(dfDistr$Parameter)=='lag time 5c beef without gi',] StorageTime.RetailSliced.Data <<- dfDistr[tolower(dfDistr$Parameter)=='home storage time, retail sliced',] StorageTime.Prepackaged.Data <<- dfDistr[tolower(dfDistr$Parameter)=='home storage time, prepackaged',] StorageTemperature.Data <<- dfDistr[tolower(dfDistr$Parameter)=='home storage temperature',] time.p2r.Data <<- dfDistr[tolower(dfDistr$Parameter)=='time, plant to retail',] temper.p2r.Data <<- dfDistr[tolower(dfDistr$Parameter)=='temperature, plant to retail',] TC.Data <<- dfDistr[tolower(dfDistr$Parameter)=='transfer coefficient',] LsLoad.Data <<- dfDistr[tolower(dfDistr$Parameter)=='lspp load during contamination',] LmLsRatio.Data <<- dfDistr[tolower(dfDistr$Parameter)=='lm lspp ratio',] FCSSwabArea.Data <<- dfDistr[tolower(dfDistr$Parameter)=='fcs swab area',] FCSArea.Data <<- dfDistr[tolower(dfDistr$Parameter)=='fcs area',] CEDurat.Data <<- dfDistr[tolower(dfDistr$Parameter)=='contamination event duration',] CEDelta.Data <<- dfDistr[tolower(dfDistr$Parameter)=='time between contamination events',] servingSize.Data <<- dfDistr[tolower(dfDistr$Parameter)=='serving size',] nEmpiric <<- sum(dfDistr$distribution=="empiric", na.rm=TRUE) if(nEmpiric >= 1) empiric <<- dfDistributions[, 9:(8 + 2 * nEmpiric)] } ############################################################# ### Function: input.dr ### Purpose: read in the dose response parameters ### Comment: variables are returned as globals ############################################################# input.dr <- function(inputFile) { drVariables <- read_excel(inputFile, sheet="dose response") for(iVar in 1:length(drVariables$Parameter)) { Name <- drVariables$Parameter[iVar] if(is.na(Name) || Name == "") next assign(Name, drVariables$Value[iVar], envir = .GlobalEnv) } nS_servings <<- frac.susceptible * nservings nH_servings <<- (1-frac.susceptible) * nservings r.healthy <<- c(r.healthy.q05, r.healthy.q50, r.healthy.q95) r.susceptible <<- c(r.susceptible.q05, r.susceptible.q50, r.susceptible.q95) } ############################################################# ### Function: input.parameters ### Purpose: read in the parameters ### Comment: variables are returned as globals ############################################################# input.parameters <- function(inputFile) { parameterVariables <- read_excel(inputFile, sheet="parameters") for(iVar in 1:length(parameterVariables$parameter)) { Name <- parameterVariables$parameter[iVar] if(is.na(Name) || Name == "") next assign(Name, parameterVariables$value[iVar], envir = .GlobalEnv) } # following need to be in global environment ReportLag <<- c(ReportLagLspp,ReportLagLm) sanitize <<- c(Sanitize.between.lots, Sanitize.between.days, Sanitize.enhanced) DMfrac <<- c(fracTurkey, fracHam, fracBeef) Lm.maxconc <<- c(Lm.maxconc.5, Lm.maxconc.5.7, Lm.maxconc.7) } ############################################################# ### Function: input.plants ### Purpose: read in the plant data ### Comment: variables are returned as globals ############################################################# input.plants <- function(inputFile) { dfPlants <- read_excel(inputFile, sheet="plants") # extract each column, convert to appropriate data type, and make global cat.ID <<- as.character(dfPlants$cat.ID[!is.na(dfPlants$cat.ID)]) # trim to just alternatives dfPlants <- dfPlants[1:length(cat.ID),] LotMassMean <<- as.numeric(dfPlants$LotMassMean) LotMassSD <<- as.numeric(dfPlants$LotMassSD) fracProduction <<- as.numeric(dfPlants$fracProduction) nFCSSamplesPerTP <<- as.numeric(dfPlants$nFCSSamplesPerTP) nLotSamplesPerTP <<- as.numeric(dfPlants$nLotSamplesPerTP) FCSSeqTrigger <<- as.numeric(dfPlants$FCSSeqTrigger) fracPP <<- as.numeric(dfPlants$fracPP) fracGIP <<- as.numeric(dfPlants$fracGI) Action.FCSEnhancedClean <<- as.logical(dfPlants$Action.FCSEnhancedClean) Action.FCSTestNextFCS <<- as.logical(dfPlants$Action.FCSTestNextFCS) Action.LotDispose <<- as.logical(dfPlants$Action.LotDispose) Action.LotTestNextLot <<- as.logical(dfPlants$Action.LotTestNextLot) Action.FCSForceLotTest <<- as.logical(dfPlants$Action.FCSForceLotTest) Action.TestAndHoldAll <<- as.logical(dfPlants$Action.TestAndHoldAll) } ############################################################# ### Function: lagtime5c ### Purpose: generate lag times at 5C ### Comment: since it is not an inactivation model: values < 0 are set to 0 ############################################################# lagtime5c <- function(dm.id, frac_GIP) { # Generate n values of EGR at 5C using the parameters # from the FDA 2003 model # Warning: since it is not an inactivation model: values < 0 are set to 0 dmt <- as.vector(table(dm.id)) if(frac_GIP==1) dmt <- c(0,dmt[1], 0, dmt[2], 0, dmt[3]) if(frac_GIP==0) dmt <- c(dmt[1], 0, dmt[2], 0, dmt[3], 0) n <- length(dm.id) res <- rep(NA,n) turkey.wo <- rGenerate(dmt[1], lagt.turkey.without.Data, empiric) turkey.w <- rGenerate(dmt[2], lagt.turkey.with.Data, empiric) ham.wo <- rGenerate(dmt[3], lagt.ham.without.Data, empiric) ham.w <- rGenerate(dmt[4], lagt.ham.with.Data, empiric) beef.wo <- rGenerate(dmt[5], lagt.beef.without.Data, empiric) beef.w <- rGenerate(dmt[6], lagt.beef.with.Data, empiric) res <- c(turkey.wo, turkey.w, ham.wo, ham.w, beef.wo, beef.w) return(res) } ############################################################# ### Function: lagtimeT ### Purpose: Convert the 5C lag time to match appropriate temperature ### Comment: ############################################################# lagtimeT <- function(lagtime5c, T) { lagtimeT <- lagtime5c/((T+1.18)/6.18)^2 } ############################################################# ### Function: LotSim ### Purpose: The main inplant model function ### Comment: ############################################################# LotSim <- function(run.ID, nlot, LotMassMean, LotMassSD, nFCSSamplesPerTP, nLotSamplesPerTP, Action.FCSEnhancedClean, Action.FCSTestNextFCS, Action.LotDispose, Action.LotTestNextLot, Action.FCSForceLotTest, Action.TestAndHoldAll, frac_PP, frac_GIP, Effect.PPMax, Effect.PPMin, CategoryID, NormalizeLotMass, FCSSeqTrigger) { nsim <- nlot+nLotsInit # Create/Initialize stochastic variables # The FCS Listeria species concentration, cfu/cm^2 FCSLs <- rep(0, nsim) # FCS test result, positive or negative FCSResult <- rep(FALSE, nsim) # Lm concentratio in the lot, cfu/gm LotLmPlant <- rep(0, nsim) # Lspp concentration in the lot, cfu/gm LotLs <- rep(0, nsim) PlantCategory<- rep(CategoryID, nsim) # Lm in lot after post processing lethality UsePP <- rbinom(nsim,1, frac_PP) # Lm concentration in the lot after PP, cfu/gm LotLmPP <- rep(0,nsim) LotLmPPResult <- rep(FALSE, nsim) InduceLot <- rep(FALSE, nsim) # Set the initial sanitization levels for each lots. # Assumes 2 lots per day # The +1 in case odd number nsim FCSSanitize <- rep(c(sanitize[1], sanitize[2]), (nsim+1) %/% 2) FCSSanitize <- FCSSanitize[1:nsim] # Should lot be tested because of FCS result LotTestbyFCS <- rep(FALSE, nsim) # Should lot be tested because of lagged previous lot positive LotTestbyLot <- rep(FALSE, nsim) # with different size plants, normalize the area # to the large plant size area FCSArea <- (LotMassMean/NormalizeLotMass) * rGenerate(nsim, FCSArea.Data, empiric) ### Generate the contamination event timing # Minor change from VB version that prevents adding extra day duration nDay <- ceiling(nsim / nLotsPerDay) CEdelta <- round(10^rGenerate(nDay, CEDelta.Data, empiric)) CEdurat <- round(10^rGenerate(nDay, CEDurat.Data, empiric)) CE <- rep(FALSE, nDay) iday <- 1 iCE <- 1 while (iday <= nDay) { iday <- iday + CEdelta[iCE] if (iday > nDay) break CE[iday] <- TRUE if (CEdurat[iCE] > 1) { for (i in 1:CEdurat[iCE]) CE[iday + i - 1] <- TRUE } iCE <- iCE+1 } CE <- rep(CE, each=nLotsPerDay) # Convert from day to lot timing CE <- CE[1:nsim] # Start Mass Balance # Set the FCS that will be tested BaseFCSTest <-SystematicTest(nsim, nDaysPerTP*nLotsPerDay, nFCSSamplesPerTP) FCSTest <- BaseFCSTest # Set the lots that will be tested LotTestSys <- SystematicTest(nsim, nDaysPerTP*nLotsPerDay, nLotSamplesPerTP) FCSSwabArea <- rGenerate(nsim, FCSSwabArea.Data, empiric) LotMass <- rnormtrunc(nsim, LotMassMean, LotMassSD, LotMassMin, LotMassMax) FCSLsAdded <- 10^rGenerate(nsim, LsLoad.Data, empiric) TC <- 10^rGenerate(nsim, TC.Data, empiric) TC <- ifelse(TC>1, 1, TC) LmLsRatio <- rGenerate(nsim, LmLsRatio.Data, empiric) LmLsRatio <- ifelse(LmLsRatio < 0, 0, LmLsRatio) LmLsRatio <- ifelse(LmLsRatio > 1, 1, LmLsRatio) # Consecutive positives as defined by USDA LsLag <- ReportLag[1] * nLotsPerDay LmLag <- ReportLag[2] * nLotsPerDay # Code for sequential positives here FCSSeqPos = 0 ## Test the lots & FCS; Listeria transfers. ## Take corrective actions if tested positive if (Action.TestAndHoldAll) { ## enable to test the corresponding lots ## after finding out the positive FCSs for(ilot in 2:nsim) { # Add Ls to FCS if ongoing contamination event FCSLs[ilot] <- FCSLs[ilot-1] + as.integer(CE[ilot]) * FCSLsAdded[ilot] # Keep FCS area constant if FCS Ls is not 0 if (FCSLs[ilot] != 0) { FCSArea[ilot] <- FCSArea[ilot-1] } # Test FCS for Ls FCSResult[ilot] <- TestPresence(1, FCSLs[ilot]*FCSSwabArea[ilot], ProbOf1) # Move the Ls from the FCS over to the lot. # The 454 converts pounds to grams # The test and hold process would not change the status # of FSCs and lots at the number of ilot LotLs[ilot] <- TC[ilot] * (FCSLs[ilot] * FCSArea[ilot] / (LotMass[ilot] * 454.0)) # Mass balance the remaining FCS Ls, both for transfer to lot and for sanitation FCSLs[ilot] <- FCSLs[ilot] * (1 - TC[ilot]) * (1-FCSSanitize[ilot]) # Don't track too low concentrations if (FCSLs[ilot] < MinFCSLs) { FCSLs[ilot] <- 0 } # Calculate the Lm in each lot at plant before any interventions LotLmPlant[ilot] <- LmLsRatio[ilot] * LotLs[ilot] LotLmPlant[ilot] <- ifelse(LotLmPlant[ilot] > Lm.maxconc[3], Lm.maxconc[3], LotLmPlant[ilot]) LotLmPP[ilot] <- ifelse(UsePP[ilot], ((1-runif(1,Effect.PPMin, Effect.PPMax)) * LotLmPlant[ilot]), LotLmPlant[ilot]) # Test each lot (although not all results actually known), return positive or negative LotLmPPResult[ilot] <- TestPresence(1, LotMassSampled * LotLmPP[ilot], ProbOf1) # If we actually tested FCS Ls and it was positive, take corrective actions if (FCSTest[ilot]) { # Prevent index for lot to be tested to fall beyond vector index because of the lag LotToTest <- min(ilot + LsLag,nsim) if (FCSResult[ilot]) { # Increment the sequential positive counter FCSSeqPos <- FCSSeqPos + 1 if (Action.FCSEnhancedClean) FCSSanitize[LotToTest:min((LotToTest+LsLag),nsim)] <- sanitize[3] # Test FCS immediately upon + news if (Action.FCSTestNextFCS) FCSTest[LotToTest] <- TRUE if (FCSSeqPos == FCSSeqTrigger) { # Set up additional FCS to be tested FCSTest[LotToTest:min(ilot+2*LsLag,nsim)] <- TRUE # Test current lot and the hold lot if required if (Action.FCSForceLotTest) { LotTestbyFCS[ilot] <- TRUE # record the tested lots at the trigger points, # i.e. the lots before held lots InduceLot[ilot] <- TRUE } ##Hold products [ilot+1:LotToTest+LsLag] HoldLotNumber <- min((ilot+1),nsim):min((LotToTest+LsLag),nsim) } if (FCSSeqPos > FCSSeqTrigger) { # Test additional FCS FCSTest[LotToTest:min(ilot+2*LsLag,nsim)] <- TRUE # Test lots being held if required if (Action.FCSForceLotTest) { LotTestbyFCS[HoldLotNumber] <- TRUE } } } else { # FCS result negative, so reset sequential counter and release held FCS tests FCSSeqPos <- 0 FCSTest[LotToTest:min(ilot+2*LsLag,nsim)] <- BaseFCSTest[LotToTest:min(ilot+2*LsLag,nsim)] } } } } else { ## according to "COMPLIANCE GUIDELINES TO CONTROL LISTERIA MONOCYTOGENES ## IN POST-LETHALITY EXPOSED READY-TO-EAT MEAT AND POULTRY PRODUCTS"(2006) ## http://www.fsis.usda.gov/oppde/rdad/FRPubs/97-013F/LM_Rule_Compliance_Guidelines_May_2006.pdf ## the FCSSeqTrigger should be 2 in this scenario for(ilot in 2:nsim) { # Add Ls to FCS if ongoing contamination event FCSLs[ilot] <- FCSLs[ilot-1] + as.integer(CE[ilot]) * FCSLsAdded[ilot] # Keep FCS area constant if FCS Ls is not 0 if (FCSLs[ilot] != 0) { FCSArea[ilot] <- FCSArea[ilot-1] } # Test FCS for Ls FCSResult[ilot] <- TestPresence(1, FCSLs[ilot]*FCSSwabArea[ilot], ProbOf1) # Move the Ls from the FCS over to the lot. The 454 converts pounds to grams # the test and hold process would not change the status of FSCs # and lots at the number of ilot LotLs[ilot] <- TC[ilot] * (FCSLs[ilot] * FCSArea[ilot] / (LotMass[ilot] * 454.0)) # Mass balance the remaining FCS Ls, both for transfer # to lot and for sanitation FCSLs[ilot] <- FCSLs[ilot] * (1 - TC[ilot]) * (1-FCSSanitize[ilot]) # Don't track too low concentrations if (FCSLs[ilot] < MinFCSLs) { FCSLs[ilot] <- 0 } # Calculate the Lm in each lot at plant before any interventions LotLmPlant[ilot] <- LmLsRatio[ilot] * LotLs[ilot] LotLmPlant[ilot] <- ifelse(LotLmPlant[ilot] > Lm.maxconc[3], Lm.maxconc[3], LotLmPlant[ilot]) LotLmPP[ilot] <- ifelse(UsePP[ilot], ((1-runif(1,Effect.PPMin, Effect.PPMax)) * LotLmPlant[ilot]), LotLmPlant[ilot]) # Test each lot (although not all results actually known), # return positive or negative LotLmPPResult[ilot] <- TestPresence(1, LotMassSampled * LotLmPP[ilot], ProbOf1) # If we actually tested FCS Ls and it was positive, take corrective actions if (FCSTest[ilot]) { # Prevent index for lot to be tested to fall beyond vector index # because of the lag LotToTest <- min(ilot + LsLag,nsim) if (FCSResult[ilot]) { # Increment the sequential positive counter FCSSeqPos <- FCSSeqPos + 1 # correct action on the following days if (Action.FCSEnhancedClean) FCSSanitize[LotToTest:min((LotToTest+LsLag),nsim)] <- sanitize[3] # Test FCS immediately upon + news if (Action.FCSTestNextFCS) FCSTest[LotToTest] <- TRUE if (FCSSeqPos == FCSSeqTrigger) { # Set up additional FCS to be tested FCSTest[LotToTest:min(ilot+2*LsLag,nsim)] <- TRUE # Test lot if required if (Action.FCSForceLotTest) { LotTestbyFCS[LotToTest] <- TRUE # record the tested lots at the trigger points, # i.e. the lots before held lots InduceLot[LotToTest] <- TRUE } ##Hold products [LotToTest+1:LotToTest+LsLag] HoldLotNumber <- min((LotToTest+1),nsim):min((LotToTest+LsLag),nsim) } if (FCSSeqPos > FCSSeqTrigger) { # Test additional FCS FCSTest[LotToTest:min(ilot+2*LsLag,nsim)] <- TRUE # Test lots being held if required if (Action.FCSForceLotTest) { # LotTestbyFCS[LotToTest:min(ilot+2*LsLag,nsim)] <- TRUE LotTestbyFCS[HoldLotNumber] <- TRUE } } } else { # FCS result negative, so reset sequential counter # and release held FCS tests FCSSeqPos <- 0 FCSTest[LotToTest:min(ilot+2*LsLag,nsim)] <- BaseFCSTest[LotToTest:min(ilot+2*LsLag,nsim)] } } } } ## for the test and hold process (recommonded scenario) ## if FCS tested negative after the trigger point and ## lot tested negative at the trigger point in the ## sequence (2 lots), release the hold products for (ilot in 1:nsim){ if (InduceLot[ilot]&!LotLmPPResult[ilot]&!FCSResult[ilot]) LotTestbyFCS[min(ilot+1,nsim):min(ilot+LsLag,nsim)] <- FALSE } ## test next lot if positive lots are found if (Action.LotTestNextLot==TRUE) { for (ilot in 1:nsim) { LotTestbyLot[min(ilot+LmLag,nsim)] <- (LotTestSys[ilot] | LotTestbyFCS[ilot]) & LotLmPPResult[ilot] } } #################################################################### ################### Lm growth from plants to consumer ############## #################################################################### # Determine which lots were sliced at retail, # 1 for retail sliced, 0 for prepackaged. RetailSliced <- double(nsim) RetailSliced <- rbinom(nsim, 1, FracRetailSlice) # Categorize the products dm.type <- sample(1:3, nsim, replace=TRUE, prob=DMfrac) UseGIP <- rbinom(nsim, 1, frac_GIP) dm.id <- dm.type * 2 + UseGIP - 1 # Calculate the egr5 for each serving. Depends on the # GIP use and deli meat type EGR5 <- egr5c(dm.id, frac_GIP) # deli meat id ############ cat("run.ID is ", run.ID,"\n") cat(" fracGIP is ", frac_GIP,"\n") cat(" fracPP is ", frac_PP,"\n") ############ # Calculate the time and temperature from plant to retail time.plant2retail <- rGenerate(nsim, time.p2r.Data, empiric) temper.plant2retail <- rGenerate(nsim, temper.p2r.Data, empiric) lagt <- lagtimeT(lagtime5c(dm.id,frac_GIP), temper.plant2retail) # Calculate the consumer storage time for this serving # based on where it was sliced # Serving time is a weibull variable, and the parameters # are calculated from the Regis data storageTime.Retail <- rGenerate(nsim, StorageTime.RetailSliced.Data, empiric) storageTime.Prepackaged <- rGenerate(nsim, StorageTime.Prepackaged.Data, empiric) serving.StorageTime <- ifelse(RetailSliced, storageTime.Retail, storageTime.Prepackaged) # Calculate total egr for time beyond lag time EGR.plant2retail <- double(nsim) EGR.plant2retail <- egrT(EGR5, temper.plant2retail) * pmax(0,(time.plant2retail-lagt)) # LotLmRetailBefore is the Lm entering retail. LotLmRetailBefore <- LotLmPP * 10^EGR.plant2retail LotLmRetailBefore <- ifelse(LotLmRetailBefore > Lm.maxconc[3], Lm.maxconc[3], LotLmRetailBefore) # After retail slicing # calculate the increase of Lm if retailsliced LotLmRetailSlicing <- ifelse(RetailSliced, 10^(mu.RetailSliced - sigma.RetailSliced*mu.Prepackaged/sigma.Prepackaged) * (LotLmRetailBefore^(sigma.RetailSliced/sigma.Prepackaged)), LotLmRetailBefore) LotLmRetailSlicing <- ifelse(LotLmRetailSlicing > Lm.maxconc[3], Lm.maxconc[3], LotLmRetailSlicing) # Test each lot at retail, although not all results known LotLmRetailAllResult <- TestPresence(nsim, LotMassSampled*LotLmRetailSlicing, ProbOf1) # Detmine which lots were tested either by systematic lot testing or # because of FCS positives and returned positive LotUse <- ifelse(rep(Action.LotDispose,nsim), !((LotTestSys | LotTestbyFCS | LotTestbyLot) & LotLmPPResult), rep(TRUE,nsim)) # If not used, set the Lm concentration to NA, i.e. removed from food supply LotLmRetailUsed <-ifelse(LotUse, LotLmRetailSlicing, NA) # Calculate test result for lots in food supply. Use False for any # removed, actual result for others not tested or did pass # note that this result is not available in real world unless all lots tested. LotLmRetailUsedResult <- ifelse(LotUse, LotLmRetailAllResult, NA) ## from retail to consumption ## time available for growth, i.e. time after lag time time.consumer <- ifelse(lagt-time.plant2retail > 0, pmax(0, serving.StorageTime-(lagt-time.plant2retail)), serving.StorageTime) serving.StorageTemperature <- double(nsim) # Adjust the values so that no storage temperature is below zero serving.StorageTemperature <- pmax(F2C(rGenerate(nsim, StorageTemperature.Data, empiric)),0) EGR.consumer <- double(nsim) # nvariab is the number of lots EGR.consumer <- egrT(EGR5, serving.StorageTemperature) * time.consumer ## Don't track too low lm (smaller than 1 cfu per lot) LotLmRetailUsed[LotLmRetailUsed < 10^(-8)] <- 0 consumption.lm <- LotLmRetailUsed * 10^EGR.consumer ######################################################################### ######################## Dose response model ############################ ######################################################################### # The following few lines check the concentration and temperature, # to see if they exceed the calculated maximum concentration at that # temperature. If they do, it replaces them with the maximum. lm.temp <- consumption.lm lm.temp <- ifelse(serving.StorageTemperature <= 5 & lm.temp > Lm.maxconc[1], Lm.maxconc[1], lm.temp) lm.temp <- ifelse(serving.StorageTemperature > 5 & serving.StorageTemperature <= 7 & lm.temp > Lm.maxconc[2], Lm.maxconc[2], lm.temp) lm.temp <- ifelse(serving.StorageTemperature > 7 & lm.temp > Lm.maxconc[3], Lm.maxconc[3], lm.temp) lm.temp[lm.temp >= max.lm.consumed] <- 0 lm.temp[is.na(lm.temp)] <- 0 conc.consumption.lm <- lm.temp # To get serving size, FDA-FSIS 2003 uses uniform random # number generator, then does linear interpolation # Use linear approximation to calculate serving sizes. # Generate a uniform variable between # Zero and one, then find out where that would be on # the quantile vector of serving # Sizes if you assumed the line was linear between them, # from FDA-FSIS 2003 serving.size <- rGenerate(nsim, servingSize.Data, empiric) # Convert the Listeria Concentration (consumption.lm is in CFU per gram) # to an actual CFU count by multiplying it by the size of # the serving, in grams. dose.lm <- conc.consumption.lm * serving.size[1:length(conc.consumption.lm)] # Calculate the chance of someone getting sick if they ate a dose. # Done for both the healthy and the susceptible # Exponential since dose.lm is a concentration (and not a count # of cfu. Otherwize: 1-(1-r)^n) # Mean risk for each population group risk.per.serving.H <- 1-exp(-r.healthy[2] * dose.lm) risk.per.serving.Hlo <- 1-exp(-r.healthy[1] * dose.lm) risk.per.serving.Hhi <- 1-exp(-r.healthy[3] * dose.lm) risk.per.serving.S <- 1-exp(-r.susceptible[2] * dose.lm) risk.per.serving.Slo <- 1-exp(-r.susceptible[1] * dose.lm) risk.per.serving.Shi <- 1-exp(-r.susceptible[3] * dose.lm) risk.per.serving.All <- risk.per.serving.H * (1-frac.susceptible) + risk.per.serving.S * frac.susceptible ####################################################################### ## write the output to a data frame ####################################################################### LotData <- data.frame(run.ID, dm.id, PlantCategory, UseGIP, UsePP, lagt, FCSArea, FCSTest, CE, FCSLsAdded, FCSLs, FCSResult, LotMass, LotLs, LotLmPlant, LotLmPP, LotLmPPResult, LotTestbyFCS, LotTestSys, LotTestbyLot, LotUse, EGR.plant2retail, LotLmRetailBefore, RetailSliced, LotLmRetailSlicing, LotLmRetailUsed, LotLmRetailUsedResult, serving.StorageTemperature, serving.StorageTime, time.consumer, EGR.consumer, serving.size, conc.consumption.lm, dose.lm, risk.per.serving.All, risk.per.serving.H, risk.per.serving.S) # Remove the initializing lots LotData <- LotData[-(1:nLotsInit),] # individual alternatives returned as global # trim initializing lots then calc mean risk for each population group risk.per.serving.H.trim <- risk.per.serving.H[-(1:nLotsInit)] risk.per.serving.Hlo.trim <- risk.per.serving.Hlo[-(1:nLotsInit)] risk.per.serving.Hhi.trim <- risk.per.serving.Hhi[-(1:nLotsInit)] risk.per.serving.S.trim <- risk.per.serving.S[-(1:nLotsInit)] risk.per.serving.Slo.trim <- risk.per.serving.Slo[-(1:nLotsInit)] risk.per.serving.Shi.trim <- risk.per.serving.Shi[-(1:nLotsInit)] riskH[run.ID] <<- Risk.Healthy <- mean(risk.per.serving.H.trim) riskHlo[run.ID] <<- Risk.HealthyLo <- mean(risk.per.serving.Hlo.trim) riskHhi[run.ID] <<- Risk.HealthyHi <- mean(risk.per.serving.Hhi.trim) riskS[run.ID] <<- Risk.Susceptible <- mean(risk.per.serving.S.trim) riskSlo[run.ID] <<- Risk.SusceptibleLo <- mean(risk.per.serving.Slo.trim) riskShi[run.ID] <<- Risk.SusceptibleHi <- mean(risk.per.serving.Shi.trim) return(LotData) } ############################################################# ### Function: LotsCSV ### Purpose: Assigns better names and writes all lots to csv ### Calls: ############################################################# LotsCSV <- function(LD) { LotDatatoCSV <- data.frame("runID"= LD$run.ID, "DeliMeatID" = LD$dm.id, "PlantCategory" = LD$PlantCategory, "UseGrowthInhibitor" = LD$UseGIP, "UsePostProcessingLethality" = LD$UsePP, "LagTimePriorToGrowth_d" = LD$lagt, "FCSArea" = LD$FCSArea, "FCSTested" = LD$FCSTest, "ContaminationEvent" = LD$CE, "LsAddedFCS" = LD$FCSLsAdded, "FCSLs" = LD$FCSLs, "FCSTestResult" = LD$FCSResult, "LotMass_g" = LD$LotMass, "LotLs_cfug" = LD$LotLs, "LotLmInPlant_cfug" = LD$LotLmPlant, "LotLmPostProc_cfug" = LD$LotLmPP, "LotLmPostProcResult" = LD$LotLmPPResult, "LotTestbyFCS" = LD$LotTestbyFCS, "LotTestSystematic" = LD$LotTestSys, "LotTestbyLot" = LD$LotTestbyLot, "LotUsed" = LD$LotUse, "ExpGrowth_PlantToRetail" = LD$EGR.plant2retail, "LotLmEnteringRetail_cfug" = LD$LotLmRetailBefore, "RetailSliced" = LD$RetailSliced, "LotLmAfterSlicing_cfug" = LD$LotLmRetailSlicing, "LotLmRetailIfUsed" = LD$LotLmRetailUsed, "LotLmRetailResultIfUsed" = LD$LotLmRetailUsedResult, "ConsumerStorageTemperature_C" = LD$serving.StorageTemperature, "ConsumerStorageTime_d" = LD$serving.StorageTime, "ConsumerTimeAfterLag_d" = LD$time.consumer, "ExpGrowth_Consumer" = LD$EGR.consumer, "ServingSize_g" = LD$serving.size, "ConcLmConsumption_cfug" = LD$conc.consumption.lm, "DoseLm_cfu" = LD$dose.lm, "RiskPerServing" = LD$risk.per.serving.All, "RiskPerServingHealthy" = LD$risk.per.serving.H, "RiskPerServingSusceptible" = LD$risk.per.serving.S) write.csv(LotDatatoCSV, file="individual lot results.csv", row.names = TRUE) } ############################################################# ### Function: PlotFractions ### Purpose: Plot the fractions of different type of Lots ### Comment: ############################################################# PlotFractions <- function(fracs, fracTitles, disk) { if (disk) { png(file="fractions.png", height=600, width=600) } barplot(fracs, names.arg=fracTitles, axis.lty=1, ylab="Fraction of lots at Retail", col="blue") if (disk) { dev.off() } par(mfrow=c(1,1), cex=1, cex.axis=1, cex.lab=1, col="black") } ############################################################# ### Function: PlotPrevalences ### Purpose: Create plot of the prevalences of Lm at different stage ### Comment: ############################################################# PlotPrevalences <- function(prevalence, prevalence.names, disk) { if (disk) { png("prevalences.png", height=600, width=600) } barplot(prevalence, names.arg=prevalence.names, axis.lty=1, ylab="Prevalence", col="red") if (disk) { dev.off() } par(mfrow=c(1,1), cex=1, cex.axis=1, cex.lab=1, col="black") } ############################################################# ### Function: PlotQuantileComparison ### Purpose: Create plot of quantile of Lm for various interventions ### Comment: ############################################################# PlotQuantileComparison <- function(run.no, ypct, c1, c2, c3, c4, c5, disk) { if (disk) { png(file="quantiles lmra.png", height=600, width=600) } par(mfrow=c(1,1), cex=1.5, cex.axis=1.2, cex.lab=1.2, lwd=2) t1<-log10(ifelse(c1 <= 0, 1e-20, c1)) t2<-log10(ifelse(c2 <= 0, 1e-20, c2)) t3<-log10(ifelse(c3 <= 0, 1e-20, c3)) t4<-log10(ifelse(c4 <= 0, 1e-20, c4)) t5<-log10(ifelse(c5 <= 0, 1e-20, c5)) cmax<-ceiling(max(t1,t2,t3,t4,t5,na.rm=TRUE)) cmin<-floor(min(t1,t2,t3,t4,t5,na.rm=TRUE)) plot(t1, ypct, col="black", type="l", lty=1, xlab="log Lm at retail (cfu/g)", ylab="Cumulative Frequency", xlim=c(cmin,cmax), ylim=c(0.8,1), col.axis="black") lines(t2, ypct, lty=2, col="green") lines(t3, ypct, lty=3, col="blue") lines(t4, ypct, lty=4, col="orange") lines(t5, ypct, lty=5, col="red") p.legend <- c("Combined", "PP&GI", "PP Only", "GI Only", "Neither") par(cex.lab=0.8) legend("bottomright", p.legend, lty=c(1,2,3,4,5), cex=0.75, col=c("black", "green", "blue", "orange", "red" ), text.col="black") if (disk) { dev.off() } par(mfrow=c(1,1), cex=1, cex.axis=1, cex.lab=1, col="black") } ############################################################# ### Function: PlotQuantileComparisonBaseline ### Purpose: for just looking at baseline (no intervention) results ### Comment: ############################################################# PlotQuantileComparisonBaseline <- function(run.no, ypct, lpi1, lri1, c1, disk) { if (disk) { png(file="quantiles lmra baseline no interventions.png", height=600, width=600) } par(mfrow=c(1,1), cex=1.5, cex.axis=1.2, cex.lab=1.2, lwd=2) t1<-log10(ifelse(lpi1 <= 0, 1e-20, lpi1)) t2<-log10(ifelse(lri1 <= 0, 1e-20, lri1)) t3<-log10(ifelse(c1 <= 0, 1e-20, c1)) cmax<-ceiling(max(t1,t2,t3,na.rm=TRUE)) cmin<-floor(min(t1,t2,t3,na.rm=TRUE)) plot(t1, ypct, col="black", type="l", lty=1, xlab="log Lm (cfu/g)", ylab="Cumulative Frequency", xlim=c(-10,cmax), ylim=c(0.75,1), col.axis="black") lines(t2, ypct, lty=2, col="green") lines(t3, ypct, lty=3, col="blue") p.legend <- c("End of Production", "End of Retail", "At Consumption") par(cex.lab=0.8) legend("bottomright", p.legend, lty=c(1,2,3), cex=0.75, col=c("black", "green", "blue" ), text.col="black") if (disk) { dev.off() } par(mfrow=c(1,1), cex=1, cex.axis=1, cex.lab=1, col="black") } ############################################################# ### Function: PlotSampleTimeSeries ### Purpose: Create the sample time series plot ### Comment: ############################################################# PlotSampleTimeSeries <-function(Lots, disk, startLot=6000, endLot=6500) { if (disk) { png(file="sample time series.png", height=600, width=600) } t1 <- log10(ifelse( Lots[,"LotLmRetailUsed"]<=0, 1e-15, Lots[,"LotLmRetailUsed"])) t2 <- log10(ifelse( Lots[,"FCSLs"]<=0, 1e-15, Lots[,"FCSLs"])) t3 <- Lots[,"FCSResult"] x <- startLot:endLot par(mfrow=c(3,1), mex=0.9) plot(x, t2[x], xlab="Time (lot)", ylab="FCS Lspp (cfu/cm^2)", col="black") plot(x, t3[x], xlab="Time (lot)", ylab="FCS Test Result", col="red", ylim=c(0,1)) plot(x, t1[x], xlab="Time (lot)", ylab="Lm Retail (cfu/g)", col="blue") par(mfrow=c(1,1)) if (disk) { dev.off() } } ############################################################# ### Function: PlotTests ### Purpose: Create plots of test results in each plant ### Comment: ############################################################# PlotTests <- function(plantID, nLots, nFCSTests, fFCSPos, nLotTests, fLotPos, disk) { if (disk) { png(file="tests.png", height=800, width=600) } par(mfrow=c(5,1)) if (max(is.na(fLotPos))==1 & min(is.na(fLotPos)) == 1) { barplot(fLotPos, ylim=c(0,0.1), names.arg=plantID, axis.lty=1, ylab="Fraction Lots Positive") } else { barplot(fLotPos, names.arg=plantID, axis.lty=1, ylab="Fraction Lots Positive") } barplot(nLotTests, names.arg=plantID, axis.lty=1, ylab="No. Lots Tested") if (max(is.na(fFCSPos))==1 & min(is.na(fFCSPos)) == 1) { barplot(fFCSPos, ylim=c(0,0.1), names.arg=plantID, axis.lty=1, ylab="Fraction FCS Positive") } else { barplot(fFCSPos, names.arg=plantID, axis.lty=1, ylab="Fraction FCS Positive") } barplot(nFCSTests, names.arg=plantID, axis.lty=1, ylab="No. FCS Tested") barplot(nLots, names.arg=plantID, axis.lty=1, ylab="No. Lots Produced") if (disk) { dev.off() } par(mfrow=c(1,1), cex=1, cex.axis=1, cex.lab=1, col="black") } ############################################################# ### Function: rGenerate ### Purpose: generate variability draws based on ### the specified distribution ### Comment: requires empiric to be global ############################################################# rGenerate <- function (n, parlist, empir) { # requires empiric to be global if (tolower(parlist$distribution) == "binomial") rbinom(n, 1, parlist$parameter1) else if (tolower(parlist$distribution) == "empiric") { probName <- paste(parlist$empiricName, "Prob", sep="") probEmpir <- empir[,probName][[1]] probEmpir <- probEmpir[!is.na(probEmpir)] probQ <- empir[,parlist$empiricName][[1]] probQ <- probQ[!is.na(probQ)] approx(probEmpir, probQ, xout=runif(n))$y } else if (tolower(parlist$distribution) == "exponential") rexp(n, parlist$parameter1) else if (tolower(parlist$distribution) == "fixed") rep(parlist$parameter1, times=n) else if (tolower(parlist$distribution) == "laplace") rlaplace(n, location=parlist$parameter1, scale=parlist$parameter2) else if (tolower(parlist$distribution) == "logistic") rlogis(n, location=parlist$parameter1, scale=parlist$parameter2) else if (tolower(parlist$distribution) == "lognormal") 10^(rnorm(n, mean=parlist$parameter1, sd=parlist$parameter2)) else if (tolower(parlist$distribution) == "normal") rnorm(n, mean=parlist$parameter1, sd=parlist$parameter2) else if (tolower(parlist$distribution) == "triangular") rTriangle(n, lo=parlist$parameter1, mid=parlist$parameter2, hi=parlist$parameter3) else if (tolower(parlist$distribution) == "uniform") runif(n, min=parlist$parameter1, max=parlist$parameter2) else if (tolower(parlist$distribution) == "weibull") rweibull(n, shape=parlist$parameter1, scale=parlist$parameter2) else warning(c("Error in rGenerate: distribution not found ", distribution), immediate.=TRUE) } ############################################################# ### Function: rTriangle ### Purpose: generate random draws from a triangular distribution ### Comment: ############################################################# rTriangle <- function(n, lo, mid, hi) { if(any(lo > mid) || any(mid > hi)) stop("Need lo < mid < hi in rTriangle") lo <- rep(lo, length=n) hi <- rep(hi, length=n) ifelse(lo==hi, lo, { cr <- (mid-lo)/(hi-lo) u <- runif(n,0,1) x <- ifelse(u<=cr, (cr*u)^0.5, 1-((1-cr)*(1-u))^ 0.5) lo+(hi-lo)*x }) } ############################################################# ### Function: rnormtrunc ### Purpose:Truncate the normal distribution. Replace with 0,1 bounds ### Comment: ############################################################# rnormtrunc <- function(n, xMean, xSD, xMin, xMax) { x<-rnorm(n, xMean, xSD) x<-ifelse(x < xMin, xMin, x) x<-ifelse(x > xMax, xMax, x) return(x) } ############################################################# ### Function: SystematicTest ### Purpose: Setup the systematic test for FCS and Lots ### Comment: ############################################################# SystematicTest <- function(n, nPossiblePerTP, nSamplesPerTP) { if (nSamplesPerTP == 0) { TestOn <-rep(FALSE, n) } else if (nSamplesPerTP < 1 && nSamplesPerTP > 0) { nPossiblePerTP <- 12 * nPossiblePerTP nSamplesPerTP <- 12 * nSamplesPerTP StartRange <- nPossiblePerTP %/% nSamplesPerTP StartTest <- sample(1:StartRange, 1) TestTP <- seq(StartTest, nPossiblePerTP, StartRange) mtest <- rep(FALSE, nPossiblePerTP) for (i in 1:nSamplesPerTP) { mtest[TestTP[i]] <- TRUE } TestOn <- rep(mtest, (n %/% nPossiblePerTP) + 1) TestOn <- TestOn[-((n+1):length(TestOn))] } else { StartRange <- nPossiblePerTP %/% nSamplesPerTP StartTest <- sample(1:StartRange,1) TestTP <- seq(StartTest, nPossiblePerTP, StartRange) mtest <- rep(FALSE, nPossiblePerTP) for (i in 1:nSamplesPerTP) { mtest[TestTP[i]] <- TRUE } TestOn <- rep(mtest, (n %/% nPossiblePerTP) + 1) TestOn <- TestOn[-((n+1):length(TestOn))] } return(TestOn) } ############################################################# ### Function: TestPresence ### Purpose: test the presence of Lss and Lm ### Comment: ############################################################# TestPresence <- function( n, MeanConc, ProbOf1) { MeanConc <- pmin(MeanConc, 1e8) nCFU <- rpois(n, MeanConc) TestResult <- ifelse(runif(n) <( (1-ProbOf1)^nCFU), FALSE, TRUE) # False negatives with the following line # sum(!TestResult & (nCFU>0)) } ############################################################# ### Function: trimSpace ### Purpose: trimSpace function taken from seqinr library. ### Removes leading and/or trailing spaces from text string ### Comment: ############################################################# trimSpace <- function (x, leading = TRUE, trailing = TRUE, space = "[:space:]") { if (leading) { pattern <- paste("^[", space, "]*", sep = "", collapse = "") x <- sub(pattern = pattern, replacement = "", x = x) } if (trailing) { pattern <- paste("[", space, "]*$", sep = "", collapse = "") x <- sub(pattern = pattern, replacement = "", x = x) } return(x) } ############################################################################# ############################################################################# ### start actual program ### ### only include alternatives that are actually produced ### ### for each category, determine number of lots to be used and run model. ### ### Combine all output into 1 data frame ### ############################################################################# ############################################################################# wd <- getwd() for (k in 1:length(FileName)){ setwd(wd) inputFile <- FileName[k] input.control(inputFile) input.parameters(inputFile) input.dr(inputFile) input.plants(inputFile) input.distributions(inputFile) folderName <- runName # paste("outputs ", runName, sep="") dir.create(folderName, showWarnings = FALSE) wdir <- paste("./", folderName, sep="") setwd(wdir) set.seed(seedRN) j <- 1 run.no <- j # number of alternatives nalt <- length(cat.ID) frac <- fracProduction nlot.total <- nlot.total nLotSim <- round(frac*nlot.total, digits=0) # define the risk (illnesses per servings) of Lm in Healthy and Susceptible population for each alternative riskH <- rep(0,nalt) riskHlo <- rep(0,nalt) riskHhi <- rep(0,nalt) riskS <- rep(0,nalt) riskSlo <- rep(0,nalt) riskShi <- rep(0,nalt) riskAll <- rep(0,nalt) riskAllLo <- rep(0,nalt) riskAllHi <- rep(0,nalt) # record the beginning time ptm <- proc.time() NormalizeLotMass <- max(LotMassMean) for (icat in 1:nalt) { x <- LotSim(icat, nLotSim[icat], LotMassMean[icat], LotMassSD[icat], nFCSSamplesPerTP[icat], nLotSamplesPerTP[icat], Action.FCSEnhancedClean[icat], Action.FCSTestNextFCS[icat], Action.LotDispose[icat], Action.LotTestNextLot[icat], Action.FCSForceLotTest[icat],Action.TestAndHoldAll[icat], fracPP[icat], fracGIP[icat], Effect.PPMax, Effect.PPMin, cat.ID[icat], NormalizeLotMass, FCSSeqTrigger[icat]) if (icat == 1) { Lots <- x } else { Lots <- rbind(Lots, x) } } print(proc.time() - ptm) ### optionally save all outputs to disk, can be large if (SaveDetail) LotsCSV(Lots) # calculate and output summary results for each plant/alternative frac.sim <- as.vector(tapply(Lots[,"LotMass"], Lots[,"PlantCategory"], sum)/sum(Lots[,"LotMass"])) altSummary <- cbind(cat.ID, frac, LotMassMean, LotMassSD, nFCSSamplesPerTP, nLotSamplesPerTP, fracPP, fracGIP, FCSSeqTrigger, Action.FCSEnhancedClean, Action.FCSTestNextFCS, Action.LotDispose, Action.LotTestNextLot, Action.FCSForceLotTest, nLotSim, frac.sim, riskH, riskS) write.csv(altSummary, file="summary results by alternative.csv", row.names = FALSE) # calculate and output summary risk statistics # the risk stats are calculated globally in the LotSim function # calculate the annual illnesses caused by Lm riskHealthy <- sum(riskH*fracProduction) riskHealthyLo <- sum(riskHlo* fracProduction) riskHealthyHi <- sum(riskHhi* fracProduction) riskSusceptible <- sum(riskS* fracProduction) riskSusceptibleLo <- sum(riskSlo* fracProduction) riskSusceptibleHi <- sum(riskShi* fracProduction) AnnualRiskH <- riskHealthy*nH_servings AnnualRiskHlo <- riskHealthyLo*nH_servings AnnualRiskHhi <- riskHealthyHi*nH_servings AnnualRiskS <- riskSusceptible*nS_servings AnnualRiskSlo <- riskSusceptibleLo*nS_servings AnnualRiskShi <- riskSusceptibleHi*nS_servings AnnualRiskAll <- AnnualRiskH + AnnualRiskS AnnualRiskAllLo <- AnnualRiskHlo + AnnualRiskSlo AnnualRiskAllHi <- AnnualRiskHhi + AnnualRiskShi AnnRiskAll <- AnnualRiskAll SummaryRisk <- c("Sum of production fractions" = sum(frac), "No. lots simulated" = length(Lots[,1]), "risk per serving for healthy population" = riskHealthy, "risk per serving for healthy population, lower CI" = riskHealthyLo, "risk per serving for healthy population, upper CI" = riskHealthyHi, "risk per serving for susceptible" = riskSusceptible, "risk per serving for susceptible, lower CI" = riskSusceptibleLo, "risk per serving for susceptible, upper CI" = riskSusceptibleHi, "illnesses per year, healthy" = AnnualRiskH, "illnesses per year, healthy, lower CI" = AnnualRiskHlo, "illnesses per year, healthy, upper CI" = AnnualRiskHhi, "illnesses per year, susceptible" = AnnualRiskS, "illnesses per year, susceptible, lower CI" = AnnualRiskSlo, "illnesses per year, susceptible, upper CI" = AnnualRiskShi, "number of illnesses per year" = AnnualRiskAll, "number of illnesses per year, lower CI" = AnnualRiskAllLo, "number of illnnesses per year, upper CI" = AnnualRiskAllHi) write.csv(SummaryRisk, file="summary risk results.csv", row.names = TRUE) # calculate and output quantiles for each intervention at plant, retail, and consumption # SWITCHED TO LOTS USED LU <- Lots[Lots[,"LotUse"]==TRUE,] # calculate average and standard deviation of concentrations at plant, retail and consumption Lavg_proc <- mean(Lots[Lots$PlantCategory=="1-L","LotLmPP"]) Savg_proc <- mean(Lots[Lots$PlantCategory=="1-S","LotLmPP"]) VSavg_proc <- mean(Lots[Lots$PlantCategory=="1-VS","LotLmPP"]) Lsd_proc <- sd(Lots[Lots$PlantCategory=="1-L","LotLmPP"]) Ssd_proc <- sd(Lots[Lots$PlantCategory=="1-S","LotLmPP"]) VSsd_proc <- sd(Lots[Lots$PlantCategory=="1-VS","LotLmPP"]) Lavg_ret <- mean(Lots[Lots$PlantCategory=="1-L" ,"LotLmRetailUsed"],na.rm=TRUE) Savg_ret <- mean(Lots[Lots$PlantCategory=="1-S","LotLmRetailUsed"],na.rm=TRUE) VSavg_ret <- mean(Lots[Lots$PlantCategory=="1-VS","LotLmRetailUsed"],na.rm=TRUE) Lsd_ret <- sd(Lots[Lots$PlantCategory=="1-L","LotLmRetailUsed"],na.rm=TRUE) Ssd_ret <- sd(Lots[Lots$PlantCategory=="1-S","LotLmRetailUsed"],na.rm=TRUE) VSsd_ret <- sd(Lots[Lots$PlantCategory=="1-VS","LotLmRetailUsed"],na.rm=TRUE) Lavg_cons <- mean(Lots[Lots$PlantCategory=="1-L","conc.consumption.lm"]) Savg_cons <- mean(Lots[Lots$PlantCategory=="1-S","conc.consumption.lm"]) VSavg_cons <- mean(Lots[Lots$PlantCategory=="1-VS","conc.consumption.lm"]) Lsd_cons <- sd(Lots[Lots$PlantCategory=="1-L","conc.consumption.lm"]) Ssd_cons <- sd(Lots[Lots$PlantCategory=="1-S","conc.consumption.lm"]) VSsd_cons <- sd(Lots[Lots$PlantCategory=="1-VS","conc.consumption.lm"]) leavingPlantConc <- cbind("leaving plant", Lavg_proc, Lsd_proc, Savg_proc, Ssd_proc, VSavg_proc, VSsd_proc) atRetailConc <- cbind("leaving retail", Lavg_ret, Lsd_ret, Savg_ret, Ssd_ret, VSavg_ret, VSsd_ret) atConsumptionConc <-cbind("at consumption", Lavg_cons, Lsd_cons, Savg_cons, Ssd_cons, VSavg_cons, VSsd_cons) dfConcAvg <- rbind(leavingPlantConc,atRetailConc,atConsumptionConc) colnames(dfConcAvg) <- c("location", "L avg", "L sd", "S avg", "S sd", "VS avg", "VS sd") write.csv(dfConcAvg, file="summary conc stats by plant category.csv", row.names = FALSE) } # output the quantiles of each intervention leaving plant lpi1 <- quantile(LU[,"LotLmPP"], LmOutputPerc, names=FALSE, na.rm=TRUE) lpi2 <- quantile(LU[LU[,"UsePP"]==TRUE & LU[,"UseGIP"]==TRUE, "LotLmPP"], LmOutputPerc, names=FALSE, na.rm=TRUE) lpi3 <- quantile(LU[LU[,"UsePP"]==TRUE & LU[,"UseGIP"]==FALSE, "LotLmPP"], LmOutputPerc, names=FALSE, na.rm=TRUE) lpi4 <- quantile(LU[LU[,"UsePP"]==FALSE & LU[,"UseGIP"]==TRUE, "LotLmPP"], LmOutputPerc, names=FALSE, na.rm=TRUE) lpi5 <- quantile(LU[LU[,"UsePP"]==FALSE & LU[,"UseGIP"]==FALSE, "LotLmPP"], LmOutputPerc, names=FALSE, na.rm=TRUE) lpi6 <- quantile(LU[LU[,"UseGIP"]==FALSE,"LotLmPP"], LmOutputPerc, names=FALSE, na.rm=TRUE) lpi7 <- quantile(LU[LU[,"UseGIP"]==TRUE, "LotLmPP"], LmOutputPerc, names=FALSE, na.rm=TRUE) qLeavingPlantByIntervenction <- cbind(rep("leaving plant", length(LmOutputPerc)), LmOutputPerc, lpi1, lpi2, lpi3, lpi4, lpi5, lpi6, lpi7) # output the quantiles of each intervention leaving retail lri1 <- quantile(LU[,"LotLmRetailUsed"], LmOutputPerc, names=FALSE, na.rm=TRUE) lri2 <- quantile(LU[LU[,"UsePP"]==TRUE & LU[,"UseGIP"]==TRUE, "LotLmRetailUsed"], LmOutputPerc, names=FALSE, na.rm=TRUE) lri3 <- quantile(LU[LU[,"UsePP"]==TRUE & LU[,"UseGIP"]==FALSE, "LotLmRetailUsed"], LmOutputPerc, names=FALSE, na.rm=TRUE) lri4 <- quantile(LU[LU[,"UsePP"]==FALSE & LU[,"UseGIP"]==TRUE, "LotLmRetailUsed"], LmOutputPerc, names=FALSE, na.rm=TRUE) lri5 <- quantile(LU[LU[,"UsePP"]==FALSE & LU[,"UseGIP"]==FALSE, "LotLmRetailUsed"], LmOutputPerc, names=FALSE, na.rm=TRUE) lri6 <- quantile(LU[LU[,"UseGIP"]==FALSE,"LotLmRetailUsed"], LmOutputPerc, names=FALSE, na.rm=TRUE) lri7 <- quantile(LU[LU[,"UseGIP"]==TRUE, "LotLmRetailUsed"], LmOutputPerc, names=FALSE, na.rm=TRUE) qLeavingRetailByIntervenction <- cbind(rep("leaving retail", length(LmOutputPerc)), LmOutputPerc, lri1, lri2, lri3, lri4, lri5, lri6, lri7) # output the quantiles of each intervention at consumption c1 <- quantile(LU[,"conc.consumption.lm"], LmOutputPerc, names=FALSE, na.rm=TRUE) c2 <- quantile(LU[LU[,"UsePP"]==TRUE & LU[,"UseGIP"]==TRUE, "conc.consumption.lm"], LmOutputPerc, names=FALSE, na.rm=TRUE) c3 <- quantile(LU[LU[,"UsePP"]==TRUE & LU[,"UseGIP"]==FALSE, "conc.consumption.lm"], LmOutputPerc, names=FALSE, na.rm=TRUE) c4 <- quantile(LU[LU[,"UsePP"]==FALSE & LU[,"UseGIP"]==TRUE, "conc.consumption.lm"], LmOutputPerc, names=FALSE, na.rm=TRUE) c5 <- quantile(LU[LU[,"UsePP"]==FALSE & LU[,"UseGIP"]==FALSE, "conc.consumption.lm"], LmOutputPerc, names=FALSE, na.rm=TRUE) c6 <- quantile(LU[LU[,"UseGIP"]==FALSE,"conc.consumption.lm"], LmOutputPerc, names=FALSE, na.rm=TRUE) c7 <- quantile(LU[LU[,"UseGIP"]==TRUE, "conc.consumption.lm"], LmOutputPerc, names=FALSE, na.rm=TRUE) qConsumptionByIntervenction <- cbind(rep("at consumption", length(LmOutputPerc)), LmOutputPerc, c1, c2, c3, c4, c5, c6, c7) dfQuantilesbyIntervention <- rbind(qLeavingPlantByIntervenction, qLeavingRetailByIntervenction, qConsumptionByIntervenction) colnames(dfQuantilesbyIntervention) <- c("location", "percentile", "combined", "PP&GI", "PPonly", "GIonly", "none", "withoutGI", "withGI") write.csv(dfQuantilesbyIntervention, file="summary quantiles by intervention.csv", row.names = FALSE) # output the quantiles of each plant alternative at plant, retail, and consumption # same results as above, just uses Alternative labels rather than specific interventions qLeavingPlantByAlternative <- cbind(rep("leaving plant", length(LmOutputPerc)), LmOutputPerc, lpi1, lpi2, lpi3, lpi4, lpi5) qLeavingRetailByAlternative <- cbind(rep("leaving retail", length(LmOutputPerc)), LmOutputPerc, lri1, lri2, lri3, lri4, lri5) qConsumptionByAlternative <- cbind(rep("at consumption", length(LmOutputPerc)), LmOutputPerc, c1, c2, c3, c4, c5) dfQuantilesbyAlternative <- rbind(qLeavingPlantByAlternative, qLeavingRetailByAlternative, qConsumptionByAlternative) colnames(dfQuantilesbyAlternative) <- c("location", "percentile", "Combined", "AlternativeI", "AlternativeIIa", "AlternativeIIb", "AlternativeIII") write.csv(dfQuantilesbyAlternative, file="summary quantiles by alternative.csv", row.names = FALSE) # calculate and output counts of specific parameters # counts the number of Lots/FCS undergoing different actions # Lots used t1 <- table(Lots[,"PlantCategory"], factor(Lots[,"LotUse"], levels=c(FALSE,TRUE))) t1 <- fixtable(cat.ID, t1) # PP&GI tmp <- ifelse(Lots[,"UsePP"]==TRUE & Lots[,"UseGIP"]==TRUE & Lots[,"LotUse"]==TRUE,TRUE,FALSE) t2 <- table(Lots[,"PlantCategory"], factor(tmp, levels=c(FALSE,TRUE))) t2 <- fixtable(cat.ID, t2) # PPonly tmp <- ifelse(Lots[,"UsePP"]==TRUE & Lots[,"UseGIP"]==FALSE & Lots[,"LotUse"]==TRUE,TRUE,FALSE) t3 <- table(Lots[,"PlantCategory"], factor(tmp, levels=c(FALSE,TRUE))) t3 <- fixtable(cat.ID, t3) # GIonly tmp <- ifelse(Lots[,"UsePP"]==FALSE & Lots[,"UseGIP"]==TRUE & Lots[,"LotUse"]==TRUE,TRUE,FALSE) t4 <- table(Lots[,"PlantCategory"], factor(tmp, levels=c(FALSE,TRUE))) t4 <- fixtable(cat.ID, t4) # None tmp <- ifelse(Lots[,"UsePP"]==FALSE & Lots[,"UseGIP"]==FALSE & Lots[,"LotUse"]==TRUE,TRUE,FALSE) t5 <- table(Lots[,"PlantCategory"], factor(tmp, levels=c(FALSE,TRUE))) t5 <- fixtable(cat.ID, t5) # FCSTests t6 <- table(Lots[,"PlantCategory"], factor(Lots[,"FCSTest"], levels=c(FALSE,TRUE))) t6 <- fixtable(cat.ID, t6) # FCSPos t7 <- table(Lots[,"PlantCategory"], factor(Lots[,"FCSResult"] & Lots[,"FCSTest"], levels=c(FALSE,TRUE))) t7 <- fixtable(cat.ID, t7) # FCSTheorPos t8 <- table(Lots[,"PlantCategory"], factor(Lots[,"FCSResult"], levels=c(FALSE,TRUE))) t8 <- fixtable(cat.ID, t8) # LotTests t9 <- table(Lots[,"PlantCategory"], factor(Lots[,"LotTestbyFCS"] | Lots[,"LotTestSys"] | Lots[,"LotTestbyLot"], levels=c(FALSE,TRUE))) t9 <- fixtable(cat.ID, t9) # LotTestSys t10 <- table(Lots[,"PlantCategory"], factor(Lots[,"LotTestSys"], levels=c(FALSE,TRUE))) t10 <- fixtable(cat.ID, t10) # LotTestFCS t11 <- table(Lots[,"PlantCategory"], factor(Lots[,"LotTestbyFCS"], levels=c(FALSE,TRUE))) t11 <- fixtable(cat.ID, t11) # LotTestLot t11a <- table(Lots[,"PlantCategory"], factor(Lots[,"LotTestbyLot"], levels=c(FALSE,TRUE))) t11a <- fixtable(cat.ID, t11a) # LotTestPosPP t12 <- table(Lots[,"PlantCategory"], factor(Lots[,"LotLmPPResult"] & (Lots[,"LotTestbyFCS"] | Lots[,"LotTestSys"]), levels=c(FALSE,TRUE))) t12 <- fixtable(cat.ID, t12) # LotTheorPosPP t13 <- table(Lots[,"PlantCategory"], factor(Lots[,"LotLmPPResult"], levels=c(FALSE,TRUE))) t13 <- fixtable(cat.ID, t13) # LotTheorRetailUsedPos t14 <- table(Lots[,"PlantCategory"], factor(Lots[,"LotLmRetailUsedResult"], levels=c(FALSE,TRUE))) t14 <- fixtable(cat.ID, t14) dfCounts <- data.frame("Plant type" = cat.ID, "LotsDesired" = nLotSim, "LotsUsed" = t1[,2], "PP&GI" = t2[,2], "PPonly" = t3[,2], "GIonly" = t4[,2], "None" = t5[,2], "FCSTests" = t6[,2], "FCSPos" = t7[,2], "FCSTheorPos" = t8[,2], "LotTests" = t9[,2], "LotTestSys" = t10[,2], "LotTestFCS" = t11[,2], "LotTestLot" = t11a[,2], "LotTestPosPP" = t12[,2], "LotTheorPosPP" = t13[,2], "LotTheorRetailUsedPos" = t14[,2], "FracFCSPos" = t7[,2]/t6[,2], "FracLotPos" = t12[,2]/t9[,2]) write.csv(dfCounts, file="results count summary.csv", row.names=FALSE) # calculate and output processing results prevalence.FCS <- sum(Lots[,"FCSResult"], na.rm=TRUE) / nlot.total prevalence.PP <- sum(Lots[,"LotLmPPResult"], na.rm=TRUE) / nlot.total prevalence.retail <- sum(Lots[,"LotLmRetailUsedResult"], na.rm=TRUE) / sum(Lots[,"LotUse"]) lots.fracPPGIP <- sum(t2[,2])/sum(t1[,2]) lots.fracPP <- sum(t3[,2])/sum(t1[,2]) lots.fracGIP <- sum(t4[,2])/sum(t1[,2]) lots.fracNone <- sum(t5[,2])/sum(t1[,2]) TotalProduction <- sum(Lots[,"LotMass"]) fracProductionWithGI <- sum(Lots[Lots[,"UseGIP"]==TRUE,"LotMass"])/TotalProduction fracProductionWithoutGI <- sum(Lots[Lots[,"UseGIP"]==FALSE,"LotMass"])/TotalProduction fracProductionWithPP <- sum(Lots[Lots[,"UsePP"]==TRUE,"LotMass"])/TotalProduction fracProductionWithoutPP <- sum(Lots[Lots[,"UsePP"]==FALSE,"LotMass"])/TotalProduction tmp <- Lots[Lots[,"UsePP"]==TRUE & Lots[,"UseGIP"]==TRUE,] fracProductionBothPPGI <- sum(tmp[,"LotMass"])/TotalProduction FoodProductionLU <- sum(LU[,"LotMass"]) fracFoodWithGILU <- sum(LU[LU[,"UseGIP"]==TRUE,"LotMass"])/FoodProductionLU fracFoodWithoutGILU <- sum(LU[LU[,"UseGIP"]==FALSE,"LotMass"])/FoodProductionLU fracFoodWithPPLU <- sum(LU[LU[,"UsePP"]==TRUE,"LotMass"])/FoodProductionLU fracFoodWithoutPPLU <- sum(LU[LU[,"UsePP"]==FALSE,"LotMass"])/FoodProductionLU tmpLU <- LU[LU[,"UsePP"]==TRUE & LU[,"UseGIP"]==TRUE,] fracFoodBothPPGILU <- sum(tmpLU[,"LotMass"])/FoodProductionLU dfProcessing <- c("FCS Prevalence" = prevalence.FCS, "PostProcessing Prevalence" = prevalence.PP, "Actual Retail Prevalence" = prevalence.retail, "Processing breakdown" = "", "Fraction of lots with both PP and GI" = lots.fracPPGIP, "Fraction of lots with only PP" = lots.fracPP, "Fraction of lots with only GI" = lots.fracGIP, "Fraction of lots with neither PP nor GI" = lots.fracNone, "Check - sum of fractions" = sum(lots.fracPPGIP, lots.fracPP, lots.fracGIP, lots.fracNone), "Fraction of lots without GI" = lots.fracPP+lots.fracNone, "Fraction of lots with GI" = lots.fracGIP+lots.fracPPGIP, "TOTAL production (ignoring sample testing)" = TotalProduction, " Fraction of production with GI" = fracProductionWithGI, " Fraction of production without GI" = fracProductionWithoutGI, " Fraction of production with PP" = fracProductionWithPP, " Fraction of production without PP" = fracProductionWithoutPP, " Fraction of production with both GIP and PP" = fracProductionBothPPGI, "Total actually used production" = FoodProductionLU, " Fraction of used production with GI" = fracFoodWithGILU, " Fraction of used production without GI" = fracFoodWithoutGILU, " Fraction of used production with PP" = fracFoodWithPPLU, " Fraction of used production without PP" = fracFoodWithoutPPLU, " Fraction of used production with both GIP and PP" = fracFoodBothPPGILU, "Total lots simulated" = nlot.total, "Total lots used after sampling" = nrow(LU), "Count of lots used above threshold" = "", " >= 0.04 cfu/g leaving plant" = length(LU[LU[,"LotLmPP"] >= 0.04, "LotLmPP"]), " >=0.04 cfu/g slicing" = length(LU[LU[,"LotLmRetailSlicing"] >= 0.04, "LotLmRetailSlicing"]), " >=0.04 cfu/g at consumption" = length(LU[LU[,"conc.consumption.lm"] >= 0.04, "consumption.lm"]), " >=100 cfu/g at consumption" = length(LU[LU[,"conc.consumption.lm"] >= 100, "consumption.lm"])) write.csv(dfProcessing, file="Processing results.csv", row.names = TRUE) # output plots #plot quantiles PlotQuantileComparison(run.no, LmOutputPerc, c1, c2, c3, c4, c5, PlotsToDisk) PlotQuantileComparisonBaseline(run.no, LmOutputPerc, lpi1, lri1, c1, PlotsToDisk) # plot prevalences prevalence <- c(prevalence.FCS, prevalence.PP, prevalence.retail) prevalence.names <- c("FCS", "Post Processing", "Retail") PlotPrevalences(prevalence, prevalence.names, PlotsToDisk) # plot fraction of alternatives fracAlternatives <- c(lots.fracPPGIP, lots.fracPP, lots.fracGIP, lots.fracNone) fracTitles <- c("PP&GIP", "PP only", "GI only", "None") PlotFractions(fracAlternatives, fracTitles, PlotsToDisk) # plot testing results graph.labels <- cat.ID[frac!=0] PlotTests(cat.ID, nLotSim, t6[,2], t7[,2]/t6[,2], t9[,2], t12[,2]/t9[,2], PlotsToDisk) # plot sample time series if desired PlotSampleTimeSeries(Lots, PlotsToDisk)