Fixed Effect Model

Das Fixed Effect Model geht davon aus, dass alle in die Metaanalyse einbezogenen Studien eine gemeinsame wahre Effektgröße aufweisen.

  • Summary Effect = Schätzer aller Effekte aus den Studien (\(\theta\))
  • Variation ergibt sich nur aus dem Schätzfehler innerhalb der genutzten Studien \(\rightarrow\) within-study sampling error \(\sigma^2\)
  • Gewichtung großer Studien stärker, als die kleiner Studien
  • \(H_0\) = Es gibt keinen Effekt in jeder Studie
Fixed Effect: Gewichtung:
\(T_i = \theta + e_i\) \(W_i = \frac{1}{\sigma_i^2}\)

Random Effects Model

Das Random Effects Model geht davon aus, dass die Fehler der Studien normalverteilt um \(\mu\) streuen und es verschiedene Effektgrößen gibt. \(\mu\) ist also der Mittelwert der Effektgrößenverteilung.

(Grafik: Borenstein et al. 2010 S.98)

  1. Schätzfehler innerhalb der Studien: \(\sigma^2\)
  2. Varianz zwischen den Studien: \(\tau^2\)
  • \(\tau^2\) ergibt sich beispielsweise durch verschiedene Charakteristika der Studie oder Einflussgrößen, die nicht ausgeschlossen werden können und immer mit einspielen wie beispielsweise das Geschlecht oder das Alter
  • Gewichtungen sind einheitlicher, große Studien bekommen weniger Gewicht als kleine Studien
  • \(H_0\) = Der durchschnittliche Effekt ist Null
  • \(\tau^2\) wird jedoch bei kleinem N ein schlechter Schätzer sein
Random Effect: Gewichtung:
\(T_i = \mu + u_i + e_i\) \(W_i^* = \frac{1}{\sigma_i^2+\tau^2}\)

Zu RE und FE gesamt:

  • Gewichtungen wirken sich bei beiden Modellen unterschiedlich aus:
  • bei Fixed Effect stärker
  • bei Random Effects schwächer
  • Die Entscheidung für oder gegen ein Modell sollte auf der Basis der Annahmen über die Verteilung der Effektgrößen geschehen

Ereignisdatenanalyse

Grundlegende Fragestellung: Wie lange dauert es bis ein bestimmtes Ereignis eintritt?

Dabei gibt es 2 Möglichkeiten:

  1. Zeit wird nur als Ordnungskriterium genommen (\(\rightarrow\) nicht-parametrische Modelle)
  2. Zeit wird als erklärende Variable genommen (\(\rightarrow\) parametrische oder semi-parametrische Modelle)

Mögliche Datenstrukturen

  • Rechtszensiert: im Beobachtungszeitraum ist das Ereignis noch nicht eingetreten
  • Linkszensiert: vor dem Beobachtungszeitraum ist das Ereignis eingetreten
  • Intervallzensiert: einzelne Perioden sind nicht beobachtet
  • Links/Rechtsstrunkiert: Eine Beobachtung außerhalb des Zeitraums und danach innerhalb des Zeitraums
Übersicht allgemeine Arbeitssschritte:
1.Explorativ/Deskriptiv (Datenstruktur)
2.Analytisch (Regression etc.)
3.Qualitativ (Infos über den Fall usw.)

Kaplan-Meier Schätzer

  • E = die Anzahl der Ereignisse in \(t_j\)
  • Z = die Anzahl der zensierten Episoden
  • R = die Anzahl der Episoden, die zum Zeitpuhnkt \(t_j\) im Risikobereich sind

Daraus ergibt sich folgender Schätzer:

\(\hat{S}(t) = \prod(1-\frac{E_t}{R_t})\)

Hazardfunktion

Misst die Eintrittswahrscheinlichkeit eines Ereignisses bei Fall i zum Zeitpunkt t, gegeben, dass dieses nicht vorher eingetreten ist

Bei diskreter Zeitvariable:

  • Hazard rate: \(h(t_i)=Pr(T_{ij}=j|T_i \geq j)\) und \(0 \leq h(t_{ij}) \leq 1\)
  • Die Grundfunktion ist: \(logit \ h(t_{ij}) = [ \alpha _{1} D_{1ij}+ \alpha _{2} D _{2ij}+...+\alpha _{J} D _{Jij}] + [ \beta _{1} X _{1ij}+ \beta _{2} X _{2ij}+...+ \beta _{P} X _{Pij}]\)
  • jeder Parameter des Y-Achsenabschnitts (\(\alpha\)) repräsentiert den Wert der logit hazard zum Zeitpunkt \(t_{ij}\)
  • der erste Teil der Funktion wird auch baseline logit hazard function genannt
  • jeder Steigungsparameter (\(\beta\)) ist als Effekt einer Steigung des Prädiktors um 1 zu interpretieren

(Grafik: Singer and Willett 2003 S.374)

Beispiel aus Broström (2012); Erstellt durch die Befehle (Die roten Striche um das Survival plot sind die Konfidenzintervalle):

par(mfrow=c(1,2)) #PLOTGRID
with(mort,plot(Surv(enter,exit,event),fn="cum"))
with(mort,plot(Surv(enter,exit,event),fn="surv"))

mort ist der Datensatz, enter,exit,event sind die EHA-Variablen

Bei kontinuierlicher Zeitvariable:

  • weiterhin Survivor Function: \(S(t{i_j}) = Pr(Ti \geq j)\)
  • Schätzung ist schwieriger als im diskreten Zeitfall (Singer and Willett 2003 S.491-497)
  • Da wir eine kontinuierliche Zeitvariable haben, ist die Wahrscheinlichkeit anders zu interpretieren. Da es unendlich viele Zeitpunkte in der Analyse gibt, müssen die Wahrscheinlichkeiten also gegen 0 gehen, je genauer das Zeitmaß wird

Nicht-parametrische Schätzung der Kaplan-Meier Schätzer und der Cox-Regression

Hier wird bei Broström (2012) eine Zufallsauswahl von 50 Fällen aus dem größeren mort Datensatz genommen. Man kann also den Schritt der deskriptiven Darstellung der KM-estimates bzw. der survival estimates in einer Regression weiterführen.

indx <- sample(NROW(mort),size=50,replace=FALSE)
# Ein neues Objekt namens indx wird erstellt, die Zufallsauswahl
#der Fälle geschieht durch NROW, wir sagen, dass wir 50 
#Fälle wollen, wollen das aber im eigentlichen Datensatz nicht so 
#verändert haben.
rsa <- mort[indx,]
fit <- coxph(Surv(enter,exit,event)~1,data=rsa)
s.fit <- survfit(fit)
# Cox regression model wird angewandt und durch survfit 
#kommen unsere survival estimates dazu
options(width=70)
summary(s.fit)
## Call: survfit(formula = fit)
## 
##   time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   1.21     38       1    0.974  0.0256        0.925        1.000
##   8.34     41       1    0.951  0.0341        0.886        1.000
##  10.52     40       1    0.927  0.0405        0.851        1.000
##  11.08     38       1    0.903  0.0461        0.817        0.998
##  11.25     37       1    0.879  0.0508        0.785        0.984
##  11.77     37       1    0.855  0.0545        0.755        0.969
##  12.44     36       1    0.832  0.0579        0.726        0.954
##  14.93     36       1    0.809  0.0606        0.699        0.937
##  17.52     36       1    0.787  0.0629        0.673        0.920
##  17.59     35       1    0.765  0.0649        0.648        0.903
##  18.33     34       1    0.743  0.0667        0.623        0.886
##  18.33     33       1    0.721  0.0683        0.598        0.868
##  18.50     32       1    0.698  0.0697        0.574        0.849
##  18.79     31       1    0.676  0.0709        0.551        0.831
##  19.85     30       1    0.654  0.0720        0.527        0.812

Parametrische Schätzung

par(mfrow=c(1,2))
fit.w <- phreg(Surv(enter,exit,event)~1,data=mort)
# Beim Befehl phreg wird per default der Weibull Schätzer ausgegeben, andere
# Schätzer kann man sich mit dem Argument "dist=" ausgeben lassen
plot(fit.w, fn=c("cum"))
plot(fit.w, fn=c("sur"))

Interpretation der Hazard-Rate

Siehe Tosun 2018:

  • HazardRatio = 1: Es gibt keinen Unterschied zwischen zwei Beobachtungen
  • HazardRatio > 1: Eine Beobachtung hat eine größere Wahrscheinlichkeit in einer bestimmten Zeit ein Ereignis zu erfahren; Beispiel: HR von 2,3 bedeutet, dass die Wahrscheinlichkeit um 130% erhöht ist [(2,3-1) *100 = 130%)
  • HazardRatio < 1: Eine Beobachtung hat eine geringere Wahrscheinlichkeit in einer bestimmten Zeit ein Ereignis zu erfahren; Beispiel: HR von 0,7 bedeutet, dass die Wahrscheinlichkeit um 30% reduziert ist [(1-0,7) *100 = 30%)

Cox Regression

Siehe Tosun 2018:

  • Bei dem Verfahren steht der Effekt der Kovariaten im Mittelpunkt
  • Zeiteffekte können mit Cox-Modellen nicht getestet werden

Wichtige Annahme der proportional hazards :

\(\frac{h_i(t)}{h_j(t)} = konstant für jeden Zeitpunkt \ t\)
  • Die Konstanz der Hazard-Ratioskann sich aber ändern, wenn wir zeitvariierende unabhängige Variablen ins Modell aufnehmen; nicht aber von der Zeit
  • Wenn sie Abhängigkeit der Hazard-Ratios von der Zeit wichtig ist, ist das Cox-Modell ungeeignet
  • Ein weiteres Problem sind ties, also wenn Ereignisse zum selben Zeitpunkt stattfinden

Das kann durch den log-rank test getestet werden, der in R bei der Durchführung einer Cox-Regression mit coxph automatisch ausgegeben wird.

require(eha) # Loads 'survival' as well.
data(oldmort)
fit <- coxph(Surv(enter, exit, event) ~ sex, data = oldmort)
summary(fit)
## Call:
## coxph(formula = Surv(enter, exit, event) ~ sex, data = oldmort)
## 
##   n= 6495, number of events= 1971 
## 
##              coef exp(coef) se(coef)      z Pr(>|z|)    
## sexfemale -0.1930    0.8245   0.0456 -4.232 2.32e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##           exp(coef) exp(-coef) lower .95 upper .95
## sexfemale    0.8245      1.213     0.754    0.9016
## 
## Concordance= 0.532  (se = 0.007 )
## Rsquare= 0.003   (max possible= 0.985 )
## Likelihood ratio test= 17.73  on 1 df,   p=2.545e-05
## Wald test            = 17.91  on 1 df,   p=2.32e-05
## Score (logrank) test = 17.96  on 1 df,   p=2.254e-05

Die Annahme der proportional hazards kann auch mithilfe einer Grafik getestet werden:

plot(Surv(oldmort$enter, oldmort$exit, oldmort$event),strata = oldmort$sex, fn = "cum")

Logistische Regression der Zeitreihendaten

Empfehlung der Arbeitsschritte nach Singer and Willett (2003) (385)

Im Beispiel schauen sich Singer und Willett eine Studie an, in der Schüler danach gefragt wurden in welcher Klasse sie das erste mal Sex hatten. Die weiteren Prädiktoren sind \(\beta_1 PT\) (Ob Eltern sich geschieden haben oder nicht) und \(\beta_2 PAS\) (ob Eltern antisoziales Verhalten an den Tag legen). Die untersuchten Modelle sind also laut Singer und Willett:

Model A: \(logit \ h(t_j)=[\alpha_7 D_7+\alpha_8 D_8+...+\alpha_{12} D_{12}]\)

Model B: \(logit \ h(t_j)=[\alpha_7 D_7+\alpha_8 D_8+...+\alpha_{12} D_{12}] \ + \beta_1 PT\)

Model C: \(logit \ h(t_j)=[\alpha_7 D_7+\alpha_8 D_8+...+\alpha_{12} D_{12}] \ + \ \beta_2 PAS\)

Model D: \(logit \ h(t_j)=[\alpha_7 D_7+\alpha_8 D_8+...+\alpha_{12} D_{12}] \ + \ beta_1 PT + \beta_2 PAS\)

Bei der Durchführung dieser sequenziellen logistischen Regression sollen folgende Arbeitsschritte beachtet werden:

  1. Die Güte verschiedener alternativer Modelle testen, sodass man weiß welche Prädiktoren man nutzen sollte und welche nicht
  2. Schätzer für die Parameter berechnen und numerische Zusammenfassungen der Effektgröße berechnen
  3. angepasste Hazard und Survivalfunktionen der gewählten Werte und Prädiktoren berechnen und grafisch darstellen
  4. Konfidenzintervalle für die interessierenden Parameter erstellen

Packages für Ereignisdatenanalyse in R

library.survival
# Dateienformat muss dementsprechend geändert werden durch Funktion:
Surv() # die wichtigsten Argumente für surv() sind 'time' und 'event'

install.packages("eha")

Plots: 
par(mfrow=c(1,2)) #PLOTGRID
with(mort,plot(Surv(enter,exit,event),fn="cum"))
with(mort,plot(Surv(enter,exit,event),fn="surv"))

Literatur

Borenstein, Michael, Larry V. Hedges, Julian P.T. Higgins, and Hannah R. Rothstein. 2010. “A Basic Introduction to Fixed-Effect and Random-Effects Models for Meta-Analysis.” Research Synthesis Methods 1 (2): 97–111.

Broström, Göran. 2012. Event History Analysis with R. CRC Press.

Singer, Judith D., and John B. Willett. 2003. Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence. Oxford University Press, USA.