#' @rdname Klassifikation
#' @title Klassifikation
#' @name Klassifikation
#' @description Classification Table classification_table
#' Richtige und falsche Klassifikationen
#' Bei 2x2 Tabellen der Kappa Test
#'
#' Sensitivity = A/(A+C)
#'
#' Specificity = D/(B+D)
#'
#' Prevalence = (A+C)/(A+B+C+D)
#'
#' PPV = (sensitivity * prevalence)/((sensitivity*prevalence) + ((1-specificity)*(1-prevalence)))
#'
#' NPV = (specificity * (1-prevalence))/(((1-sensitivity)*prevalence) + ((specificity)*(1-prevalence)))
#'
#' Detection Rate = A/(A+B+C+D)
#'
#' Detection Prevalence = (A+B)/(A+B+C+D)
#'
#' Balanced Accuracy = (sensitivity+specificity)/2
#'
#' Precision = A/(A+B)
#'
#' Recall = A/(A+C)
#'
#' F1 = (1+beta^2)*precision*recall/((beta^2 * precision)+recall)
#'
#'
#' Klassifikation fuer Binominal-GLM
#' und zeigt die ...
#' @param x glm oder xtab Objekt
#' @param thresh Klassifikation auf Basis der Vorhersage Schwelle bei P=0.5
#' @param caption an Output
#' @param ... weitere Objekte nicht benutzt
#' @return A data.frame Objekt.
#' @export
#' @examples
#'
#' ##https://www.hranalytics101.com/how-to-assess-model-accuracy-the-basics/
#'
#'
#'
#' Kappa2(x<-xtabs(~gruppe+lai, hkarz))
#' head(hkarz)
#'
#' fit1<- glm(gruppe~lai, hkarz, family = binomial)
#' #thkarz <- as.data.frame(xtabs(~gruppe+lai, hkarz))
#' #fit2<- glm(Freq ~ gruppe*lai, thkarz, family = poisson())
#'
#' tab <- Klassifikation(x)$statistic[c(1,7,8),]
#' x2 <- Klassifikation(fit1)
#' tab$fit <- x2$statistic[c(1,7,8), 2]
#' tab
#'
#' require(pROC)
#' roc_curve <- roc(x2$response, x2$predictor)
#' windows(8,8)
#' plot(roc_curve, print.thres = "best",
#' print.auc=TRUE)
#' abline(v = 1, lty = 2)
#' abline(h = 1, lty = 2)
#' #text(.90, .97, labels = "Ideal Model")
#' #points(1,1, pch = "O", cex = 0.5)
Klassifikation <- function(x, ...) {
UseMethod("Klassifikation")
}
#' @rdname Klassifikation
#' @export
Klassifikation2<- function(x, ..., caption="",
include.xtab=TRUE,
include.statistics=TRUE) {
res <- Klassifikation(x, ...)
if(include.xtab)
Output(res$xtab,
caption=caption)
if(include.statistics & !is.null(res$statistic))
Output(res$statistic,
caption=caption)
invisible(res)
}
#' @rdname Klassifikation
#' @description Klassifikation.glm
#' @export
#' @examples
#' #fit1<- glm(gruppe~lai, hkarz, family = binomial)
#' #thkarz <- as.data.frame(xtabs(~gruppe+lai, hkarz))
#' #fit2<- glm(Freq ~ gruppe*lai, thkarz, family = poisson())
#'
Klassifikation.glm<- function(x,
thresh = 0.5,
caption="Klassifikation",
...){
response <- all.vars(formula(formula(x)))[1]
data <- x$model
predictor <- fitted(x) # vorhergesagte Wahrscheinlichkeit
data$Response <- data[, response]
# cat("\nKontrolle:\n")
# print(head(data))
#response predictor
mylevels <- if(is.factor(data$Response))
levels(data$Response)
else 0:1
data$Predictor <- cut(predictor,
breaks=c(-Inf, thresh, Inf), labels=mylevels )
# print(head(data))
# Kontingenztafel: tatsaechliche vs. vorhergesagte Kategorie
cTab <- stats::xtabs( ~ Response + Predictor, data = data)
if (length(cTab)==4){
res <- Klassifikation.xtabs(cTab)
res$response=data$Response
res$predictor=predictor
# res
}
else res <- list(xtab=cTab,
statistic=NULL,
response=response,
predictor=predictor)
return(res)
}
#' @rdname Klassifikation
#' @export
Klassifikation.table <- function(...) {
Klassifikation.xtabs(...)
}
#' @rdname Klassifikation
#' @description xtabs-Objekt
#' @export
#' @examples
#' hkarz$LAI<- factor(hkarz$lai, 0:1, c("pos", "neg"))
#' APA2(xtabs(~gruppe+LAI, hkarz), test=TRUE, type="fischer")
Klassifikation.xtabs<- function(x,
lvs = c("positiv", "negativ"),
digits = 2,
prevalence = NULL,
# caption="Klassifikation",
# note="",
...){
if(!length(x)==4) return("Nur mit 2x2 Tabellen moeglich!")
Positive_Class <-
paste(attr(x, "dimnames")[[1]][1],
attr(x, "dimnames")[[2]][1], sep ="/")
attr(x, "dimnames")[[1]] <- lvs
attr(x, "dimnames")[[2]] <- lvs
xtab<-x
x <- caret::confusionMatrix(x, prevalence = prevalence)
out <- as.character(c(
rndr_(x$overall["Accuracy"], digits),
rndr_CI(x$overall[c("AccuracyLower", "AccuracyUpper")]),
rndr_(x$overall["AccuracyNull"], digits),
rndr_P(x$overall["AccuracyPValue"]),
rndr_(x$overall["Kappa"], digits),
rndr_P(x$overall["McnemarPValue"]),
rndr_(x$byClass, digits)
,Positive_Class))
list(xtab=xtab,
statistic=
data.frame(Statistic =
c("Accuracy",
"95% CI",
"No Information Rate",
"P-Value [Acc > NIR]",
"Kappa",
"Mcnemar's Test P-Value",
"Sensitivity",
"Specificity",
"Pos Pred Value" ,
"Neg Pred Value",
"Precision",
"Recall",
"F1",
"Prevalence",
"Detection Rate",
"Detection Prevalence" ,
"Balanced Accuracy",
"Positive Class"
)
, Value = out,
stringsAsFactors = FALSE),
stat=x,
response=NULL,
predictor=NULL
)
}
#' APA Methode für epi.tests
#'
#' @param x epi.tests Objekt
#' @param ... an Output
#'
#' @rdname APA2
#'
#' @return data.frame
#' @export
#'
#' @examples
#'
#' #' library(caret)
#' require(epiR)
#' DF<- GetData("
#' GoldStandart RT.qPCR Anzahl
#' positiv positiv 111
#' positiv negativ 12
#' negativ positiv 1
#' negativ negativ 62 ", Tabel_Expand =TRUE, id.vars=1:2, output=FALSE)
#' DF$GoldStandart<- factor( DF$GoldStandart, rev(levels( DF$GoldStandart)))
#' DF$RT.qPCR<- factor( DF$RT.qPCR, rev(levels( DF$RT.qPCR)))
#'
#'
#'
#'
#' dat <- as.table(matrix(c(111, 12, 1, 62), nrow = 2, byrow = TRUE))
#' colnames(dat) <- c("RT.qPCR +", "RT.qPCR -")
#' rownames(dat) <- c("Gold Standart +", "GoldStandart -")
#' rval <- epi.tests(dat, conf.level = 0.95)
#' APA2(rval)
#'
#'
APA2.epi.tests <-
function (x,
caption=paste("Point estimates and", x$conf.level * 100, "%","CIs:"),
...)
{
# Output(cbind(Test= row.names( x$tab), x$tab), ...)
stat <- with(x$rval, {
data.frame(
"Diagnostic Parameters"
= c(
"Apparent prevalence",
"True prevalence",
"Sensitivity",
"Specificity",
"Positive predictive value",
"Negative predictive value",
# "Positive likelihood ratio",
# "Negative likelihood ratio",
"Diagnostic Accuracy"
) ,
"Point Estimates" = c(
sprintf("%.2f (%.2f, %.2f)", aprev$est, aprev$lower, aprev$upper),
sprintf("%.2f (%.2f, %.2f)", tprev$est, tprev$lower, tprev$upper),
sprintf("%.2f (%.2f, %.2f)", se$est, se$lower, se$upper),
sprintf("%.2f (%.2f, %.2f)", sp$est, sp$lower, sp$upper),
sprintf("%.2f (%.2f, %.2f)", ppv$est, ppv$lower, ppv$upper),
sprintf("%.2f (%.2f, %.2f)", npv$est, npv$lower, npv$upper),
# sprintf("%.2f (%.2f, %.2f)", plr$est, plr$lower, plr$upper),
# sprintf("%.2f (%.2f, %.2f)", nlr$est, nlr$lower, nlr$upper),
sprintf("%.2f (%.2f, %.2f)", diag.acc$est, diag.acc$lower, nlr$upper))
)
})
Output(stat,
caption=caption,
...)
invisible(stat)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.