####################### #This document is a material for the FDA/FSIS Assessment of the risk of salmonellosis linked #to the consumption of liquid egg products made from internally contaminated shell eggs #initially stored at 18°C (65°F) compared to 7°C (45°F). # #This information does not represent and should no be construed to represent any FDA or FSIS #determination or policy or recommended code. ####################### #Note: please set the Working Directory to the Source File Location library(readxl) file <- "ShellEggModel2019.xlsm" #Version of sample without the behavior if x is a single value. sample <- function (x, size, replace = FALSE, prob = NULL) { x[sample.int(length(x), size, replace, prob)] } #Define Output File #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv fileOutput <-"Simul.csv" #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #Parameters #Load the set of 100,000 simulations: each value is the contamination level of one egg DataEggs <- read_xlsx(path = file, sheet = "Results", range = "R19C9:R100018C9", col_names = FALSE)[[1]] DataEggs <- DataEggs[!is.na(DataEggs)] dim(DataEggs) #Load the type of contamination for each eggs #ote: the same for all simulations as the random seed is controlled in the Egg SHell VBA code DataTypeEggs <- read_xlsx(path = file, sheet = "Results", range = "R19C2:R100018C2", col_names = FALSE)[[1]] DataTypeEggs <- DataTypeEggs[!is.na(DataTypeEggs)] #Quantiles of the output distribution (correspond to the one used in the Liquid Egg model) quant <- c(seq(0,0.9,length=2501), seq(0.9,0.95,length=2500)[-1], seq(0.95,0.99,length=2500)[-1], seq(0.99,1,length=2501)) #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #Major function #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv Simul <- function(Data, Prevalence, VatMin = 50000, VatMax = 350000, test=FALSE, vtest=100, salt = FALSE, nDist = 1E6){ # Data: a list # $EggType is the type of contamination (vector): #"Alb C G", "Shell", "Alb F N","Alb F G", "Alb C N","VM Low","VM High","Yolk Low","Yolk High" # $nBact is the matrix of contamination. Three columns: "Whole", "White", "Yolk" # Prevalence: prevalence of contaminated eggs # VatMin = 50000, VatMax = 350000, : minimal and maximal vat size (in lbs). Will be sed in a Uniform distribution # test: is there any bacteriological test applied # vtest: size of the sampling test (in g) # salt = FALSE: is it salted or non salted product (will change the pasteurization value for the test) # nDist: size of the sample used to derive the empirical distribution. ##### Function ### # This function is launched for each "TypeEgg" products # It provides a vector of nDist random values of the nb of bacteria in the vat. dist <- function(Data, type){ #Data: same as Data for the function #type: "White", "Yolk" OR "Whole" #find the index of eggs that are of the "type" type quel <- Data$EggType %in% Type[[type]] #extract the contamination data for these eggs BactPreProcess <- Data$nBact[quel] #number of values nValues <- sum(quel) #The prevalence of contaminated "White", "Yolk" or "Whole" is the product of the prevalence of contaminated #eggs by the probability that the egg was contaminated on the "White", the "Yolk" or the "Whole" Prev <- Prevalence * nValues / length(quel) #Set seed (always the same) set.seed(1234) #How many parts are contaminated: binomial draw from NProdVatWeighted by Prev howMany <- rbinom(nDist, round(NProdVatWeighted), Prev) #nb bact Tot: for each simulations, sample among the BactPreProcess the number of contaminated eggs # and sum nbBactTot <- sapply(1:nDist, function(x){ sum(sample(BactPreProcess,howMany[x],replace=TRUE)) } ) return(nbBactTot) } ##### Some Variables TypeEgg <- c("White","Yolk","Whole") #Weight of egg parts Weight <- list(White=2/3*50,Yolk=1/3*50,Whole=50) #which type of egg contamination is used for which parts. Type <- list(White = c("Alb C G", "Shell", "Alb F N","Alb F G", "Alb C N"), Yolk = c("VM Low","VM High","Yolk Low","Yolk High") ) #Whole is all Type[["Whole"]] <- c(Type[["White"]],Type[["Yolk"]]) ##### set.seed(5678) #Size of the vat gVat <- runif(nDist, min=VatMin, max=VatMax) * 453.592 # (g) #the larger the vat, the more it should be represented for one serving. #we resample the size proportionnaly to ... the size gVatWeighted <- sample(gVat, nDist, replace=TRUE, prob = gVat) Res <- BactConc <- NULL for(type in TypeEgg){ cat("Type:",type,"\n") #Number of product parts in the vat NProdVatWeighted <- gVatWeighted / Weight[[type]] #The bact concentration is equal to the number of bacteria in the vat divided by the vat size (g) nbBactConc <- dist(Data, type) / gVatWeighted if(test) { # If test, resample proportionnaly to the prob of being negative in a test # and proportionnaly to the vat size (so that large vat are proportionnaly more frequently sampled) # fud the pasteurization log reduction past <- switch(type, "White"= -5, "Yolk" = ifelse(salt,-7.2,-5.5), #-7.2,# -5.5, "Whole" = ifelse(salt,-6.0, -5.9))# -6.0) # nbBactConc <- sample(nbBactConc, nDist, replace =TRUE, prob = gVatWeighted * exp(-vtest * nbBactConc * 10^past)) } Conc <- as.vector(quantile(nbBactConc,probs=quant)) Res <- cbind(Res, Conc) BactConc <- cbind(BactConc, nbBactConc) } # end loop type colnames(Res) <- TypeEgg return(list(Res = Res, BactConc =BactConc)) } Resi <- quant # Data <- list(EggType = DataTypeEggs,nBact=DataEggs[,1]) # Res <- Simul(Data, Prevalence=0.03/100) # Resi <- cbind(Resi, Res$Res) # Res <- Simul(Data, Prevalence=0.03/8.706115215/100) # 130,017 / 14,934 = 8.706115215 # Resi <- cbind(Resi, Res$Res) # write.csv(Resi, fileOutput) #Run the simulation Data <- list(EggType = DataTypeEggs,nBact=DataEggs) Res <- Simul(Data, Prevalence= 0.03/100/8.706115215) Resi <- cbind(Resi, Res$Res) #Write the output write.csv(Resi, fileOutput)