# FUNCTION to compute the a-posteriori error considering a population stratified into risk geroups. For a sampling scheme designed to substantiate freedom from disease the function computes the a-posteriori alpha-error, i.e., the actual alpha-error based on the drawn sample, when the population is stratified into risk groups for which the relative infection risk is known. ```r computeAposterioriErrorRiskGroups(alphaErrorVector, groupVec, groupLevels, nPopulationVec, nRelRiskVec, nDiseased, method = "default") ``` ## Arguments - `alphaErrorVector`: Numeric vector. Alpha-error (between 0 and 1) of each herd in the sample. - `groupVec`: Character vector. Risk group to which each of the herds in the sample belongs. Must have the same length (and order) as `alphaErrorVector`. - `groupLevels`: Character vector. (Unique) levels/names of the risk group. Defines the order of the values in `nPopulationVec` and `nRelRiskVec`. - `nPopulationVec`: Integer vector. Population sizes of the risk groups. Must have the same length (and order) as `groupLevels`. - `nRelRiskVec`: Numeric vector. (Relative) infection risks of the risk groups. Must have the same length (and order) as `groupLevels`. - `nDiseased`: Integer. Number of diseased herds in the population according to the design prevalence. - `method`: Character string. "exact" for exact error, "approx" for approximation (recommended for nDiseased > 7). ## Details The exact evaluation of the alpha-error is computationally complex, due to combinatirical issues. In order to increase effectivity parts of the code were implemented in C. Still, for `nDiseased` > 3 the computation may take very long and it is generally not recommended to use the exact calculation. Rather the approximation should be used for `nDiseased` > 3. ## Returns The return value is the a-posteriori alpha-error based on the sample at hand (numeric scalar). ## Author(s) Ian Kopacka ## Examples ```r data(sheepData) sheepData$size <- ifelse(sheepData$nSheep < 30, "small", "large") riskValueData <- data.frame(riskGroup = c("small", "large"), riskValues = c(1,2)) mySurvey <- surveyData(nAnimalVec = sheepData$nSheep, riskGroupVec = sheepData$size, riskValueData = riskValueData, populationData = sheepData, designPrevalence = 0.002, alpha = 0.05, intraHerdPrevalence = 0.13, diagSensitivity = 0.9, costHerd = 30, costAnimal = 7.1) ## Limited sampling with risk groups: myLtdSamplingRG <- ltdSampling(survey.Data = mySurvey, sampleSizeLtd = 20, nSampleFixVec = NULL, probVec = c(1,2)) ## Draw sample manually: set.seed(20110708) x <- myLtdSamplingRG indexSampleRG <- sapply(seq(along = x@surveyData@riskValueData$riskGroup), function(ii){ riskGroup <- as.character(x@surveyData@riskValueData$riskGroup[ii]) nSampleRG <- x@nHerdsPerRiskGroup[riskGroup] indexRG <- which(x@surveyData@riskGroupVec == riskGroup) indexOut <- sample(x = indexRG, size = nSampleRG, replace = FALSE) }) indexSample <- sort(Reduce(function(x,y) c(x,y), indexSampleRG)) ## Compute the a-posteriori alpha error: alphaErrorVector <- computeAlpha(nAnimalVec = x@surveyData@nAnimalVec[indexSample], method = "limited", sampleSizeLtd = x@sampleSizeLtd, intraHerdPrevalence = x@surveyData@intraHerdPrevalence, diagSensitivity = x@surveyData@diagSensitivity, diagSpecificity = 1) ## Determine the number of herds in each risk group: riskValueDf <- mySurvey@riskValueData[,1:2] names(riskValueDf) <- c("riskGroup", "riskValues") riskValueDf$riskGroup <- as.character(riskValueDf$riskGroup) riskValueDf$id <- seq(along = riskValueDf[,1]) riskGroupTab <- table(mySurvey@riskGroupVec) riskGroupDf <- data.frame(riskGroup = as.character(names(riskGroupTab)), nPopulation = as.vector(riskGroupTab)) riskValueDf <- merge(x = riskValueDf, y = riskGroupDf, by = "riskGroup", sort = FALSE) riskValueDf <- riskValueDf[order(riskValueDf$id, decreasing = FALSE),] aPostAlphaRG <- computeAposterioriErrorRiskGroups(alphaErrorVector = alphaErrorVector, groupVec = x@surveyData@riskGroupVec[indexSample], groupLevels = riskValueDf$riskGroup, nPopulationVec = riskValueDf$nPopulation, nRelRiskVec = riskValueDf$riskValues, nDiseased = max(round(length(x@surveyData@nAnimalVec)*x@surveyData@designPrevalence),1), method = "approx") ```

computeAposterioriErrorRiskGroups function