Korrelierte Daten simulieren für PROC GEE

[sauer]

Hallo und guten Abend,

mit PROC GEE möchte ich verschiedene Szenarien durchspielen und brauche dazu simulierte Datensätze korrelierter Daten.

Es sollen dichotome Datenpaare sein, für die ich sowohl die Verteilung als auch die Kovarianzstruktur vorgeben möchte. Die Cluster sollen nicht alle gleich groß sein. Die Wertepaare mit x not eq y sollen selten sein.

Hat denn jemand eine Idee, wie man so etwas mit SAS hinbekommen kann?

Vor allen scheitere ich an der Kovarianz 😞

Danke und viele Grüße

Frank

[FreelanceReinh]

Hallo Frank,

ich denke, es wäre hilfreich, wenn Du die Frage noch weiter konkretisieren würdest.

Vereinfachtes Beispiel: "Es sollen drei Bernoulli-verteilte Zufallsvariablen X1, X2, X3 mit vorgegebenen Parametern p1, p2, p3 simuliert werden, für deren Kovarianzstruktur gelten möge: Corr(X1,X2)=Corr(X1,X3)=Corr(X2,X3)=a mit vorgegebenem a."

Viele Grüße

Reinhard

[sauer]

Hallo Reinhard, vielen Dank für Deinen Hinweis.

Allerdings bin ich mit den Begrifflichkeiten nicht so sattelfest, daher versuche ich es "andersherum".

Mein code ist:

proc genmod data = test;

class pat_id;

model y= / dist=Poisson link=log;

repeated subject=pat_id / type=ar(1);

where x=1;

run;

Ich will nun für x,y (welche beide dichotom sein sollen) Datensätze simulieren, wo ich mit DIST, LINK und TYPE variieren kann.

Dann will ich die Interzept-Schätz- und die QIC Werte untersuchen.

Z.B. was passiert mit dem Beta und QIC, wenn ich bei Poissonverteilten Daten bei DIST Normal vorgebe usw.

Bei meinen real-Daten bestimmt nur TYPE, was für Beta herauskommt ...

Hoffentlich konnte ich mich verständlich machen 😉

Danke und Grüße

Frank

[FreelanceReinh]

Hallo Frank,

danke für die ergänzenden Informationen.

Wenn es sich um eine wissenschaftliche Auswertung von Patientendaten handelt ("pat_id"), sollte eigentlich ein statistischer Analyseplan vorliegen, aus dem insbesondere das für die zu erhebenden Daten angenommene mathematisch-statistische Modell hervorgeht. Danach richtet sich dann der zu verwendende SAS-Code.

Hier scheint es aber eher um die Simulation geeigneter Daten zu vorgegebenem SAS-Code zu gehen.

Sowohl im Programmcode als auch im Text ist von der Poisson-Verteilung die Rede. Welche der beteiligten Variablen wird denn speziell als (näherungsweise) Poisson-verteilt angenommen, wenn doch x und y ausdrücklich dichotom sein sollen?

Ist x eine Prädiktorvariable oder werden darüber bloß Datensätze aus dem Analyse-Dataset ausgewählt, wie es im Beispielprogramm der Fall ist?

Wie viele Messwiederholungen sind denn maximal zu erwarten?

[sauer]

Hallo Reinhard, vielen Dank für Deinen Kommentar.

Jetzt schreibe ich doch etwas mehr über den Hintergrund.

Ich habe 1700 Probanden, von denen jeweils ein 20-min-EKG genommen wurde. In der Summe sind es ca. 2,3 Millionen Herzschläge. Ein Algorithmus hat jeden Schlag identifiziert als normal, supraventrikulär, ventrikulär oder Artefakt.

Ein Arzt hat danach jeden Schlagtyp überprüft und bestätigt oder geändert (das ist mein Gold Standard, so dass ich für jeden Schlag und Typ ein Pärchen 0,0 oder 0,1 oder  1,0 oder 1,1 habe).

Die Anzahl der Schläge pro Proband schwankt zwischen 1000 und 2000.

Das ist alles lange vor meiner Zeit passiert ;-).

Ich will die Performance des Algorithmus beurteilen.

Bisher habe ich mit PROC FREQ die globale Sensitivität / Spezifität je Schlagtyp und über die Einzelwerte jedes EKG die gemittelte (jedes EKG gleich gewichtet).

Meine Idee war nun, mittels GEE die Korrelation der Schläge innerhalb der einzelnen EKGs (=pat_id) zu modellieren (mit AR).

Poisson habe ich gewählt, weil es sehr seltene Ereignisse sind (ca. 12 Tsd. bei 2,3 Mio.).

Verwunderlich ist nun, egal was ich als DIST oder LINK wähle, es bestimmt immer TYPE, was als Ergebnis erscheint, auf die Kommastelle genau.

Und die Unterschiede sind beträchtlich, obwohl meiner Meinung nach selbst eine falsche Wahl der Kovarianzstruktur TYPE das Ergebnis nicht so dramatisch verändern dürfte.

Auch die Betrachtung von QIC/QICu bringt kein Licht ins Dunkel.

In Ermangelung anderer (geclusterter) Datensätze wollte ich mit simulierten Daten experimentieren, um vielleicht ein Gefühl dafür zu bekommen, was hier eigentlich passiert.

Hoffentlich wird mein Ansinnen so klarer 🙂

Viele Grüße

Frank

[FreelanceReinh]

Hallo Frank,

besten Dank für die nähere Erläuterung des Hintergrundes.

Die Variablen X und Y beschreiben also nach meinem Verständnis jeweils für einen Schlagtyp (z. B. "ventrikulär"), ob er bei einem gegebenen Herzschlag nach Einschätzung des Arztes bzw. des Algorithmus vorliegt (1) oder nicht (0).

Die nicht-normalen Schlagtypen sind vermutlich selten (wenige Einsen). An die Poisson-Verteilung, die häufig als "Verteilung seltener Ereignisse" bezeichnet wird, könnte man denken, wenn etwa die absolute Häufigkeit eines der seltenen Schlagtypen für einen Probanden modelliert werden sollte. Für die Modellierung einer dichotomen Variable wie Y taugt sie jedoch nicht, da ja keine Werte >1 vorkommen.

Nur die Binomialverteilung mit n=1 (=Bernoulli-Verteilung) kommt dafür sinnvoll in Frage. Wenn man hier also überhaupt mit GEE modellieren will, sollte DIST=BIN verwendet werden. Das von Dir anvisierte Intercept-only-Modell würde dann, soweit ich sehe, den mit der Linkfunktion transformierten Erwartungswert von Y schätzen. Die Wahl der Linkfunktion dürfte in diesem Fall belanglos sein.

Mir fiel zu dem Thema der Kappa-Koeffizient ein (Cohen's Kappa, "inter-rater agreement"). Dieser kann leicht mit der AGREE-Option des TABLES-Statements von PROC FREQ berechnet werden. Das wäre vielleicht eine Alternative. Mit dieser Statistik habe ich allerdings noch nie gearbeitet und weiß auch nicht, inwieweit man hier die Korrelationen in die Betrachtung einbeziehen könnte.

Zu der ursprünglichen Frage nach Simulationen: Dafür wäre mein Ansatz, von der gemeinsamen Verteilung der korrelierten dichotomen Zufallsvariablen auszugehen und für die 2^k Parameter dieser Verteilung (wenn k die Anzahl der Variablen ist) Gleichungen anhand der gewünschten Vorgaben (also Randverteilungen und Korrelationen) aufzustellen. Jede Lösung des sich so ergebenden (zum Glück letztlich) linearen Gleichungssystems liefert dann die Parameter der Bernoulli-Verteilungen, mit denen man die RAND-Funktion von SAS füttern muss.

Für kleine Werte von k ist das gut machbar. Wenn k dagegen, wie es sich hier andeutet, "zwischen 1000 und 2000" liegt, muss man sich freilich mehr Gedanken über die praktische Lösung des Gleichungssystems machen (vielleicht SAS/IML einsetzen?). Zum bloßen "Experimentieren" genügen aber womöglich kleine k.

[sauer]

Danke für Deine Antwort, Reinhard.

Mit Cohen's Kappa war ich auch schon unterwegs, aber da hat mein Professor nur abgewinkt und gemeint, das wäre ein Thema für sich ;-).

Natürlich habe ich auch DIST=B gerechnet, aber wie schon gesagt, das Ergebnis (also der Intercept-Schätzer) ist bei B, N, P immer das gleiche und zwar kommastellengenau. Nur TYPE wirkt sich aus ..... um mal Zahlen zu nennen für die Sensitivität  supraventrikulär:

IND = 38 (wie PROC FREQ)

AR = 42

CS = 60

Und gerade TYPE sollte nach meinem Verständnis lediglich die Effizienz der Schätzung beeinflussen.

Die QIC Werte sagen: DIST=N, LINK=Identity und TYPE=IND bietet die beste Anpassung, es folgt (N, Identity, AR).

Irgendwie muss ich das Ganze interpretieren und möglichst auch erklären 😞

Herzlich!

Frank

[FreelanceReinh]

Hallo Frank,

eine einzelne Zufallsvariable Y, die nur die Werte 0 und 1 annehmen kann, hat nun einmal zwangsläufig eine Bernoulli-Verteilung und deren Parameter p ist zugleich der Erwartungswert E(Y). Wenn nun wie hier der Funktionswert g(E(Y)) einer Link-Funktion g mit Hilfe eines Modells geschätzt wird, ist es absehbar, dass selbst eine fehlerhaft spezifizierte Verteilung von Y (etwa "DIST=N") nichts am Resultat ändert.

Dagegen sollte die Wahl der Link-Funktion (z. B. log oder logit) sehr wohl einen Einfluss auf diese Schätzung haben, denn log(p), logit(p) usw. sind ja im allgemeinen verschieden.

Wie sich unterschiedliche Kovarianzstrukturen, insbesondere Korrelationen, in diesem speziellen statistischen Modell auf geschätzte Kenngrößen wie die Sensitivität auswirken, kann man vielleicht am besten nachvollziehen, indem man ein kleines Beispiel (etwas mühsam) "per Hand" durchrechnet. Die Tatsache, dass anscheinend keine Prädiktorvariablen vorkommen, macht das leichter. Eine dazu passende Simulation, wie in meiner vorigen Antwort skizziert, kann man zu Kontrollzwecken parallel durchführen. Wenn dabei Schwierigkeiten auftreten, helfe ich gern weiter.

Viel Erfolg!

[sauer]

Hallo Reinhard,

nun habe ich in den vergangenen Tagen mit verschiedenen Varianten herumprobiert und (dank Deines Hinweises) zumindest eine Idee davon entwickelten können, warum DIST und LINK (ich transformiere natürlich zurück) keinen Einfluss haben.

Die Kovarianz ist bei meinen Simulationen immer Null, ich weiß also nicht, wie ich die hinein"simulieren" kann 😞

Aber das ist jetzt auch nicht mehr wichtig, denn die großen Unterschiede der Schätzer bei Wechsel von TYPE sind wohl durch Ausreißer zu erklären. Wenn ich von meinen 1700 Clustern die 10 mit den meisten "FalschNegativen" herausnehme, scheint wieder die Sonne 🙂

Ein Kollege hatte mir gesagt, weil ich keine Kovariablen habe, kann die Prozedur nur am Interzept "anpassen". Und wenn sich das dann mit den Ausreißern relativiert ... dann ist alles gut.

Nochmal vielen Dank für Deine "Wegweisung".