stat-modeling-gee.qmd

```{r echo = FALSE}
pacman::p_load(tidyverse, readxl, knitr, kableExtra, Hmisc)
```

# Generalized Estimating Equations (GEE) {#sec-gee}

*Letzte Änderung am `r format(fs::file_info("stat-modeling-gee.qmd")$modification_time, '%d. %B %Y um %H:%M:%S')`*

Verallgemeinerte Schätzgleichungen (eng. *Generalized Estimating Equations*, abk. *GEE*) sind eine Methode zur Modellierung von Längsschnitt- oder Clusterdaten (eng. *longitudinal* bzw. *clustered*). Ich nutrze nur die Abkürzung GEE im weiteren Text, sonst wird mir das hier zu lang. Unter dem deutschen Begriff sind die GEE's eigentlich nicht bekannt. Jedenfalls nicht bei Anwendern. Die GEE's werden in der Regel bei nicht-normalen Daten wie binären oder Zähldaten verwendet. Damit siehst du auch schon, warum wir eigentlich nicht so oft GEE's in der Anwendung finden. Wir haben in den Agrarwissenschaften meist ein normalverteiltes Outcome $y$ und so nutzen wir dann häufig eben lineare gemischte Modelle. Damit ist das GEE auch eine Alternative für das lineare gemischte Modell. In *beiden* Modellklassen lassen sich normalverteilte, binäre und auch Zähldaten auswerten. Es ist eher eine Frage, was wir für eine Aussage über die Daten treffen wollen. Wollen wir eine beobachtungsbezogene Aussage (eng. *subject-specific*) treffen, dann nutzen wir lineare gemischte Modelle. Wollen wir eine populationsbezogene Aussage (eng. *population average*) treffen, dann nutzen wir das GEE.

::: column-margin
Das Tutorial [Getting Started with Generalized Estimating Equations](https://data.library.virginia.edu/getting-started-with-generalized-estimating-equations/) gibt nochmal einen guten Überblick über GEE's.

Ebenso liefert das Tutorial [Generalized Estimating Equations (GEE)](https://rlbarter.github.io/Practical-Statistics/2017/05/10/generalized-estimating-equations-gee/) einen schnellen Überblick über das Thema.
:::

## Annahmen an die Daten

Im folgenden Kapitel zu den Generalized Estimating Equations (GEE) gehen wir davon aus, dass die Daten in der vorliegenden Form *ideal* sind. Das heißt wir haben weder fehlende Werte vorliegen, noch haben wir mögliche Ausreißer in den Daten. Auch wollen wir keine Variablen selektieren. Wir nehmen alles was wir haben mit ins Modell. Sollte eine oder mehre Bedingungen nicht zutreffen, dann schaue dir einfach die folgenden Kapitel an.

-   Wenn du fehlende Werte in deinen Daten vorliegen hast, dann schaue bitte nochmal in das @sec-missing zu Imputation von fehlenden Werten.
-   Wenn du denkst, dass du Ausreißer oder auffälige Werte in deinen Daten hast, dann schaue doch bitte nochmal in das @sec-outlier zu Ausreißer in den Daten.
-   Wenn du denkst, dass du zu viele Variablen in deinem Modell hast, dann hilft dir das @sec-variable-selection bei der Variablenselektion.

Grundsätzlich ist das Thema GEE eher ein stiefmütterliches statistisches Thema. Ich selber habe gar nicht so viel zu GEE's gefunden, so dass wie immer gilt: *Augen auf im statistischen Straßenverkehr*! Besonders die Variablenselektion, die ja an die Modellklasse gebunden ist, mag nicht so funktionieren wie gewollt. Bitte bei GEE Fragestellungen keine automatisierte Selektion anwenden. Dann lieber über `compare_models()` aus dem R Paket `{parameters}` die Modellvergleiche direkt vergleichen.

## Genutzte R Pakete

Wir wollen folgende R Pakete in diesem Kapitel nutzen.

```{r echo = TRUE}
#| message: false
pacman::p_load(tidyverse, magrittr, broom,
               parameters, performance, geepack, gee,
               geesmv, multcomp, emmeans, scales, conflicted)
conflicts_prefer(dplyr::select)
conflicts_prefer(magrittr::set_names)
conflicts_prefer(dplyr::filter)
#cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", 
#                "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
```

An der Seite des Kapitels findest du den Link *Quellcode anzeigen*, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.

## Daten

::: callout-important
## Du musst deine Daten nach der ID sortieren

Ganz wichtig, sonst funktioniert das GEE nicht und du kriegst auch keine Warnmeldung! Du musst die Daten mit `arrange` für deine ID Spalte sortieren.
:::

Ich habe hier einmal zwei Datenbeispiel mitgebracht. Wir werden uns aber im folgenden Abschnitt dann nur die Schweine anschauen, das Kuhbeispiel können wir dann nochmal anderweitig nutzen oder aber du rechnest nochmal selber mit den Kühen. Wichtig hierbei ist, dass wir sicher sind, dass wir die Daten nach der ID Spalte der Tiere sortiert haben. Das heist, dass alle Tier ID's Zeilen wirklich beieinander stehen. Das ist wichtig, sonst schafft GEE nur eine sehr seltsame Ausgaben zu produzieren. Leider ohne eine Warnung auszugeben. Deshalb nutzen wir die Funktion `arrange()` um nochmal nach der Spalte `pig_id` zu sortieren.

```{r}
pig_gain_tbl <- read_excel("data/pig_feed_data.xlsx") |> 
  mutate(weight_gain = round(weight_gain, 2)) |> 
  arrange(pig_id)

```

In @tbl-gee-pigs sehen wir nochmal einen Auszug aus den Daten. Wir haben unsere wiederholte Messung `time`. Das heißt wir haben unsere Schweine wiederholt gemessen. Jedes Schwein für jede Behandlung fünfmal. Wir brauchen die `pig_id` um zu wissen, welche Werte der Geichtszunahme dann auch immer zu einem Ferkel gehören. Im Weiteren haben wir noch die Bucht, in der die Ferkel gehalten wurden, notiert. Die Information lassen wir aber hier erstmal im späteren Modell weg.

```{r}
#| echo: false
#| message: false
#| warning: false
#| label: tbl-gee-pigs
#| tbl-cap: Auszug aus dem Daten zu den kranken Ferkeln. Jedes Ferkel wurde wiederholt gemessen.


rbind(head(pig_gain_tbl, n = 4),
      rep("...", times = ncol(pig_gain_tbl)),
      tail(pig_gain_tbl, n = 4)) |> 
  kable(align = "c", "pipe")

```

Das zweite Datenbeispiel dient zur Veranschaulichung eines weiteres Messwiederholungsbeispiels. Wir haben drei Kühe wiederholt an drei Zeitpunkten gemessen. JEde Kuh hat immer nur die gleiche Behandlung erhalten. Das Outcome ist einmal die Anzahl an Zellen in der Milch pro ml und einmal der Fettgehalt in %. Die Daten sind in der Form relativ übersichtlich. Wir haben leider sehr wenige Messwiederholungen, so dass hier ein GEE oder aber auch ein lineares gemischtes Modell fraglich ist. Wir wollen eigentlich mindesnten fünf Level für den Clusterfaktor. Wir gehen wieder sicher, dass die Daten auch richtig nach ID sortiert sind.

```{r}
#| message: false
#| warning: false

milk_tbl <- read_csv2("data/milk_feeding.csv") |> 
  rename(cow_id = id_cow) |> 
  arrange(cow_id)

```

In @tbl-gee-milk sehen wir nochmal den Ausschnitt aus den Milchdaten. Wir haben insgesamt auch nur vierzehn Kühe gemessen, was auch nicht so viele Tiere sind. Im Ferkelbeispiel hatten wir uns 120 Ferkel angeschaut. Deshal ist dieser Datensatz sehr klein für ein komplexes Modell wie GEE.

```{r}
#| echo: false
#| message: false
#| warning: false
#| label: tbl-gee-milk
#| tbl-cap: Auszug aus Daten zu Milchkühen. Jede Kuh wurde wiederholt gemessen.

rbind(head(milk_tbl, n = 4),
      rep("...", times = ncol(milk_tbl)),
      tail(milk_tbl, n = 4)) |> 
  kable(align = "c", "pipe")
```

Gehen wir einmal auf den theoretischen Hintergrund zu GEE ein und schauen wir mal, wie wir da unser Datenbeispiel zu passt.

## Theoretischer Hintergrund

Ganz wichtig, wir gehen jetzt nicht auf den mathematischen Hintergurnd ein. Das ist auch zu *schräg*. Das will heißen, dass der mathematische Hintergrund von GEE's wirklich vieles übersteigt. Ich kann GEE's anwenden, aber ich weis nicht, wie ein GEE mathematisch funktioniert. Das muss man ja auch nicht. Deshalb hier nur die Theorie, was ein GEE macht und in welchen Hintergründen wir das GEE anwenden. Zuerst schätzt das GEE die durchschnittlichen Auswirkungen auf die Population (eng. *population average*). Betrachten wir dabei die folgenden zwei Szenarien nach @allison2009fixed:

-   *Szenario 1:* Du bist ein Arzt. Du möchtest wissen, um wie viel ein Cholesterinmedikament die Wahrscheinlichkeit eines Herzinfarkts bei deinem Patienten senkt.
-   *Szenario 2:* Du bist ein staatlicher Gesundheitsbeamter. Du möchtest wissen, wie sich die Zahl der Menschen, die an einem Herzinfarkt sterben, verändern würde, wenn alle Menschen in der Risikogruppe das Cholesterinmedikament einnehmen würden.

Im ersten Szenario wollen wir die subjektspezifischen (eng. *subject-specific*) Chancen wissen. Im zweiten Fall sind wir an der Vorhersage für die gesamte Bevölkerung interessiert. GEE kann uns Schätzungen für das zweite, aber nicht für das erste Szenario liefern. Dami sind wir schon recht weit. Wir wollen also nichts über die *einzelnen* Ferkel wissen, sondern nur über die Gesamtzahl an Ferklen mitteln. Das ist natürlich manachmal gewollt und manchmal eher nicht. In der Zucht kommt es drauf an, ob du individuelle Effekte haben möchtest, also für einen Eber oder eben die Leistung der gesamten Rasse bewerten willst. Je nachdem kannst du dan ein GEE einsetzen oder nicht. GEE's sind somit für einfaches Clustering oder wiederholte Messungen gedacht. Komplexere Designs wie verschachtelte oder gekreuzte Gruppen, z. B. verschachtelte Messwiederholungen innerhalb eines Probanden oder einer Gruppe, können nicht ohne weiteres berücksichtigt werden. Hier nutzen wir dann wieder gemischte lineare Modelle.

Ein großer Vorteil der GEE ist, dass wir eine Korrelation zwischen den wiederholten Messungen, also Ferkeln, annehmen können. Das heist, wir können die Verwandtschaft oder den zeitlichen Zusammenhang zwischend den Messwiederholungen abbilden. Dafür brauchen wir dann natürlich auch Fallzahl, die schnell mal über die hundert Beobachtungen geht. Wir können dann zwischen folgenden Zusammenhängen der Korrelation entscheiden.

-   **independence**, daher sind die Beobachtungen im Zeitverlauf sind unabhängig.
-   **exchangeable**, daher haben alle Beobachtungen im Zeitverlauf die gleiche Korrelation $\rho_{const.}$.
-   **ar1**, die Korrelation $\rho$ nimmt als Potenz der Anzahl $p$ der Zeitpunkte, die zwischen zwei Beobachtungen liegen, ab. Daher rechnen wir mit $\rho, \rho^2, \rho^3,..., \rho^p$ über die Zeitpunkte.
-   **unstructured**, daher kann die Korrelation zwischen allen Zeitpunkten unterschiedlich sein.

Leider gibt es keine automatische Auswahl. Wir müssen also überlegen, welche Korrelationmatrix am besten passen würde. Da *unstructured* sehr viel Fallzahl benötigt um valide zu sein, nehme ich meistens *exchangeable*, wenn ich ein GEE rechne. Eine unabhänige Korrealtion anzunehmen macht wenig Sinn, dann brauche ich auch kein GEE rechnen. Die Korrelation ist ja die Stärke von einem GEE.

Wir haben nun zwei R Pakete, die beide das gleiche tun, nämlich ein GEE rechnen. Wir haben die Wahl zwischen dem R Paket `gee` und der Funktion `gee()` sowie der Funktion `geeglm()` aus dem R Paket `{geepack}`. Ich neige zu dem letzteren Paket. Das R Paket `{geepack}` ist etwas neueren Ursprungs und funktioniert bisher recht reibungslos.

## Modellieren mit `gee()`

Leider ist es so, dass wir kaum kontrollieren können, was alles aus den Funktionen in die R Console geschrieben wird. Die Funktionen sind schon recht alt und es gab mal einen Trend, dass eine Funktion immer schön was wiedergeben soll. Das ist natürlich immer etwas nervig, wenn man das nicht will. Wir erhalten also bei der Funktion `gee()` immer die Koeffizienten des Modells ausgegeben, ob wir es nun in einem Objekt speichern oder auch nicht. Ich finde sowas immer ziemlich nervig.

Also wir bauen wir uns unser `gee` Modell? Zuerst kommt wie immer die `formula`, da ändert sich nichts. Wir nehmen in unser Modell als Outcome die Gewichtszunahme und als Einflusvariablen dann die Behandlung sowie die Zeit und die Interaktion zwischen der Behandlung und der Zeit. Die Daten sind auch gleich. Erst mit der Option `id =` ändert sich was. Hier geben wir die Spalte ein, in der die ID's der Ferkel bzw. der Beobachtungen stehen. Das war es auch schon für den CLustereffekt. Dann nach die Verteilungsfamilie, wir können hier auch für nicht normalverteilte Daten ein GEE schätzen. Zum Abschluss noch die Korrelationsstruktur definiert. Wir nehmen hier *exchangeable*, diese Korrelationsstruktur ist für den Anfang immer ganz gut und macht auch häufig Sinn.

```{r}
#| message: false
#| warning: false

gee_fit <- gee(weight_gain ~ treatment + treatment * time,
               data = pig_gain_tbl, 
               id = pig_id, 
               family = gaussian,
               corstr = "exchangeable")
```

Nachdem wir das Modell gefittet haben, können wir uns einmal die Korrelationsstruktur anschauen. Da ist die Funktion `gee` wirklich gut. Die Korrelationsstuktur können wir uns einfach so rausziehen.

```{r}
#| message: false
#| warning: false

pluck(gee_fit, "working.correlation") |> 
  round(3)
```

Was sehen wir? Natürlich muss auf der Diagonalen eine 1 stehen, den untereinander sind die Variablen ja identisch und damit mit 1 korreliert. Auf der Nicht-Diagonalen finden wir dann die Korrelation untereinander. Da wir *exchangeable* für die Korrelationsstruktur gewählt haben, haben wir überall die gleiche Korrelation. Alle Ferkel sind untereinander über die Zeitpunkte gleich mit $\rho = 0.82$ korreliert.

Wir lassen uns jetzt noch die Modellparameter ausgeben und schauen uns einmal an, ob wir was signifikantes gefunden haben.

```{r}
#| message: false
#| warning: false

gee_fit |> model_parameters()
```

Wir sehen, dass es einen signifikanten Unterschied in der Zeit gibt, das war ja auch zu erwarten, denn mit der Zeit werden die Ferkel schwerer. Wir haben aber auch einen schwach signifikanten Effekt zwischen `feed_10` und `feed_20` mit einem $p$-Wert von $0.041$. Hier machen wir kurz Stop, dann geht es aber in dem Abschnitt zu den Posthoc Tests mit dem Modell weiter. Wir wollen ja noch für alle Behanlungslevel einen paarweisen Vergleich rechnen.

## Modellieren mit `geeglm()`

Der eigentlcihe Unetrschied zwischen der Funktion `gee()` und `geeglm()` ist, dass sich im Hintergrund eine Menge anders abspielt, das wir nicht sehen. Für mich war `geeglm()` immer schneller und stabiler. Der einzige Grund war immer nochmal ein `gee()` laufen zu lassen, da sich die Korrelationsmatrix so einfach aus dem `gee()` Objekt ziehen lassen lässt.

```{r}
#| message: false
#| warning: false

geeglm_fit <- geeglm(weight_gain ~ treatment + treatment * time,
                     data = pig_gain_tbl, 
                     id = pig_id, 
                     family = gaussian,
                     corstr = "exchangeable")

```

Wir erhalten zwar auch die geschätzte Korrelation, aber nicht in so einer schönen Matrix. Also ist es dann Geschmackssache. Du weist dann ja, das wir mit *exchangeable* überall die gleiche Korrelation angenommen haben.

```{r}
#| message: false
#| warning: false

pluck(geeglm_fit, "geese", "alpha")
```

Am Ende schauen wir uns dann nochmal den Fit aus dem `geeglm()` Modell an. Und stellen fest, dass das Modell numerisch fast identisch ist. Wir haben also nur dir Wahl in der Darstellungsform und in der Geschwindigkeit.

```{r}
#| message: false
#| warning: false

geeglm_fit |> model_parameters()
```

Wir haben also dann zwei Funktionen, die wir nutzen können. Am Ende kannst du dann beide ausprobieren. Machmal hat man mit `geeglm()` etwas mehr Glück, wenn die Daten einfach mal nicht wollen.

## Multipler Vergleich mit `emmeans`

Gut, soweit sind wir dann gekommen. Wir haben unser Modell gefittet und meistens wollen wir dann noch einen all-pair Vergleich bzw. den paarweisen Vergleich rechnen. Das machen wir erst einmal mit der Funktionalität aus dem R Paket `{emmeans}`, das uns erlaubt auch das *compact letter display* wiederzugeben. Wenn dich mehr zum Prozess des Codes für die Nutzung von `{emmeans}` interessiert, dann schaue doch einfach nochmal ins @sec-posthoc. In dem Kapitel zu den multiplen Vergleichen erkläre ich dir nochmal genauer den Funktionsablauf.

Wichtig ist, dass wir unsere Vergleiche mit Bonferroni adjustieren. Wenn du das nicht möchtest, dann musst du `adjust = "none"` auswählen. Sonst machen wir die Ausageb nochmal `tidy()` und dann runden wir noch. Wir erhalten dann das *compact letter display* wieder.

```{r}
#| message: false
#| warning: false

res_gee <- geeglm_fit |> 
  emmeans(~ treatment) 

res_gee_cld <- res_gee |> 
  cld(adjust = "bonferroni", Letters = letters) |> 
  tidy() |> 
  select(treatment, estimate, conf.low, conf.high, .group) |> 
  mutate(across(where(is.numeric), round, 2))

res_gee_cld 
```

Wenn dich die Abbildungen und weiteres interessieren, dann schaue einfach nochmal ins Kapitel zu den multiplen vergleichen. Dort zeige ich dann wie wir das *compact letter display* in eine Abbildung ergänzen. Der Ablauf ist auch im @sec-mixed zu den linearen gemischten Modellen gezeigt.

Wir sehen, dass sich die Gruppe `feed_10` von der Gruppe `feed_20` unterscheidet. Beide Gruppen haben nicht den gleichen Buchstaben. Die Gruppe `feed_10+10` unterscheidet sich weder von der Gruppe `feed_10` noch von der Gruppe `feed_20`. Wir können uns im folgenden Codeblock dann auch die $p$-Werte für die Vergleiche wiedergeben lassen. Die Aussagen sind die selben.

```{r}
#| message: false
#| warning: false

res_gee_tbl <- res_gee |> 
  contrast(method = "pairwise", adjust = "bonferroni") |> 
  tidy(conf.int = TRUE) |> 
  mutate(p.value = pvalue(adj.p.value),
         across(where(is.numeric), round, 2)) |> 
  select(contrast, estimate, 
         conf.low, conf.high, p.value) 

res_gee_tbl
```

Die Funktion `emmeans()` hätten wir auch mit dem Modell aus dem `gee()` Fit nutzen können. Als letzten Abschnitt wollen wir uns jetzt noch eine Besonderheit der GEE Varianzschätzung anschauen.

## Multipler Vergleich mit `multcomp` und `geesmv`

Natürlich haben wir uns nicht in den Details wie ein GEE funktioniert verloren. Es ist nun aber so, dass ein GEE auf sehr unterschiedliche Art und Weise die Korrelationsstruktur und die Varianzen dahinter schätzen kann. Je nach geschätzter Varianz kommen natürlich auch eventuell ganz andere Signifikanzen aus dem Modell. Deshalb hat sich wirklich eine Horde an *mathematischen* Statistikern an der Varianzschätzung im GEE abgearbeitet.

::: column-margin
[geesmv: Modified Variance Estimators for Generalized Estimating Equations](https://rdrr.io/cran/geesmv/man/geesmv.html)
:::

Das R Paket `{geesmv}` bietet ganze neun Implementierungen von Schätzern für die Varianz/Covarianzstruktur der Daten an. Jetzt stellt sich die Frage, welche Implementierung für den Varianzschätzer denn nun nutzen? Zum einen hat natürlich die geschätzte Varianz einen nicht zu unterschätzenden Effekt auf die Signifikanz der Koeffizienten des GEE Models. Zum anderen ist aber der multiple Vergleich nach dem Schätzen des Modells und dem *getrennten* Schätzen der Varianz sehr mühselig. Leider helfen uns auch unsere Standardpakete nicht so richtig weiter. Die Funktionalität ist nicht für `{geesmv}` implementiert. Was wiederum dafür spricht, dass der Bedarf von Anwendern sehr eingeschränkt zu seien scheint. Nun müssen wir folgende epischen Schritte durchführen um einen multiplen Vergleich rechnen zu können.

1)  Wir fitten unser `geeglm()` Modell in der *mean parametrization*, dass heist wir entfernen den Intercept aus dem Modell und lassen unser Modell somit durch den Urspung laufen. Im Prinzip setzen wir den Intercept auf 0 und erhalten so die Mittelwerte jedes Levels des Faktors `treatment`.
2)  Wir speichern die $\beta$-Koeffizienten von dem `treatment` aus unserem GEE Modell in einem Objekt ab.
3)  Wir rechnen mit der gleichen Angabe wie vorher das `geeglm()` Modell eine der neun Funktion. Ich habe hier zufällig die Funktion `GEE.var.lz()` gewählt. Wir speichern die Ausgabe der Varianz der Koeffizienten in einem Objekt.
4)  Wir kombinieren die $\beta$-Koeffizienten und die Varianz in einem Objekt mit der Funktion `left_join()`.
5)  Wir bauen uns unsere eigene Kontrastmatrix in der steht welches Level der Behandlung mit welchen anderen Level verglichen werden soll.
6)  Wir übergeben alle Einzelteile an die Funktion `glht()` aus dem R Paket `{multcomp}` und rechnen unseren multiplen Vergleich.

Na dann mal auf. Gehen wir die Schritte einmal nacheinander durch und schauen, was wir da so alles gemacht haben. Nochmal Achtung, hier musst du wirklich schauen, ob sich der Aufwand lohnt. Ich zeige es hier einmal, den in bestimmten Fällen kann sich eine andere Implementierung für die Schätzung der Varianz durchaus lohnen. Denn aus Erfahrung weiß ich, dass der Standardvarianzschätzer nicht immer der beste Schätzer sein muss [@kruppa2021comparison].

Im Folgenden schätzen wir einmal ein ganz normales GEE Modell mit der Funktion `geeglm()`. Wir werden aber nur die Koeffizienten brauchen. Die Varianz der Koeffizienten nutzen wir nicht. Ebenso brauchen wir die *mean* Parametrisierung, dass heißt wir setzen den Intercept auf 0.

```{r}
#| message: false
#| warning: false

geeglm_fit <- geeglm(weight_gain ~ 0 + treatment + treatment * time,
                     data = pig_gain_tbl, 
                     id = pig_id, 
                     family = gaussian,
                     corstr = "exchangeable")
```

Wir speichern einmal die Koeffizienten in dem Objekt `beta_tbl`. Die brauchen wir später um die paarweisen Vergleiche zu rechnen.

```{r}
#| message: false
#| warning: false

beta_tbl <- coef(geeglm_fit) |> 
  enframe()
```

Um die Varianz der Koeffizienten zu schätzen nutzen wir jetzt *eine* der Implementierungen in `geesmv`. Ich habe mich etwas zufällig für die Implementierung `GEE.var.lz()` entschieden. Diese Funktion liefer *nur* die Varianz der Koeffizienten. Leider aber nicht auch gleich noch die Koeffizienten dazu... deshalb der blöde doppelte Schritt. Wir speichern dann die Varianzen in dem Objekt `vbeta_tbl`.

```{r}
#| message: false
#| warning: false

gee_lz_vcov <- GEE.var.lz(weight_gain ~ 0 + treatment + treatment * time,
                          data = as.data.frame(pig_gain_tbl), 
                          id = "pig_id",
                          family = gaussian,
                          corstr = "independence") 

vbeta_tbl <- pluck(gee_lz_vcov, "cov.beta") |> 
  enframe()
```

Jetzt verbinden wir noch die beiden Objekte `beta_tbl` und `vbeta_tbl` über die Funktion `left_join()`. Wir können mit der Funktion zwei Datensätze nach einer gemeinsamen Spalte zusammenführen. Dann müssen zwar die Einträge in der Spalte gleich sein, aber die Sortierung kann anders sein. Dann müssen wir noch die Zeilen rausfiltern in denen die Behandlungsmittelwerte sind. Am Ende benennen wir die Spalten noch sauber nach dem was die Spalten sind.

```{r}
#| message: false
#| warning: false

coef_tbl <- left_join(beta_tbl, vbeta_tbl, by = "name") |> 
  filter(str_detect(name, "time", negate = TRUE)) |> 
  set_names(c("parameter", "beta", "vbeta"))

```

Das war jetzt ein Angang. Leider geht es nicht so einfach weiter. Wir müssen uns für die Vergleiche die Kontrastmatrix selberbauen. Wir machen einen paarweisen Vergleich, also wählen wir den Tukey Kontrast aus.

```{r}
#| message: false
#| warning: false

contrMat_n <- setNames(rep(1, length(coef_tbl$parameter)),
                       coef_tbl$parameter) |> 
  contrMat(type = "Tukey")

contrMat_n 
```

Nun können wir alles zusammenbringen. Wir nutzen die Helferfunktion `parm()` aus dem R Paket `{multcomp}` um diei Koeffizienten richtig in `glht()` zuzuordnen. Dann noch der Kontrast mit rein in die Funktion und wir können unseren Vergleich rechnen. Leider fehlen noch die Freiheitsgrade, die wären dann in unserem Fall null, das ist aber Unsinn. Wir ergänzen die Freiheitsgrade aus unserem ursprünglichen Modell für die Koeffizienten.

```{r}
#| message: false
#| warning: false

mult_gee <- glht(parm(coef = coef_tbl$beta, 
                      vcov = diag(coef_tbl$vbeta)), 
                 linfct = contrMat_n)
mult_gee$df <- geeglm_fit$df.residual
```

Jetzt können wir uns die $p$-Werte und die 95% Konfidenzintervalle wiedergeben lassen. Du musst echt überlegen, ob sich der Aufwand lohnt. Wir erhalten hier jetzt kein signifikanten Unterschied mehr. Das liegt daran, dass wir in diesem Fall höhere Varianzen geschätzt haben als das `geeglm()` normalerweise tun würde. Höhere Varianzen der Koeffizienten, weniger signifikante Koeffizienten. Und dann auch weniger signifikante paarweise Unterschiede.

```{r}
#| message: false
#| warning: false

mult_gee |> 
  tidy(conf.int = TRUE) |> 
  select(contrast, estimate, conf.low, conf.high, adj.p.value)
```

Als Fazit nehmen wir mit, dass wir noch die Möglichkeit haben auf andere Art und Weise die Varianz in einem GEE zu schätzen. Ob uns das hilft steht auf einen anderem Blatt, aber wir haben die Möglichkeit hier noch nachzuadjustieren, wenn es mit dem Varianzschätzer klemmen sollte. Großartig unterstützt wird das Paket nicht, dass sieht man ja schon daran wie Oldschool die Analyse ist.

## Referenzen {.unnumbered}