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")54 Generalized Estimating Equations (GEE)
Letzte Änderung am 02. April 2024 um 09:52:43
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.
Das Tutorial Getting Started with Generalized Estimating Equations gibt nochmal einen guten Überblick über GEE’s.
Ebenso liefert das Tutorial Generalized Estimating Equations (GEE) einen schnellen Überblick über das Thema.
54.1 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 Kapitel 43 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 Kapitel 42 zu Ausreißer in den Daten.
- Wenn du denkst, dass du zu viele Variablen in deinem Modell hast, dann hilft dir das Kapitel 42 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.
54.2 Genutzte R Pakete
Wir wollen folgende R Pakete in diesem Kapitel nutzen.
An der Seite des Kapitels findest du den Link Quellcode anzeigen, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.
54.3 Daten
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.
pig_gain_tbl <- read_excel("data/pig_feed_data.xlsx") |>
mutate(weight_gain = round(weight_gain, 2)) |>
arrange(pig_id)In Tabelle 54.1 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.
| time | pig_id | cove | treatment | weight_gain |
|---|---|---|---|---|
| 1 | 1 | 1 | feed_10 | 44.29 |
| 2 | 1 | 1 | feed_10 | 59 |
| 3 | 1 | 1 | feed_10 | 49.96 |
| 4 | 1 | 1 | feed_10 | 66.17 |
| … | … | … | … | … |
| 2 | 120 | 10 | feed_20 | 56.17 |
| 3 | 120 | 10 | feed_20 | 58.13 |
| 4 | 120 | 10 | feed_20 | 68.48 |
| 5 | 120 | 10 | feed_20 | 34.99 |
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.
milk_tbl <- read_csv2("data/milk_feeding.csv") |>
rename(cow_id = id_cow) |>
arrange(cow_id)In Tabelle 54.2 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.
| cow_id | trt | time_point | cell_count | fat_perc |
|---|---|---|---|---|
| 1 | 1 | 1 | 1932 | 0.69 |
| 1 | 1 | 2 | 6771 | 0.94 |
| 1 | 1 | 3 | 2225 | 0.01 |
| 2 | 0 | 1 | 2572 | 0.15 |
| … | … | … | … | … |
| 13 | 1 | 3 | 8445 | 0.06 |
| 14 | 1 | 1 | 19707 | 0.95 |
| 14 | 1 | 2 | 9428 | 0.5 |
| 14 | 1 | 3 | 8184 | 0.96 |
Gehen wir einmal auf den theoretischen Hintergrund zu GEE ein und schauen wir mal, wie wir da unser Datenbeispiel zu passt.
54.4 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 Allison (2009):
- 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.
54.5 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.
gee_fit <- gee(weight_gain ~ treatment + treatment * time,
data = pig_gain_tbl,
id = pig_id,
family = gaussian,
corstr = "exchangeable") (Intercept) treatmentfeed_10+10 treatmentfeed_20
56.441300 1.573500 3.378825
time treatmentfeed_10+10:time treatmentfeed_20:time
-2.358900 0.207000 0.182875 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.
pluck(gee_fit, "working.correlation") |>
round(3) [,1] [,2] [,3] [,4] [,5]
[1,] 1.000 0.082 0.082 0.082 0.082
[2,] 0.082 1.000 0.082 0.082 0.082
[3,] 0.082 0.082 1.000 0.082 0.082
[4,] 0.082 0.082 0.082 1.000 0.082
[5,] 0.082 0.082 0.082 0.082 1.000Was 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.
gee_fit |> model_parameters()Parameter | Coefficient | SE | 95% CI | z | p
------------------------------------------------------------------------------------
(Intercept) | 56.44 | 1.23 | [52.53, 60.35] | 28.27 | < .001
treatment [feed_10+10] | 1.57 | 1.65 | [-3.96, 7.11] | 0.56 | 0.341
treatment [feed_20] | 3.38 | 1.66 | [-2.16, 8.91] | 1.20 | 0.041
time | -2.36 | 0.19 | [-3.49, -1.22] | -4.07 | < .001
treatment [feed_10+10] * time | 0.21 | 0.24 | [-1.40, 1.81] | 0.25 | 0.386
treatment [feed_20] * time | 0.18 | 0.25 | [-1.42, 1.79] | 0.22 | 0.456 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.
54.6 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.
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.
pluck(geeglm_fit, "geese", "alpha") alpha
0.082597 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.
geeglm_fit |> model_parameters()Parameter | Coefficient | SE | 95% CI | Chi2(1) | p
--------------------------------------------------------------------------------------
(Intercept) | 56.44 | 1.23 | [54.04, 58.84] | 2120.76 | < .001
treatment [feed_10+10] | 1.57 | 1.65 | [-1.66, 4.81] | 0.91 | 0.341
treatment [feed_20] | 3.38 | 1.66 | [ 0.13, 6.62] | 4.17 | 0.041
time | -2.36 | 0.19 | [-2.73, -1.99] | 153.65 | < .001
treatment [feed_10+10] * time | 0.21 | 0.24 | [-0.26, 0.67] | 0.75 | 0.386
treatment [feed_20] * time | 0.18 | 0.25 | [-0.30, 0.66] | 0.56 | 0.456 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.
54.7 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 Kapitel 34. 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.
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 # A tibble: 3 x 5
treatment estimate conf.low conf.high .group
<chr> <dbl> <dbl> <dbl> <chr>
1 feed_10 49.4 47.0 51.7 " a "
2 feed_10+10 51.6 49.1 54.0 " ab"
3 feed_20 53.3 51.0 55.6 " b" 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 Kapitel 52 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.
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# A tibble: 3 x 5
contrast estimate conf.low conf.high p.value
<chr> <dbl> <dbl> <dbl> <chr>
1 feed_10 - (feed_10+10) -2.19 -5.56 1.17 0.355
2 feed_10 - feed_20 -3.93 -7.2 -0.66 0.012
3 (feed_10+10) - feed_20 -1.73 -5.06 1.59 0.636 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.
54.8 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.
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.
- 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 Faktorstreatment. - Wir speichern die \(\beta\)-Koeffizienten von dem
treatmentaus unserem GEE Modell in einem Objekt ab. - Wir rechnen mit der gleichen Angabe wie vorher das
geeglm()Modell eine der neun Funktion. Ich habe hier zufällig die FunktionGEE.var.lz()gewählt. Wir speichern die Ausgabe der Varianz der Koeffizienten in einem Objekt. - Wir kombinieren die \(\beta\)-Koeffizienten und die Varianz in einem Objekt mit der Funktion
left_join(). - Wir bauen uns unsere eigene Kontrastmatrix in der steht welches Level der Behandlung mit welchen anderen Level verglichen werden soll.
- 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 (Kruppa und Hothorn 2021).
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.
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.
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.
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") treatmentfeed_10 treatmentfeed_10+10 treatmentfeed_20
56.441300 58.014800 59.820125
time treatmentfeed_10+10:time treatmentfeed_20:time
-2.358900 0.207000 0.182875 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.
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.
contrMat_n <- setNames(rep(1, length(coef_tbl$parameter)),
coef_tbl$parameter) |>
contrMat(type = "Tukey")
contrMat_n
Multiple Comparisons of Means: Tukey Contrasts
treatmentfeed_10 treatmentfeed_10+10
treatmentfeed_10+10 - treatmentfeed_10 -1 1
treatmentfeed_20 - treatmentfeed_10 -1 0
treatmentfeed_20 - treatmentfeed_10+10 0 -1
treatmentfeed_20
treatmentfeed_10+10 - treatmentfeed_10 0
treatmentfeed_20 - treatmentfeed_10 1
treatmentfeed_20 - treatmentfeed_10+10 1Nun 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.
mult_gee <- glht(parm(coef = coef_tbl$beta,
vcov = diag(coef_tbl$vbeta)),
linfct = contrMat_n)
mult_gee$df <- geeglm_fit$df.residualJetzt 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.
mult_gee |>
tidy(conf.int = TRUE) |>
select(contrast, estimate, conf.low, conf.high, adj.p.value)# A tibble: 3 x 5
contrast estimate conf.low conf.high adj.p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 treatmentfeed_10+10 - treatmentfeed_10 1.57 -2.30 5.45 0.607
2 treatmentfeed_20 - treatmentfeed_10 3.38 -0.510 7.27 0.103
3 treatmentfeed_20 - treatmentfeed_10+10 1.81 -1.88 5.49 0.483Als 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
Allison PD. 2009. Fixed effects regression models. SAGE publications.
Kruppa J, Hothorn L. 2021. A comparison study on modeling of clustered and overdispersed count data for multiple comparisons. Journal of Applied Statistics 48: 3220–3232.