Mundlak-Zerlegung in R: Within- und Between-Effekt in Paneldaten trennen

Ein einfaches Random-Effects-Modell vermengt bei Paneldaten zwei analytisch zu trennende Fragestellungen: die Veränderung innerhalb einer Einheit über die Zeit einerseits, die Unterschiede zwischen den Einheiten andererseits. Die Mundlak-Korrektur trennt beide Varianzquellen, mit einer einzigen zusätzlichen Kovariate und ohne dass zeitkonstante Merkmale aus dem Modell fallen, wie es beim reinen Fixed-Effects-Ansatz der Fall ist.

Das Problem: Random-Effects-Modelle vermischen Within- und Between-Variation

Modelliert sei die Häufigkeit von Streckenstörungen im Schienenverkehr als Funktion der Instandhaltungsintensität, beobachtet über mehrere Jahre je Streckenabschnitt (strecke_id × jahr). Zu klären ist, ob eine höhere Instandhaltungsintensität die Störungen senkt. Ein naives Random-Effects-Modell liefert hier einen Koeffizienten, der Within- und Between-Variation vermischt. Verzerrt ist dieser Schätzer immer dann, wenn stärker beanspruchte Strecken zugleich intensiver gewartet werden und häufiger stören: Die Selektion schlägt sich dann als vermeintlicher Effekt nieder.

Die Mundlak-Zerlegung in einem Satz

In das Modell wird zusätzlich der Einheiten-Mittelwert jeder zeitveränderlichen Kovariate aufgenommen. Der ursprüngliche Koeffizient misst dann den reinen Within-Effekt, also die Veränderung innerhalb einer Strecke. Der Koeffizient des Mittelwerts trägt die Between-Korrektur. Der Unterschied ist derselbe wie zwischen zwei Fragen an eine Belegschaft: Verdienen Beschäftigte in Abteilung A mehr als in Abteilung B? Und: Verdient dieselbe Person nach einer Fortbildung mehr als vorher? Die erste Frage vergleicht zwischen den Einheiten, die zweite innerhalb einer Einheit. Nur die zweite sagt etwas über Wirkung. Fällt die Between-Korrektur signifikant aus (bei mehreren zeitveränderlichen Kovariaten ist der gemeinsame Wald-Test auf alle Mittelwerte das richtige Kriterium), dann unterscheiden sich Within- und Between-Steigung. Das ist der übliche Anlass, dem RE-Koeffizienten nicht als Wirkung zu trauen, und zugleich die Mundlak-Form des Hausman-Tests. Ein Nachweis von Endogenität ist es nicht: Auch echte Heterogenität kann die beiden Steigungen auseinandertreiben, etwa wenn kurz- und langfristige Reaktionen verschieden ausfallen. Ein nicht signifikanter Test ist ebenso wenig ein Freibrief; bei wenigen Einheiten fehlt häufig schlicht die Power. Zu beachten ist die Ablesung: In dieser Schreibweise ist der Koeffizient auf instand_bar die Differenz zwischen Between- und Within-Effekt, nicht der Between-Effekt selbst. Der Between-Effekt ergibt sich als Summe beider Koeffizienten. Wer ihn direkt ablesen will, zentriert die Kovariate:

library(dplyr)

# Einheiten-Mittelwert der zeitveränderlichen Kovariate
d <- d |> group_by(strecke_id) |>
       mutate(instand_bar = mean(instandhaltung)) |> ungroup()

# Mundlak-RE-Modell: Within-Effekt + Between-Korrektur
library(plm)
m <- plm(stoerungen ~ instandhaltung + instand_bar,
         data = d, index = c("strecke_id", "jahr"),
         model = "random")

library(lmtest); library(sandwich)
coeftest(m, vcov. = vcovHC(m, method = "arellano",
                             type = "HC3", cluster = "group"))
# geclustert auf strecke_id; instand_bar signifikant => Within ≠ Between

# Hybrid-/REWB-Schreibweise: beide Effekte direkt ablesbar
d <- d |> mutate(instand_dev = instandhaltung - instand_bar)
m_hyb <- plm(stoerungen ~ instand_dev + instand_bar,
             data = d, index = c("strecke_id", "jahr"),
             model = "random")
# instand_dev = Within-Effekt, instand_bar = Between-Effekt

Der Koeffizient auf instandhaltung ist damit der reine Within-Effekt; im balancierten Panel stimmt er numerisch mit dem Fixed-Effects-Schätzer überein. Das hier gerechnete Beispiel ist ein balanciertes Panel mit 60 Streckenabschnitten × 8 Jahren, 480 Beobachtungen. Dort unterscheiden sich die beiden Schätzer um 1,2e−15. Der Gewinn liegt nicht in höherer Präzision dieses Koeffizienten, sondern darin, dass zeitkonstante Merkmale wie Streckentyp oder Region im Modell erhalten bleiben. Deren Koeffizienten sind allerdings als deskriptive Between-Vergleiche zu lesen und nicht als Wirkungen: Die Mundlak-Korrektur bindet nur den Teil des Streckeneffekts, der linear auf den Kovariaten-Mittelwerten liegt. Der Rest bleibt mit ihnen konfundiert. Der Koeffizient auf instand_bar bindet jene Selektion, die den naiven Schätzer verzerrt hätte, also die zeitkonstanten Unterschiede zwischen den Strecken. Selektion, die sich über die Zeit verändert, bleibt davon unberührt.

Der Gruppenmittelwert muss über genau die Beobachtungen laufen, die anschließend auch geschätzt werden. Im Beispiel oben ist das Panel vollständig. Liegen fehlende Werte vor, gehört die Beschränkung auf vollständige Fälle vor die Mittelwertbildung; ein nachträgliches Wegrechnen der fehlenden Werte in mean() zerstört die Mundlak-Eigenschaft, weil der Mittelwert dann über mehr Zeilen gebildet wird als in die Schätzung eingehen.

Eine Annahme bleibt und wird durch die Zerlegung nicht erledigt: strikte Exogenität. Der Within-Effekt ist nur dann als Wirkung lesbar, wenn die Instandhaltung eines Jahres nicht selbst auf die Störungen früherer oder späterer Jahre reagiert und keine zeitveränderliche Größe beides zugleich treibt. Gerade bei Instandhaltung ist die Rückkopplung plausibel: Nach einem Störungsjahr wird mehr gewartet. Wo das zu erwarten ist, gehört die Zeitstruktur ins Modell (Lags, dynamische Panelschätzer) oder die Formulierung heißt „Zusammenhang innerhalb einer Strecke" statt „Wirkung".

Mundlak-Korrektur für Zähldaten: negativ-binomiales Panelmodell mit glmmTMB

Das lineare Modell bildet die Mechanik am klarsten ab. Streckenstörungen sind jedoch Zähldaten und häufig überdispers; methodisch gehört die Mundlak-Korrektur daher in ein Zähldatenmodell mit zufälligem Streckenintercept. Die Wahl der Varianzfunktion ist dabei keine Formsache: nbinom2 unterstellt Var = μ + μ²/θ, nbinom1 dagegen Var = φμ. Letzteres ist das Mixed-Model-Analogon zur Quasi-Poisson-Schätzung, die selbst keine Likelihood besitzt und in glmmTMB nicht zur Verfügung steht. Beide gewichten große und kleine Erwartungswerte unterschiedlich. Die Zerlegung funktioniert unverändert. Hinzu kommen die Verteilungsannahme und die Voraussetzung, dass der Streckeneffekt mit dem Gruppenmittelwert linear zusammenhängt. Die exakte Übereinstimmung mit dem Fixed-Effects-Schätzer gilt nur im linearen Fall:

library(glmmTMB)
mc <- glmmTMB(stoerungen ~ instandhaltung + instand_bar +
                offset(log(laenge_km)) + (1 | strecke_id),
              family = nbinom2, data = d)   # Overdispersion abgedeckt

laenge_km ist die Länge des Streckenabschnitts. Der Offset darauf gehört dazu, weil Störungen eine Zählung sind, die mit der Länge des Abschnitts mitwächst: Ohne ihn schätzt das Modell Störungen je Abschnitt statt Störungen je Kilometer. Ob er die Koeffizienten tatsächlich verschiebt, hängt davon ab, ob die Länge mit der Instandhaltungsintensität zusammenhängt. Im hier gerechneten Beispiel ist sie unabhängig gezogen. In echten Netzdaten ist das selten so. Zu beachten ist außerdem, dass ein Offset einen Koeffizienten von exakt 1 unterstellt; wer das nicht voraussetzen will, nimmt log(laenge_km) als Kovariate auf und testet gegen 1. Im gerechneten Beispiel liegt der Within-Koeffizient bei −0,135 (SE 0,052), der Koeffizient auf instand_bar bei +0,328. Der Between-Effekt ist die Summe beider Werte und liegt damit bei +0,193, auf der Log-Rate-Skala. Die Trennung ist dieselbe wie im linearen Fall: innerhalb der Strecke sinkende, zwischen den Strecken steigende Störungsraten.

Was das Modell in dieser Form noch nicht enthält, sind Zeiteffekte. Gemeinsame Schocks, ein Witterungsjahr etwa oder eine Regelwerksänderung, treffen alle Strecken zugleich und gehen sonst in den Within-Effekt ein. Jahresdummies gehören deshalb in aller Regel dazu.

Für die Fachabteilung ändert das die Aussage und nicht nur die Methode: Der Rohvergleich zwischen Strecken legt nahe, dass mehr Instandhaltung mit mehr Störungen einhergeht. Innerhalb derselben Strecke dreht sich das Vorzeichen. Unter strikter Exogenität ist dieser Within-Zusammenhang als Wirkung lesbar, und erst dann trägt die Zahl eine Budgetentscheidung.

Fazit: Within-Effekt berichten, Between-Korrektur ausweisen

Die Mundlak-Zerlegung kostet je zeitveränderlicher Kovariate eine Zeile Code und trennt die beiden Varianzquellen, an denen naive Panelmodelle scheitern: die Veränderung innerhalb einer Einheit und die Unterschiede zwischen den Einheiten. Was davon Wirkung ist, entscheidet die Identifikationsannahme und nicht die Zerlegung. Wer Effekte berichtet, ohne beide Quellen getrennt zu haben, berichtet im Zweifel das Falsche.

Drei Anlässe für ein Gespräch über Paneldaten

• Eine Kennzahl wird je Einheit über mehrere Jahre geführt (Störungen je Streckenabschnitt, Fälle je Standort, Kosten je Vertrag) und der Vergleich zwischen den Einheiten soll in eine Steuerungsentscheidung münden. Dann ist zu trennen, was Wirkung ist und was Ausgangslage.
• Eine Maßnahme wird evaluiert, in die sich die Einheiten nicht zufällig hineinbewegt haben: Wer stärker beansprucht ist, bekommt mehr Maßnahme. In einem Ex-post-Design mit Treatment- und Vergleichsgruppe schlägt diese Selektion sonst als vermeintlicher Effekt durch.
• Eine Auswertung läuft über fortlaufende Erhebungswellen und wird bei jeder Lieferung neu gerechnet. Springt ein Koeffizient zwischen den Wellen, ist zuerst zu klären, ob sich die Wirkung geändert hat oder die Zusammensetzung der Einheiten.

Die Koeffizienten +0,37 (rohe Between-Beziehung) und −0,21 (Within-Effekt) stammen aus einem tatsächlich gelaufenen Beispiel mit simulierten Paneldaten: keine Projektergebnisse und keine illustrativen Koeffizienten. Der Datenerzeugungsprozess ist so gewählt, dass stärker beanspruchte Strecken intensiver gewartet werden und häufiger stören; das Panel ist balanciert mit 60 Strecken × 8 Jahren, set.seed(42), gerechnet unter R 4.3.3 mit plm 2.6.3 und glmmTMB 1.1.8. Die Standardfehler der linearen Panelmodelle sind cluster-robust nach Arellano (HC3, geclustert auf strecke_id); der ausgewiesene SE 0,052 des Zähldatenmodells stammt aus glmmTMB und ist modellbasiert. Im Zähldatenmodell sind −0,135 der Within-Koeffizient und +0,328 der Koeffizient auf instand_bar, der in dieser Schreibweise die Differenz Between minus Within trägt; der Between-Effekt +0,193 ist die Summe der beiden. Das vollständige Skript steht unter mundlak-zerlegung.R; alle im Beitrag genannten Zahlen stammen aus diesem Lauf.

Quellen

• Mundlak, Y. (1978): On the Pooling of Time Series and Cross Section Data. Econometrica 46(1), 69–85.
• Bell, A. & Jones, K. (2015): Explaining Fixed Effects: Random Effects Modeling of Time-Series Cross-Sectional and Panel Data. Political Science Research and Methods 3(1), 133–153.
• Wooldridge, J. M. (2010): Econometric Analysis of Cross Section and Panel Data. 2. Aufl., MIT Press, Kap. 10.
• Hausman, J. A. (1978): Specification Tests in Econometrics. Econometrica 46(6), 1251–1271.
• Brooks, M. E. et al. (2017): glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling. The R Journal 9(2), 378–400.