Die Rohdaten der Studie mit dem Titel “Test-retest reliability and minimal detectable change scores for the short physical performance battery, one-legged standing test and timed up and go test in patients undergoing hemodialysis” finden Sie auf Moodle. Auf den Volltext können Sie hier zugreifen.

Übung 1

Aufgaben

  1. Laden Sie die Daten in R. Splitten Sie den Datensatz auf nach Test. Erstellen Sie dazu jeweils neue R-Objekte. Bringen Sie dann die drei Teildatensätze in ein long-format. Sie können für den zweiten Schritt folgenden Code verwenden:
library(tidyverse)
df.long <- df %>% 
  pivot_longer(cols = c(X,Y), values_to = "value_name", names_to = "time_point")

  1. Reproduzieren Sie die angegebenen SEM in Tabelle 4. Verwenden Sie das Paket SimplyAgree (siehe Slides für den R-Code).

  2. Ohne auf die Details des ICCs (Tabelle 3) einzugehen: Haben die Autor:innen den SEM mit oder ohne Berücksichtigung eines systematischen Fehlers berechnet?

Lösungen

Hier wird der Datensatz importiert, gesplitted nach Test und dann jeweils ins long-format gebracht. Natürlich gibt es auch andere Möglichkeiten, dies zu tun.

library(rio)
df.tests <- import("../data/data-villar.csv")

library(tidyverse)
df.sppb <- df.tests %>% 
  select(ID, SPPB_A, SPPB_B)

df.sppb.long <- df.sppb %>% 
  pivot_longer(cols = c(SPPB_A, SPPB_B), values_to = "value", names_to = "time")

df.olst <- df.tests %>% 
  select(ID, OLST_A, OLST_B)

df.olst.long <- df.olst %>% 
  pivot_longer(cols = c(OLST_A, OLST_B), values_to = "value", names_to = "time")

df.tug <- df.tests %>% 
  select(ID, TUG_A, TUG_B)

df.tug.long <- df.tug %>% 
  pivot_longer(cols = c(TUG_A, TUG_B), values_to = "value", names_to = "time")

Die drei SEM erhalten wir wie folgt:

SPPB:

library(SimplyAgree)

reli.sppb <- reli_stats(
  measure = "value",   # name of the column with ratings
  item = "time",      # rater or time identifier
  id = "ID",      # subject identifier
  data = df.sppb.long, # your data frame
  se_type = "ICC2",   # this is important!
  cv_calc = "MSE",     # how to calculate CV
  conf.level = 0.95,
  other_ci = TRUE
)

reli.sppb
## 
## Coefficient of Variation (%): 7.04 95% C.I. [6.04, 8.45]
## Standard Error of Measurement (SEM): 0.71 95% C.I. [0.61, 0.851]
## Standard Error of the Estimate (SEE): 0.69 95% C.I. [0.592, 0.827]
## Standard Error of Prediction (SEP): 0.99 95% C.I. [0.85, 1.19]
## 
## Intraclass Correlation Coefficients
##            Model         Measures  Type    ICC Lower CI Upper CI
## 1 one-way random        Agreement  ICC1 0.9442   0.9183   0.9620
## 2 two-way random        Agreement  ICC2 0.9442   0.9178   0.9622
## 3  two-way fixed      Consistency  ICC3 0.9459   0.9206   0.9632
## 4 one-way random   Avg. Agreement ICC1k 0.9713   0.9574   0.9806
## 5 two-way random   Avg. Agreement ICC2k 0.9713   0.9571   0.9808
## 6  two-way fixed Avg. Consistency ICC3k 0.9722   0.9587   0.9813

Auf Grund der Beschreibung in der Methode ist anzunhemen, dass die Autor:innen den systematischen Fehler berücksichtigt haben.

“two-way random-effects model were used to assess the test–retest reliability of the data”

Mit “ICC2” erhalten Sie den SEM mit Berücksichtigung des systematischen Fehlers, mit “ICC3” ohne. Die Resultate zwischen der Studie und hier sind ähnlich, aber nicht gleich.

Die manuelle Berechnung mit der Formel (wie in der Methodik angegeben)

\[ SEM = SD * \sqrt{1-R} \]

führt zu einem ähhnlichen Resultat:

sd(c(df.sppb$SPPB_A, df.sppb$SPPB_B), na.rm = TRUE) * sqrt(1-0.9442)
## [1] 0.6904926

Diese Formel sollte jedoch für den \(SEM_{agreement}\) nicht verwendet werden.

Die Varianzkomponenten (beste Methode) bekommen wir mit

reli.sppb$var_comp
##            variance    percent
## ID       8.23756523 0.94420903
## Items    0.01521492 0.00174397
## Residual 0.47152236 0.05404700
## Total    8.72430251 1.00000000

Hier wäre der SEM

sqrt(0.01521492 + 0.47152236)
## [1] 0.6976656


OLST:

library(SimplyAgree)

reli.olst <- reli_stats(
  measure = "value",   # name of the column with ratings
  item = "time",      # rater or time identifier
  id = "ID",      # subject identifier
  data = df.olst.long, # your data frame
  se_type = "ICC2",   # this is important!
  cv_calc = "MSE",     # how to calculate CV
  conf.level = 0.95,
  other_ci = TRUE
)

reli.olst
## 
## Coefficient of Variation (%): 34.3 95% C.I. [28.9, 42.4]
## Standard Error of Measurement (SEM): 4.99 95% C.I. [4.27, 6.01]
## Standard Error of the Estimate (SEE): 4.71 95% C.I. [4.03, 5.67]
## Standard Error of Prediction (SEP): 6.86 95% C.I. [5.87, 8.26]
## 
## Intraclass Correlation Coefficients
##            Model         Measures  Type    ICC Lower CI Upper CI
## 1 one-way random        Agreement  ICC1 0.8905   0.8405   0.9255
## 2 two-way random        Agreement  ICC2 0.8905   0.8405   0.9255
## 3  two-way fixed      Consistency  ICC3 0.8905   0.8403   0.9256
## 4 one-way random   Avg. Agreement ICC1k 0.9421   0.9134   0.9613
## 5 two-way random   Avg. Agreement ICC2k 0.9421   0.9134   0.9613
## 6  two-way fixed Avg. Consistency ICC3k 0.9421   0.9132   0.9614

Gleiche Interpretation wie oben.

Mit der angegebenen Formel (für \(SEM_{agreement}\) nicht ideal):

sd(c(df.olst$OLST_A, df.olst$OLST_B), na.rm = TRUE) * sqrt(1-0.8905)
## [1] 4.928476

Varianzkomponenten:

reli.olst$var_comp
##           variance  percent
## ID       194.95939 0.890523
## Items      0.00000 0.000000
## Residual  23.96746 0.109477
## Total    218.92685 1.000000
sqrt(23.96746)
## [1] 4.895657

Es gibt hier überhaupt keinen systematischen Fehler (Between-Rater Varianz = 0). Daher sind \(SEM_{agreement}\) und \(SEM_{consistency}\) gleich.


TUG:

library(SimplyAgree)

reli.tug <- reli_stats(
  measure = "value",   # name of the column with ratings
  item = "time",      # rater or time identifier
  id = "ID",      # subject identifier
  data = df.tug.long, # your data frame
  se_type = "ICC2",   # this is important!
  cv_calc = "MSE",     # how to calculate CV
  conf.level = 0.95,
  other_ci = TRUE
)

reli.tug
## 
## Coefficient of Variation (%): 10.6 95% C.I. [9.11, 12.8]
## Standard Error of Measurement (SEM): 1.18 95% C.I. [1.01, 1.41]
## Standard Error of the Estimate (SEE): 1.16 95% C.I. [0.995, 1.39]
## Standard Error of Prediction (SEP): 1.65 95% C.I. [1.42, 1.98]
## 
## Intraclass Correlation Coefficients
##            Model         Measures  Type    ICC Lower CI Upper CI
## 1 one-way random        Agreement  ICC1 0.9644   0.9476   0.9758
## 2 two-way random        Agreement  ICC2 0.9644   0.9476   0.9759
## 3  two-way fixed      Consistency  ICC3 0.9645   0.9477   0.9760
## 4 one-way random   Avg. Agreement ICC1k 0.9819   0.9731   0.9878
## 5 two-way random   Avg. Agreement ICC2k 0.9819   0.9731   0.9878
## 6  two-way fixed Avg. Consistency ICC3k 0.9819   0.9731   0.9878

Gleiche Interpretation wie oben.

Mit der angegebenen Formel (für \(SEM_{agreement}\) nicht ideal):

sd(c(df.tug$TUG_A, df.tug$TUG_B), na.rm = TRUE) * sqrt(1-0.9644)
## [1] 1.14371

Varianzkomponenten:

reli.tug$var_comp
##            variance      percent
## ID       36.6500366 0.9643688872
## Items     0.0049599 0.0001305094
## Residual  1.3491709 0.0355006034
## Total    38.0041674 1.0000000000
sqrt(0.0049599 + 1.3491709)
## [1] 1.163671

Hier ist die Abweichung vom Resultat in der Studie eher gross.

PS: Das 95% CI in der Publikation sieht etwas seltsam aus…

Übung 2

Aufgaben

  1. Reproduzieren Sie auch die angegebenen SDCs.
  2. Rechnen Sie zusätzlich einen \(SDC_{95}\)
  3. Interpretieren Sie die Werte.

Lösungen

Wir bestimmen zuerst die z-Werte:

z.9 <- qnorm(0.95)
z.95 <- qnorm(0.975)

Damit können wir nun die SDCs bestimmen. Hier werden die SEM-Werte aus der Studie verwendet, man könnte aber auch jene von oben nehmen:

(SDC.9.sppb <- z.9 * sqrt(2) * 0.72)
## [1] 1.674846
(SDC.95.sppb <- z.95 * sqrt(2) * 0.72)
## [1] 1.995702
(SDC.9.olst <- z.9 * sqrt(2) * 4.82)
## [1] 11.21216
(SDC.95.olst <- z.95 * sqrt(2) * 4.82)
## [1] 13.36011
(SDC.9.tug <- z.9 * sqrt(2) * 1.24)
## [1] 2.884456
(SDC.95.tug <- z.95 * sqrt(2) * 1.24)
## [1] 3.437041

Wenn wir mit dem entsprechenden Assessment bei einer Person der entsprechenden Population eine Differenz von z.B. 1.7 (SPPB) messen, können wir zu 90% sicher sein, dass es eine wahre Veränderung und nicht nur Messfehelr ist. (bzw. 95% sicher sein, wenn die Differenz \(\ge 2\) ist. Analog für die anderen Assessments.

Übung 3

Aufgaben

  1. Untersuchen Sie für das Assessment SPPB den systematischen Fehler. Berechnen Sie dafür zuerst die durchschnittliche Differenz zwischen den zwei Messzeitpunkten.

  2. Erstellen Sie einen Bland-Altman Plot mit mit einem Agreement-Level von 95% (Syntax, siehe Slides). Vergleichen Sie mit Abbildung 2 aus der Publikation.

  3. Interpretieren Sie den Bias und die LoA.

Lösungen

Die durchschnittliche Differenz berechnet sich wie folgt:

mean(df.sppb$SPPB_A - df.sppb$SPPB_B, na.rm = TRUE)
## [1] -0.2

Für die Bland-Altman-Analyse verwenden wir das Paket “SimplyAgree”:

ba.sppb <- agree_test(x = df.sppb$SPPB_A,
                 y = df.sppb$SPPB_B,
                agree.level = .95)

plot(ba.sppb)

Der Plot sieht genau gleich aus (schöner 🤓) wie jener in der Studie.

Der Analyse entnehmen wir den oben schon manuell berechneten Bias. Da das 95% CI die 0 knapp beinhaltet, unterscheidet sich dieser nicht signifikant von 0.

print(ba.sppb)
## Limit of Agreement = 95%
## 
## ###- Shieh Results -###
## Exact 90% C.I.  [-2.3423, 1.9423]
## Hypothesis Test: No Hypothesis Test
## 
## ###- Bland-Altman Limits of Agreement (LoA) -###
##           Estimate Lower CI Upper CI CI Level
## Bias        -0.200  -0.4407  0.04072     0.95
## Lower LoA   -2.104  -2.4495 -1.75858     0.90
## Upper LoA    1.704   1.3586  2.04950     0.90
## 
## ###- Concordance Correlation Coefficient (CCC) -###
## CCC: 0.9416, 95% C.I. [0.9063, 0.9638]

Das könnte man auch mit einem einfachen t-Test erfahren:

t.test(df.sppb$SPPB_A - df.sppb$SPPB_B)
## 
##  One Sample t-test
## 
## data:  df.sppb$SPPB_A - df.sppb$SPPB_B
## t = -1.6598, df = 64, p-value = 0.1018
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.44071797  0.04071797
## sample estimates:
## mean of x 
##      -0.2

Warum ist in der Studie ein P-Wert von 0.942 angegeben? Was ist in dieser Tabelle schief gelaufen? (siehe hier).