I. Einführung und Zielsetzung

In der klinischen Forschung ist es oft notwendig, nicht direkt messbare Konstrukte (z. B. „Schweregrad der Fatigue“) zu quantifizieren.
Solche latenten Variablen werden in Studien häufig über mehrere Fragebogenskalen oder klinische Tests erfasst.

Problemstellung:
In einer Pilotstudie zu einer seltenen Erkrankung wurden nur 21 Patient:innen untersucht.
Es liegen vier klinische Scores vor, die allesamt den gleichen zugrunde liegenden Schweregrad (Faktor Disease Severity) messen sollen.
Die klassische Maximum-Likelihood-CFA scheitert hier oft wegen der kleinen Stichprobe.

Ziel dieser Fallstudie:
Wir wollen untersuchen, ob eine Bayesianische CFA/SEM in der Lage ist, bei sehr kleiner Stichprobe stabile Schätzungen zu liefern,
und vergleichen dafür die beiden Implementierungen bcfa() und bsem() aus dem R-Paket blavaan.


II. Methodik

II.1. Studiendesign und Datensimulation

Obwohl der Kunde uns die Erlaubnis zur methodischen Veröffentlichung seines Projekts erteilt hat, wurden die Daten zu Demonstrationszwecken simuliert. Sie orientieren sich jedoch an typischen klinischen Studien zur Fatigue bei chronischen Erkrankungen.
Jede der vier Variablen entspricht einem validierten Subscore (z. B. körperliche Fatigue, mentale Fatigue, allgemeiner Gesundheitszustand, Aktivitätseinschränkung).

Wir simulieren:

  • einen latenten Faktor (Disease Severity), normalverteilt

  • vier klinische Scores mit hoher, aber nicht perfekter Korrelation zum Faktor

  • Messfehler, um realistische klinische Datenvariabilität abzubilden

library(blavaan)
library(bayesplot)
library(dplyr)

set.seed(1234)
n <- 21  # Kleine Pilotstudie
DiseaseSeverity <- rnorm(n, mean = 0, sd = 1)

# Faktorladungen: realistische Werte aus klinischen Skalen
lambda <- c(0.8, 0.75, 0.85, 0.7)
measurement_error <- matrix(rnorm(n * 4, mean = 0, sd = 0.5), nrow = n)

# Beobachtete Scores generieren
scores <- sweep(measurement_error, 2, lambda * DiseaseSeverity, '+')
colnames(scores) <- c("Fatigue_Physical", "Fatigue_Mental", "Health_Status", "Activity_Limitation")
df <- as.data.frame(scores)

II.2. Modell-Spezifikation

Wir gehen davon aus, dass alle vier Scores denselben latenten Faktor Disease Severity messen.

model <- '
  DiseaseSeverity =~ Fatigue_Physical + Fatigue_Mental + Health_Status + Activity_Limitation
'

II.3. Modellschätzung

Wir vergleichen zwei Ansätze:

  1. bcfa() – optimiert für CFA-Modelle

  2. bsem() – flexiblere SEM-Variante

fit_bcfa <- bcfa(model, data = df,
                 burnin = 500, sample = 1000, n.chains = 2,
                 target = "stan", seed = 1)
## 
## SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.016963 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 169.63 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 1500 [  0%]  (Warmup)
## Chain 1: Iteration:  150 / 1500 [ 10%]  (Warmup)
## Chain 1: Iteration:  300 / 1500 [ 20%]  (Warmup)
## Chain 1: Iteration:  450 / 1500 [ 30%]  (Warmup)
## Chain 1: Iteration:  501 / 1500 [ 33%]  (Sampling)
## Chain 1: Iteration:  650 / 1500 [ 43%]  (Sampling)
## Chain 1: Iteration:  800 / 1500 [ 53%]  (Sampling)
## Chain 1: Iteration:  950 / 1500 [ 63%]  (Sampling)
## Chain 1: Iteration: 1100 / 1500 [ 73%]  (Sampling)
## Chain 1: Iteration: 1250 / 1500 [ 83%]  (Sampling)
## Chain 1: Iteration: 1400 / 1500 [ 93%]  (Sampling)
## Chain 1: Iteration: 1500 / 1500 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.445 seconds (Warm-up)
## Chain 1:                1.42 seconds (Sampling)
## Chain 1:                2.865 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000122 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.22 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 1500 [  0%]  (Warmup)
## Chain 2: Iteration:  150 / 1500 [ 10%]  (Warmup)
## Chain 2: Iteration:  300 / 1500 [ 20%]  (Warmup)
## Chain 2: Iteration:  450 / 1500 [ 30%]  (Warmup)
## Chain 2: Iteration:  501 / 1500 [ 33%]  (Sampling)
## Chain 2: Iteration:  650 / 1500 [ 43%]  (Sampling)
## Chain 2: Iteration:  800 / 1500 [ 53%]  (Sampling)
## Chain 2: Iteration:  950 / 1500 [ 63%]  (Sampling)
## Chain 2: Iteration: 1100 / 1500 [ 73%]  (Sampling)
## Chain 2: Iteration: 1250 / 1500 [ 83%]  (Sampling)
## Chain 2: Iteration: 1400 / 1500 [ 93%]  (Sampling)
## Chain 2: Iteration: 1500 / 1500 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.6 seconds (Warm-up)
## Chain 2:                2.2 seconds (Sampling)
## Chain 2:                3.8 seconds (Total)
## Chain 2:
## Computing post-estimation metrics (including lvs if requested)...
fit_bsem <- bsem(model, data = df,
                 burnin = 500, sample = 1000, n.chains = 2,
                 target = "stan", seed = 2)
## 
## SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000109 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.09 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 1500 [  0%]  (Warmup)
## Chain 1: Iteration:  150 / 1500 [ 10%]  (Warmup)
## Chain 1: Iteration:  300 / 1500 [ 20%]  (Warmup)
## Chain 1: Iteration:  450 / 1500 [ 30%]  (Warmup)
## Chain 1: Iteration:  501 / 1500 [ 33%]  (Sampling)
## Chain 1: Iteration:  650 / 1500 [ 43%]  (Sampling)
## Chain 1: Iteration:  800 / 1500 [ 53%]  (Sampling)
## Chain 1: Iteration:  950 / 1500 [ 63%]  (Sampling)
## Chain 1: Iteration: 1100 / 1500 [ 73%]  (Sampling)
## Chain 1: Iteration: 1250 / 1500 [ 83%]  (Sampling)
## Chain 1: Iteration: 1400 / 1500 [ 93%]  (Sampling)
## Chain 1: Iteration: 1500 / 1500 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.354 seconds (Warm-up)
## Chain 1:                1.519 seconds (Sampling)
## Chain 1:                2.873 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000106 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 1500 [  0%]  (Warmup)
## Chain 2: Iteration:  150 / 1500 [ 10%]  (Warmup)
## Chain 2: Iteration:  300 / 1500 [ 20%]  (Warmup)
## Chain 2: Iteration:  450 / 1500 [ 30%]  (Warmup)
## Chain 2: Iteration:  501 / 1500 [ 33%]  (Sampling)
## Chain 2: Iteration:  650 / 1500 [ 43%]  (Sampling)
## Chain 2: Iteration:  800 / 1500 [ 53%]  (Sampling)
## Chain 2: Iteration:  950 / 1500 [ 63%]  (Sampling)
## Chain 2: Iteration: 1100 / 1500 [ 73%]  (Sampling)
## Chain 2: Iteration: 1250 / 1500 [ 83%]  (Sampling)
## Chain 2: Iteration: 1400 / 1500 [ 93%]  (Sampling)
## Chain 2: Iteration: 1500 / 1500 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.382 seconds (Warm-up)
## Chain 2:                2.232 seconds (Sampling)
## Chain 2:                3.614 seconds (Total)
## Chain 2:
## Computing post-estimation metrics (including lvs if requested)...

III. Analyse & Visualisierung

III.1. Parameterübersicht

mcmc_bcfa <- blavInspect(fit_bcfa, "mcmc")[[1]]
mcmc_bsem <- blavInspect(fit_bsem, "mcmc")[[1]]

params_bcfa <- colnames(mcmc_bcfa)
params_bsem <- colnames(mcmc_bsem)

cat("bcfa Parameters:\n")
## bcfa Parameters:
print(params_bcfa)
## [1] "DiseaseSeverity=~Fatigue_Mental"          "DiseaseSeverity=~Health_Status"          
## [3] "DiseaseSeverity=~Activity_Limitation"     "Fatigue_Physical~~Fatigue_Physical"      
## [5] "Fatigue_Mental~~Fatigue_Mental"           "Health_Status~~Health_Status"            
## [7] "Activity_Limitation~~Activity_Limitation" "DiseaseSeverity~~DiseaseSeverity"
cat("bsem Parameters:\n")
## bsem Parameters:
print(params_bsem)
## [1] "DiseaseSeverity=~Fatigue_Mental"          "DiseaseSeverity=~Health_Status"          
## [3] "DiseaseSeverity=~Activity_Limitation"     "Fatigue_Physical~~Fatigue_Physical"      
## [5] "Fatigue_Mental~~Fatigue_Mental"           "Health_Status~~Health_Status"            
## [7] "Activity_Limitation~~Activity_Limitation" "DiseaseSeverity~~DiseaseSeverity"

III.2. Traceplots der Faktorladungen

color_scheme_set("brightblue")
pars_lambda_bcfa <- params_bcfa[grepl("^lambda", params_bcfa)]
mcmc_trace(mcmc_bcfa, pars = pars_lambda_bcfa)

plot of chunk plot-all-predictorssummaryday2newplot1

color_scheme_set("viridisC")
pars_lambda_bsem <- params_bsem[grepl("^lambda", params_bsem)]
mcmc_trace(mcmc_bsem, pars = pars_lambda_bsem)

plot of chunk plot-all-predictorssummaryday2newplot1


IV. Interpretation

  • bcfa(): konzentriert sich auf reine CFA-Parameter, liefert kompaktere Posterior-Ausgaben.

  • bsem(): flexibler, schätzt zusätzliche Parameter, was zu mehr Posterior-Ausgaben führt.

  • Beide Methoden zeigen stabile Traceplots ohne Anzeichen für (sehr) schlechte Konvergenz, trotz der kleinen Stichprobe.

Klinische Relevanz:

Die latente Variable Disease Severity kann hier zuverlässig geschätzt werden, was für
Score-Validierung und klinische Entscheidungsunterstützung wichtig ist – auch wenn nur wenige Patient:innen verfügbar sind.


V. Zusammenfassung & Empfehlungen

  • Bayesianische CFA/SEM kann bei sehr kleinen Stichproben stabilere Ergebnisse liefern als klassische Methoden.

  • Für Pilotstudien in der Epidemiologie oder klinischen Forschung kann diese Methode helfen,
    latente Konstrukte wie Symptomschwere, Lebensqualität oder Funktionsfähigkeit zuverlässig zu quantifizieren.

  • bcfa() ist schlanker und fokussierter, bsem() flexibler und liefert erweiterte Modelloptionen.

  • Empfehlung: Bei n < 30 unbedingt MCMC-Diagnosen prüfen (Traceplots, R-hat, ESS)
    und wenn möglich informative Priors setzen, um klinische Vorinformationen einzubeziehen.