Kapitel 38 Physiotherapeutische Diagnosen stellen

Ein Beispiel für eine physiotherapeutische Diagnose ist das Bestimmen eines Sturzrisikos.

Physiotherapeut:innen stellen keine medizinischen Diagnosen, sondern physiotherapeutische. Wir stellen zum Beispiel fest, ob jemand ein Gleichgewichtsproblem hat, ob die Patient:innen eine erhöhte Sturzgefahr aufweisen, oder ob ein Schmerz eher vom unteren Rücken kommt oder von der Hüfte. Hier werden jetzt viele Mediziner sagen, dass das letztere ihr Job ist. Doch das stimmt nicht, da viele ärztlichen Verordnungen nur grobe “Diagnosen” aufweisen. Ausserdem müssen wir Physiotherapeut:innen fähig sein, bei einem Patienten mit Rückenschmerzen im Verlaufe der Therapie zu erkennen, ob sich eine neurologische Beteiligung entwickelt. Auch sollten wir Red Flags beachten, die manchmal erst im Verlaufe der Behandlung - also nach der ärztlichen Untersuchung - auftreten.

Wenn wir eine Diagnose stellen, tendieren wir dazu dichotom zu denken: “Die Patientin ist Sturzgefährdert” oder “Die Patientin ist nicht Sturzgefährdet”. Das Argument ist: “Wir müssen ja entscheiden, ob wir etwas tun oder nicht tun, das ist halt in [Hier einsetzen was ihr möchtet] Namen dichotom”. Ja, wir müssen entscheiden etwas zu tun oder nicht, aber wir müssen lernen, diese Entscheidung auf Wahrscheinlichkeiten zu basieren - mit dem Wissen, dass wir eben auch falsch liegen können. Wir müssen die Unsicherheit akzeptieren und aus ethischen Gründen auch mit den Patient:innen diskutieren.

Idealerweise würde uns ein diagnostischer Test eine Wahrscheinlichkeit geben, dass ein:e Patient:in das gesuchte Problem hat. Wie könnten wir das tun: Wenn wir aus grossen Studien für ein Set von diagnostisch Wertvollen Variablen die Koeffizienten aus einer logistischen Regression haben, so können wir direkt die Wahrscheinlichkeit berechnen. Für bestimmte medizinischen Diagnosen gibt es dies schon (hier klicken für ein Beispiel), in der Physiotherapie gibt es nur wenige Beispiele. Siehe hier ein Beispiel (hier klicken) aus Australien zur Sturzrisikoberechnung oder hier ein europäisches Tool.

Wir benutzen jedoch oft einen Umweg. Wir haben Tests, die an einem Schwellenwert dichotomisiert werden und uns ein positives oder ein negatives Testresultat geben. Mit der aus Studien bekannten Sensitivität und Spezifität können wir die Likelihood Ratio für einen positiven Test oder für einen negativen Test berechnen. Anhand der Likelihood Ratio und der Vortestwahrscheinlichkeit berechnen wir nun - oder lesen es aus einem Fagan-Nomogramm ab - die Nachtestwahrscheinlichkeit berechnen.

Sie finden ein online Fagan-Nomogramm hier.

Beim Gebrauch von Sensitivität und Spezifität machen wir oft den Fehler, aus der Wahrscheinlichkeit eines positiven Tests bei einer Krankheit A direkt auf die Wahrscheinlichkeit dieser Krankheit bei einem positiven Test zu schliessen. Wir werden in diesem Kapitel auch lernen, warum dies sehr problematisch ist.

38.1 Was lernen wir in diesem Kapitel

Wir lernen in diesem Kapitel, dass wir in der Physiotherapie oft Sensitivität und Spezifität benutzen, respektive die Likelihood-Ratios für positive oder negative Tests. Wir lernen aber auch, dass dieser Ansatz eigentlich problematisch ist, da er auf einer Dichotomisierung beruht. Um die Sensitivität und die Spezifität zu berechnen, müssen wir ein Testresultat in positiv und negativ einteilen. Das bestimmen der dazu nötigen Schwellenwerte oder Cut Offs ist jedoch fast etwas eine Glücksache. Besser wäre es, mit Wahrscheinlichkeiten zu arbeiten. Diese Wahrscheinlichkeiten kommen, genau wie die Sensitivitäten und Spezifitäten aus grössern Studien und werden zum Beispiel mit logistischen Regressionen berechnet. Würden die Studien die Koeffizienten aus den logistischen Regressionen veröffentlichen, könnten wir diese für das Berechnen der Wahrscheinlichket einer Krankheit benutzen. Leider werden diese Koeffizienten zu wenig oft angegeben oder stammen aus zu kleinen Studien.

Begriffe, die Sie nach diesem Kapitel interpretieren können:

• Wahrscheinlichkeit einer Krankheit / Problematik
• Logistische Formel um die Wahrscheinlichkeit zu berechnen
• Sensitivität und Spezifität
• Likelihood Ratio für einen positiven Test
• Likelihood Ratio für einen negativen Test
• Cut Offs und warum man sie eigentlich nicht benutzen sollte
• Unterschied Probabilität Test+ | Disease+ und Probabilität Disease+ | Test+ (das liest sich übrigens: Wahrscheinlichket, dass der Test positive ist, gegeben dass die Krankheit wirklich vorhanden ist. oder (2.Teil): Wahrscheinlichkeit, dass die Krankheit wirklich vorhanden ist, gegeben, dass der Test positiv ist).
• Einfluss der Vortestwahrscheinlichkeit auf die Nachtestwahrscheinlichkeit

38.2 Vorhersage Sturzrisiko

Im ersten Beispiel möchten wir ein erhöhtes Sturzrisiko diagnostizieren.

Idealerweise hätten wir ein Tool, in das wir die Testresultate - und zwar als Rohdaten, das heisst, nicht dichotomisiert, sonder den effektiven Wert des Tests eingeben. Das Tool berechnent dann mit den Koeffizienten aus einer logistischen Regression die Wahrscheinlichkeit, dass ein:e Patient:in in einer bestimmten Zeit stürzt. Ein Tool, dass dies schon fast erfüllt finden Sie hier, leider muss man dort die für viele Tests die dichotomisierten Werte eingeben, was nicht optimal ist, aber es ist trotzdem ein gutes Tool. (Das Tool müsste jetzt noch extern evaluiert werden. Eine solche externe Validation ist jedoch sehr teuer und wird deswegen zu wenig durchgeführt.) Den Artikel zur Entwicklung und internen Validität finden Sie hier

Obschon ja jede Person ein Sturzrisiko hat, wird heute immer noch meistens dichotomisiert in Sturzrisiko ja versus Sturzrisiko nein. Das macht eigentlich keinen Sinn. Das zeigen wir hier auf.

In der Physiotherapie wird leider immer noch häufig der Timed-up-and-go Test zur Bestimmung eines Sturzrisikos benutzt, deswegen brauchen wir den auch für das erste Beispiel. Warum leider? Der Test ist sehr gut für Verlaufsmessungen zu machen, er taugt jedoch nicht viel zur Bestimmung des Sturzrisikos Beauchet et al 2011 (hier klicken) oder diese Meta-Analyse.

df<-data_script[,c("id", "Faller", "TUG")]

Zuerst erstellen wir eine Graphik, die die Wahrscheinlichkeit zu stürzen über die benötigte Zeit für den TUG graphisch darstellt.

Dafür müssen wir mit einer logistischen Regression die Koeffizienten berechnen. Danach können wir für jede:n Patient:in die Wahrscheinlichkeit in den nächsten zwölf Monaten zu stürzen berechnen. Wir können das mit der logistischen Formel tun, oder mit der R Funktion predict.

Logistische Regression

lr<-glm(Faller~TUG, data=df, family = binomial)
summary(lr)

Call:
glm(formula = Faller ~ TUG, family = binomial, data = df)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.364  -1.127  -1.079   1.228   1.292  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.39974    0.27522  -1.452    0.146
TUG          0.02372    0.01994   1.190    0.234

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 517.61  on 373  degrees of freedom
Residual deviance: 516.18  on 372  degrees of freedom
AIC: 520.18

Number of Fisher Scoring iterations: 3

Wir berechnen zu Anschaungszwecken die Wahrscheinlichkeiten manuell.

Im ersten Schritt berechnen wir die Log Odds für jeden Patienten:

df<-df %>% 
  mutate(manuell_berechnete_LogOdds=-0.39974+TUG*0.02372)

Im zweiten Schritt wandeln wir die LogOdds in Odds um:

df<-df %>% 
  mutate(manuell_berechnete_Odds=exp(manuell_berechnete_LogOdds))

Im dritten Schritt wandeln wir die Odds in Wahrscheinlichkeiten um:

df<-df %>% 
  mutate(manuell_berechnete_Wahrscheinlichkeit=manuell_berechnete_Odds/(1+manuell_berechnete_Odds))

Jetzt benutzen wir die R Funktion predict, die Variable mit den Wahrscheinlichkeiten heisst hier dann fit.

predicted.data<- as.data.frame(predict(lr, type="response", se=TRUE))
df <- cbind(df, predicted.data)

Jetzt schauen wir, ob die manuell berechneten Wahrscheinlichkeiten mit denen von der R Funktion übereinstimmen. Diese Graphik ist nur zur Kontrolle 😄, da wir ja in der Realität nicht beide Varianten berechnen, sondern einfach nur die predict Funktion benutzen.

ggplot(df, aes(manuell_berechnete_Wahrscheinlichkeit, fit))+
  geom_point()

Abbildung 22.1: Wir sind zufrieden, unsere Berechnung ergibt genau die gleichen Wahrscheinlichkeiten, wie die R Funktion predict.

Wir plotten die TUG-Zeit gegen die Wahrscheinlichkeit zu stürzen. Dazu müssen wir zuerst die beobachteten (empirischen) Stürze in die beobachteten Wahrscheinlichkeiten umrechnen.

Zuerst erstellen wir Gruppen basiert auf dem TUG, wir erstellen alle 5 Sekunden im TUG einen neue Gruppe

df<-df %>% 
    mutate(
        tug_cut_width = cut_width(TUG, width=5),
        ) 

plot(df$TUG, df$tug_cut_width)

Abbildung 1.2: Zeit im TUG aufgeteilt in Gruppen.

Wir sehen, dass in den drei langsamsten Gruppen nur sehr wenige Personen sind, deswegen werden wir die letzten vier Gruppen kombinieren.

df<-df %>% 
    mutate(
        tug_cut_width = cut(TUG,c(0,2.5,7.5,12.5,17.5,22.5,100)),
        ) 

plot(df$TUG, df$tug_cut_width)

Abbildung 1.3: Zeit im TUG aufgeteilt in Gruppen.

Mit dem folgenden Code berechnen wir die Wahrscheinlichkeit, dass jemand innerhalb einer Gruppe gestürzt ist.

Dank an den Autor vom Cookbook-r

# New version of length which can handle NA's: if na.rm==T, don't count them
    length2 <- function (x, na.rm=FALSE) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
    }



df_summary<-df %>% 
  group_by(tug_cut_width) %>% 
  summarise(Freq = length2(tug_cut_width, na.rm=TRUE),
            empirical_prob_fall=sum(as.numeric(as.character(Faller)), na.rm=TRUE)/sum(Freq))

ggplot(df_summary, aes(x=tug_cut_width, y=empirical_prob_fall))+
  geom_point(size=4, color="grey")+
  scale_y_continuous(limits = c(0.0, 1), breaks = seq(0, 1, by = 0.1))+
  theme_classic()

Abbildung 21.6: Empirische Wahrscheinlichkeit in den nächsten 12 Monaten zu stürzen in Gruppen nach TUG Zeit. Wir sehen, dass die Wahrscheinlichkeit nach 17.5 Sekunden leicht ansteigt. Das ist jetzt aber kein Argument, dort einen Cut-Off zu setzen, da danach das Risiko konstant ansteigt. Würden wir einen Cut-off (Schwellenwert) einsetzen, würden wir ja allen, die über dem Schwellenwert sind, das gleichde Risiko zuordnen, was definitiv nicht der Fall ist. Wir sollten lernen, mit Wahrscheinlichkeiten zu arbeiten und nicht auf täuschende Positiv / Negativ Klassifikationen zu setzen.

Wir haben jetzt gesehen, dass es keinen Sinn macht, die Werte des TUGs zu dichotomisieren, um jemanden als Sturzrisiko vorhanden oder Sturzrisiko nicht vorhanden zu dichotomisieren. Aber für diese Übung werden wir es trotzdem tun. In der Literatur wird oft ein Cut Off von grösser gleich 12 Sekunden empfohlen Empfehlung CDC, hier klicken für das PDF.

Wir können also jetzt Stürzer mit einem Wert von 12 oder mehr Sekunden im Tug als Test-Positive klassifizieren.

df<-df %>% 
  mutate(Test_positive=ifelse(TUG>=12,1,0))

summarytools::ctable(df$Test_positive,df$Faller)

Error in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, : invalid value 0 for 'digits' argument

Nun erstellen wir noch eine klassische diagnostische Tabelle (auf englisch auch confusion matrix genannt).

(#fig:Confusion_Matrix)Diagnostische Tabelle, auch Confusion Matrix genannt.

Abbildung 38.1: Berechnung der Sensitivität.

Abbildung 38.2: Berechnung der Spezifität

Abbildung 38.3: Berechnung des positiven prädiktiven Wertes

Abbildung 38.4: Berechnung des negativen prädiktiven Wertes

Hier möchten wir in der ersten Zeile die Test positiven Personen haben und in der ersten Spalte die, die wirklich Stürzen werden.

df<-df %>% 
  mutate(Test_positive=factor(Test_positive, levels=c(1,0), labels=c("Test +", "Test -")),
         Faller=factor(Faller, levels=c("1","0"), labels=c("Stürzer", "Nicht-Stürzer")))

summarytools::ctable(df$Test_positive,df$Faller)

Cross-Tabulation, Row Proportions  
Test_positive * Faller  
Data Frame: df  

--------------- -------- ------------- --------------- --------------
                  Faller       Stürzer   Nicht-Stürzer          Total
  Test_positive                                                      
         Test +             81 (47.6%)      89 (52.4%)   170 (100.0%)
         Test -             97 (47.5%)     107 (52.5%)   204 (100.0%)
          Total            178 (47.6%)     196 (52.4%)   374 (100.0%)
--------------- -------- ------------- --------------- --------------

Wir könnten jetzt selber die Sensitivität und Spezifität ausrechnen, doch da das manuelle ausrechnen der Konfidenzinveralle etwas mühsam wäre, benutzen wir das epiR Paket.

tabelle_diag <- table(df$Test_positive, df$Faller)

results <- epi.tests(tabelle_diag, method = "exact", digits = 2, 
   conf.level = 0.95)
print(results)

          Outcome +    Outcome -      Total
Test +           81           89        170
Test -           97          107        204
Total           178          196        374

Point estimates and 95% CIs:
--------------------------------------------------------------
Apparent prevalence *                  0.45 (0.40, 0.51)
True prevalence *                      0.48 (0.42, 0.53)
Sensitivity *                          0.46 (0.38, 0.53)
Specificity *                          0.55 (0.47, 0.62)
Positive predictive value *            0.48 (0.40, 0.55)
Negative predictive value *            0.52 (0.45, 0.59)
Positive likelihood ratio              1.00 (0.80, 1.25)
Negative likelihood ratio              1.00 (0.83, 1.20)
False T+ proportion for true D- *      0.45 (0.38, 0.53)
False T- proportion for true D+ *      0.54 (0.47, 0.62)
False T+ proportion for T+ *           0.52 (0.45, 0.60)
False T- proportion for T- *           0.48 (0.41, 0.55)
Correctly classified proportion *      0.50 (0.45, 0.55)
--------------------------------------------------------------
* Exact CIs

Eine Sensitivität oder eine Spezifität von 0.5 bedeutet eigentlich, dass wir auch eine Münze werfen könnten. wir sehen also, dass der dichotomisierte TUG nicht eine wahnsinnig gute Information bietet.

Wir können jetzt aber auch direkt von unserer logistischen Regression eine Receiver Operating Characteristic Curve erstellen und so einen optimalen Cut Off bestimmen (obschon wir ja eigentlich wissen, dass es diesen nicht wirklich gibt).

Wie meistens in R gibt es tausend Lösungen und fast so viele Pakete für ein Problem. Wir werden hier ein paar Varianten aufzeigen.

Zum Beispiel können wir das Paket ROCnReg.

Da die Variable Faller eine Faktor-Variable ist, müssen wir die zuerst zu einer numerischen Variable umwandeln für dieses Paket.

df<-data_script[,c("id", "Faller", "TUG")]
df<-df %>% 
  mutate(Faller=as.numeric(as.character(Faller)))

m1_emp <- pooledROC.emp(marker = "TUG", group = "Faller",
tag.h = 0, data = df, p = seq(0,1,l=101), B = 10,
method = "coutcome", pauc = pauccontrol(compute = TRUE, value = 0.5, focus = "FPF"))
summary(m1_emp)

Call:
pooledROC.emp(marker = "TUG", group = "Faller", tag.h = 0, data = df, 
    p = seq(0, 1, l = 101), B = 10, method = "coutcome", pauc = pauccontrol(compute = TRUE, 
        value = 0.5, focus = "FPF"))

Approach: Pooled ROC curve - Empirical
----------------------------------------------
Area under the pooled ROC curve: 0.524 (0.46, 0.578)*
Partial area under the pooled ROC curve (FPF = 0.5): 0.281 (0.216, 0.351)*
 * Confidence level:  0.95

Sample sizes:
                           Group H     Group D
Number of observations         196         178
Number of missing data           0           0

plot(m1_emp)

Abbildung 38.5: Receiver Operating Characteristic Curve (ROC-Curve). AUC = Area Under the Curve, Fläche unter der Kurve. AUC von 0.5 bedeutet, dass man eigentlich eine Münze werfen könnte. Die AUC kann als Wahrscheinlichkeit interpretiert werden, dass bei einem zufällig ausgewählten Stürzer und einem zufällig ausgewählten Nicht-Stürzer, der Stürzer eine höhere vorhergesagte Sturzwahrscheinlichkeit hätte. TPF = True Positive Fraction, FPF = False Positive Fraction.

Ja, die ROC-Kurve für den TUG sieht nicht sehr gut aus. Doch wie würde eine besser ROC-Kurve aussehen? Ein sehr guter Prädiktor, ist die Variable, ob die Patient:innen in den letzten zwölf Monaten vor der Prädiktion schon gestürzt waren. Wir benutzen diesmal ein anderes Paket (zur Illustration der Diversität von R), nämlich das pRoc Paket.

df<-data_script[,c("id", "Faller", "N_Falls_Past_12Months_prior_t0")]
df<-df %>% 
  mutate(Faller=as.numeric(as.character(Faller)))

roc2<-pROC::roc(df$Faller, df$N_Falls_Past_12Months_prior_t0)

Setting levels: control = 0, case = 1

Setting direction: controls < cases

ggroc(roc2, legacy.axes = TRUE)+
  geom_abline(intercept=0, slope=1)+
  annotate("text", x = 0.25, y = 0.55, label = paste0("AUC: ", round(roc2$auc,3)))

Abbildung 11.1: ROC-Kurve für den Prädiktor Anzahl Stürze in der Vergangenheit.

Ein fast perfekter Prädiktor hätte eine ROC-Kurve, die in etwa so aussehen würde (hier müssen wir etwas schummeln und Daten erfinden): Wir erstellen eine Variable fast_perfekter_prädiktor, die nur in ungefähr 5% der Fälle falsch liegt.

df<-df %>% 
  mutate(fast_perfekter_prädiktor=case_when(
    Faller==0~sample(c(0,1),size=length(id), replace=TRUE, prob=c(0.95,0.05)), 
    Faller==1~sample(c(0,1),size=length(id), replace=TRUE, prob=c(0.05,0.95))))

roc3<-pROC::roc(df$Faller, df$fast_perfekter_prädiktor,  direction="<")

Setting levels: control = 0, case = 1

ggroc(roc3, legacy.axes = TRUE)+
  geom_abline(intercept=0, slope=1)+
  annotate("text", x = 0.25, y = 0.55, label = paste0("AUC: ", round(roc3$auc,3)))

So zum Spass können wir ja auch noch einen Prädiktor anschauen, der schlechter ist als eine Münze zu werfen. Wir erfinden einen Prädiktor, der nur in etwa 5% der Fälle richtig liegt:

df<-df %>% 
  mutate(nutloser_prädiktor=case_when(
    Faller==1~sample(c(0,1),size=length(id), replace=TRUE, prob=c(0.95,0.05)), 
    Faller==0~sample(c(0,1),size=length(id), replace=TRUE, prob=c(0.05,0.95))))

roc4<-pROC::roc(df$Faller, df$nutloser_prädiktor, direction="<")

Setting levels: control = 0, case = 1

ggroc(roc4, legacy.axes = TRUE)+
  geom_abline(intercept=0, slope=1)+
  annotate("text", x = 0.25, y = 0.55, label = paste0("AUC: ", round(roc4$auc,3)))

Achtung: Wenn Sie in einer Analyse oder inem veröffentlichten Artikel eine solche Kurve sehen, ist meistend etwas schief gelaufen, zum Beispiel, dass die Forscher:innen die Richtung des Testes falsch definiert haben. Wir müssen ja angeben, ob hohe Werte zu einem höheren Risiko führen oder tiefere Werte ein höheres Risiko aufweisen. Macht man dies falsch, so geht die ROC-Kurve gerade in die falsche Richtung (d.h. nach unten).

Gehen wir wieder zurück zum TUG. Wir könnten jetzt mit unserer ROC-Kurve noch versuchen, einen optimalen Cut Off zu bestimmen. Dazu gibt es folgende Überlegunen:

• Ist Sensitivität und Spezifität - respektive falsch negative Testresultate und falsch positive Testresultate - gleich wichtig?
• Ist Sensitivität wichtiger?
• Ist Spezifität wichtiger?
• Andere Kriterien?

Mit dem Paket OptimalCutpoints können wir Kriterien auswählen, mit denen wir den Cut Off bestimmen.

Ein weiteres Paket, ist das cutpointr Paket

df<-data_script[,c("Faller", "TUG")]
df<-df %>% 
  mutate(Faller=as.numeric(as.character(Faller)))

Paket optimalCutpoints

• Achtung: Bei mir lief der Befehl optimalCutpoints nicht, ich bekam folgenden Fehler (siehe unterhalb des Code-Chunks –> und vor allem die Lösung etwas weiter unten! 😄)

optimal.cutpoint.Youden<-optimal.cutpoints(X = "TUG", status = "Faller", tag.healthy = 0, 
                                           methods = "Youden", data = df)

Warning in xtfrm.data.frame(x): cannot xtfrm data frames

Error in `x[order(x, na.last = na.last, decreasing = decreasing)]`:
! Can't subset columns past the end.
ℹ Locations 229, 194, 215, …, 301, and 221 don't exist.
ℹ There is only 1 column.

Google hatte auch hier eine Lösung, manchmal gibt es Probleme mit dem Datenformat. Unser df ist ein tibble (da wir es mit dem Paket dplyr erstellt haben). Nun müssen wir es ganz einfach in ein echtes data.frame umwandeln. Das geht mit folgendem Befehl: df <- as.data.frame(df). Wir berechnen zuerst den oft benutzten Youden-Index, um den Cut-Off zu bestimmen (der Wert, bei dem Sens + Spez -1 am höchsten ist).

df <- as.data.frame(df)
optimal.cutpoint.Youden<-optimal.cutpoints(X = "TUG", status = "Faller", tag.healthy = 0, 
                                           methods = "Youden", data = df)

summary(optimal.cutpoint.Youden)

Call:
optimal.cutpoints.default(X = "TUG", status = "Faller", tag.healthy = 0, 
    methods = "Youden", data = df)

Area under the ROC curve (AUC):  0.524 (0.466, 0.583) 

CRITERION: Youden
Number of optimal cutoffs: 1

                      Estimate
cutoff             15.29000000
Se                  0.27528090
Sp                  0.81122449
PPV                 0.56976744
NPV                 0.55208333
DLR.Positive        1.45824476
DLR.Negative        0.89336443
FP                 37.00000000
FN                129.00000000
Optimal criterion   0.08650539

Wir könnten jedoch auch sagen, dass wir den Cut Off möchten, bei dem wir eine bestimmte Likelihood-Ratio für einen negativen Test haben.

control.cutpoints(valueDLR.Negative = 0.5)

$costs.ratio
[1] 1

$CFP
[1] 1

$CFN
[1] 1

$valueSp
[1] 0.85

$valueSe
[1] 0.85

$maxSp
[1] TRUE

$generalized.Youden
[1] FALSE

$costs.benefits.Youden
[1] FALSE

$costs.benefits.Efficiency
[1] FALSE

$weighted.Kappa
[1] FALSE

$standard.deviation.accuracy
[1] FALSE

$valueNPV
[1] 0.85

$valuePPV
[1] 0.85

$maxNPV
[1] TRUE

$valueDLR.Positive
[1] 2

$valueDLR.Negative
[1] 0.5

$adjusted.pvalue
[1] "PADJMS"

$ci.SeSp
[1] "Exact"

$ci.PV
[1] "Exact"

$ci.DLR
[1] "Transformed"

optimal.cutpoint.minimiseLR_negative<-optimal.cutpoints(X = "TUG", status = "Faller", tag.healthy = 0, 
                                           methods = "ValueDLR.Negative", control = control.cutpoints(), data = df)

Warning: There is no cutpoint that yields the exact Diagnostic Negative 
 Likelihood Ratio designated. The cutpoint having the closest value to the 
 designated Diagnostic Negative Likelihood Ratio has therefore been selected.

summary(optimal.cutpoint.minimiseLR_negative)

Call:
optimal.cutpoints.default(X = "TUG", status = "Faller", tag.healthy = 0, 
    methods = "ValueDLR.Negative", data = df, control = control.cutpoints())

Area under the ROC curve (AUC):  0.524 (0.466, 0.583) 

CRITERION: ValueDLR.Negative
Number of optimal cutoffs: 1

                 Estimate
cutoff         5.85000000
Se             0.99438202
Sp             0.01020408
PPV            0.47708895
NPV            0.66666667
DLR.Positive   1.00463338
DLR.Negative   0.55056180
FP           194.00000000
FN             1.00000000

Zur Illustration hier noch ein weitere Paket, das oben erwähnte cutpointr. Paket cutpointr

cp <- cutpointr(df, x=TUG, class=Faller, 
                method = maximize_metric, metric = youden)

Assuming the positive class is 1

Assuming the positive class has higher x values

summary(cp)

Method: maximize_metric 
Predictor: TUG 
Outcome: Faller 
Direction: >= 

    AUC   n n_pos n_neg
 0.5243 374   178   196

 optimal_cutpoint youden    acc sensitivity specificity tp  fn fp  tn
            15.29 0.0865 0.5561      0.2753      0.8112 49 129 37 159

Predictor summary: 
    Data  Min.     5% 1st Qu.  Median     Mean  3rd Qu.      95%   Max.
 Overall 5.345 6.9030  9.2775 11.6400 12.78533 14.44250 23.40225 40.805
       0 5.345 7.0075  9.1400 11.6100 12.47714 14.19250 21.65125 34.960
       1 5.620 6.8000  9.5275 11.6925 13.12469 15.75375 24.11575 40.805
       SD NAs
 5.241816   0
 4.878775   0
 5.608807   0

38.3 Beispiel eines etwas besseren Prädiktors, nämlich der Anzahl Stürze in den letzten zwölf Monaten vor dem Assessment.

roc2<-pROC::roc(data_script$Faller, data_script$N_Falls_Past_12Months_prior_t0)

Setting levels: control = 0, case = 1

Setting direction: controls < cases

ggroc(roc2, legacy.axes = TRUE)+
  geom_abline(intercept=0, slope=1)+
  annotate("text", x = 0.25, y = 0.55, label = paste0("AUC: ", round(roc2$auc,3)))+
  theme_classic()

Abbildung 4.8: ROC-Kurve der Anzahl Stürze in den zwölf Monaten vor dem Assessment. Im Vergleich zum Tug schneidet diese simple Frage “Wie oft sind sie in den letzten zwölf Monaten gestürzt” gar nicht so schlecht ab.

Zuerst erstellen wir Gruppen basiert auf dem TUG, wir erstellen alle 5 Sekunden im TUG einen neue Gruppe

df<-data_script %>% 
    mutate(
        group_N_Falls_prior_t0 = cut_width(N_Falls_Past_12Months_prior_t0, width=3),
        ) 

plot(df$N_Falls_Past_12Months_prior_t0, df$group_N_Falls_prior_t0)

Abbildung 38.6: Anzahl Stürze in den letzten zwölf Monaten aufgeteilt in Gruppen.

Wahrscheinlichkeit, dass jemand innerhalb einer Gruppe gestürzt ist:

Dank an den Autor vom Cookbook-r

# New version of length which can handle NA's: if na.rm==T, don't count them
    length2 <- function (x, na.rm=FALSE) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
    }

df_summary<-df %>% 
  group_by(group_N_Falls_prior_t0) %>% 
  summarise(Freq = length2(group_N_Falls_prior_t0, na.rm=TRUE),
            empirical_prob_fall=sum(as.numeric(as.character(Faller)), na.rm=TRUE)/sum(Freq))

ggplot(df_summary, aes(x=group_N_Falls_prior_t0, y=empirical_prob_fall))+
  geom_point(size=4, color="#830CF5")+
  scale_y_continuous(limits = c(0.0, 1), breaks = seq(0, 1, by = 0.1))+
  theme_classic()

Abbildung 38.7: Empirische Wahrscheinlichkeit in den nächsten 12 Monaten zu stürzen in Gruppen nach Anzahl der Stürze in den letzten zwölf Monaten. Wir sehen, dass die Wahrscheinlichkeit stetigt steigt, bis 8 Stürze in der Vergangenheit. Das ist jetzt aber kein Argument, dort einen Cut-Off zu setzen, da das Risiko zuvor konstant ansteigt. Würden wir einen Cut-off (Schwellenwert) einsetzen, würden wir ja allen, die unter dem Schwellenwert sind, kein Risiko zuordnen, was definitiv falsch negativ wäre.

38.3.1 smoothed plot

Wir haben oben die Graphik für das Sturzrisiko über Gruppen gerechnet. Gruppen anhand einer kontinuierlichen Variable erstellen, ist eigentlich nicht die optimale Vorgehensweise. Eigentlich sollten wir eine geglättete (smoothed) Linie zwischen dem Sturzrisiko und der kontinuierlichen Variable zeichnen. Die nächste Graphik zeigt dies. Normalerweise benutzen wir für die logistische Regression den Befehl glm, damit wir jedoch die folgende Graphik erstellen können, benutzen wir den Befehl lrm des Pakets rms von Frank Harrell. Siehe auch hier

lr.fit_wie_wir_es_normalerweise_tun_würden<-glm(Faller~N_Falls_Past_12Months_prior_t0, data=data_script, family="binomial") # So würden wir sonst eine logistische Regression in R rechnen. 

lr.fit<-rms::lrm(Faller~N_Falls_Past_12Months_prior_t0, data=data_script) # Diese Variante führen wir durch, damit wir die folgende Graphik einfach zeichnen können. 

ggplot(Predict(lr.fit,N_Falls_Past_12Months_prior_t0=0:30, fun=plogis))+
  labs(y="Probability of a Fall")+
  theme_classic()+
  scale_x_continuous(breaks=seq(0,30, by=2))

Error in Predict(lr.fit, N_Falls_Past_12Months_prior_t0 = 0:30, fun = plogis): could not find function "Predict"

Ein anderes Beispiel ist die Diagnose von Schulterproblemen.