Herbstworkshop „Flexible Regressionsmodelle“Magdeburg, 22./23. November 2007

Schätzung von ein- und zwei-dimensionalen Perzentilkurven mit

der LMS-Methode

Siegfried Kropf1), Brigitte Peters1), Karl-Otto Dubowy2)

1) Institut f. Biometrie u. Medizinische Informatik, Otto-von-Guericke-Universität Magdeburg

2) Klinik für Angeborene Herzfehler, Herzzentrum NRW Bad Oeynhausen

ProblemstellungGesucht: altersabhängige Referenzwerte für diagnostische Parameter aus großen Stichproben

⇒ Detektion von pathologischen Befunden (häufig Perzentile 90 und 97)

⇒ Umrechnen von Rohwerten in SDS-Werte (Standard Deviation Scores)

HR_peak (1/min)

100

120

140

160

180

200

220

4 9 14 19 24 29 34 39 44 49 54 59 64 69 74Alter in Jahren

HR_peak (1/min)

120

140

160

180

200

220

3 8 13 18 23 28 33 38 43 48 53 58 63 68 73

Alter in Jahren

männlichweiblich Perz.

979075502510

3

Zu berücksichtigende Probleme:

• Parameterwerte werden in der Regel nicht normalverteilt sein.• Vollständig nichtparametrische Perzentilschätzung könnte zu große

Stichproben erfordern.

• Altersabhängigkeit ist vorhanden, aber zumeist nicht linear (z.B. Wachstumsschübe bei Körpergröße).

• Altersgruppen wären andererseits irgendwie künstlich; die Natur hat keine Sprünge.

Vorschlag von Cole und Green (1992):• Erweiterung des Normalmodells: es wird unterstellt, dass für die unter-

suchte Zufallsgröße y = y(t) bei fixem t die Größe yλ normalverteilt ist (Box-Cox-Potenz λ; damit schiefe Verteilungen modellierbar).

• Für die Parameter der resultierenden Verteilung wird angenommen, dass sie stetig und glatt über die Zeit variieren; Schätzung über einen penali-sierten Likelihoodansatz.

Verteilungsmodell bei fixiertem Alter• Zielgröße y soll nur positive Werte annehmen

• Median µ

• yλ soll normalverteilt sein (bzw. im Grenzfall λ = 0 soll ln(y) normal-verteilt sein), ⇒ Einbeziehung ursprünglich schiefer Verteilungen

• betrachten dann Box-Cox-transformierte Größe

bzw.

• Abbildung µ → 0 (zentrierte Normalverteilung für x)

• Standardabweichung von x sei σ; für λ =1 ist σ Variationskoeffizient von y

0λ,λ

1)µ/( λ

≠−= yx

0λ,)µ/ln( == yx

Anwendungen:

• Standard Deviation Score für einen Wert y der Zielgröße:

zur Charakterisierung z.B. eines medizinischen Parameterwertes in einem „Standardformat“

• Berechnung von Perzentilen:

(Zα = α-Quantil der Standardnormalverteilung)

)0λ(σ/)µ/ln(σ/bzw.)0λ(λσ

1)µ/(σ

λ

===≠−== yxzyxz

)0λ(µbzw.)0λ()λσ1(µ )σ(α100

λ/1αα100

α ==≠+= ZeCZC

Berücksichtigen der Altersabhängigkeit

• Die drei Parameter λ, µ und σ werden jetzt als altersabhängige Größen betrachtet (bzw. abhängig von anderer Einflussgröße).

• Man unterstellt einen stetigen und glatten Verlauf.

• Die Zeitverläufe der drei Parameter werden als L(t), M(t) und S(t) bezeichnet und geben der Methoden den Namen „LMS-Methode“.

• Aus der Standardisierungsformel und der für z

angenommenen Standardnormalverteilung lässt sich die Dichte für y und damit die Loglikelihood-Funktion der Stichprobe ableiten:

( ))()(

1)(/ )(

tStLtMyz

tL −=

∑=

⎟⎟⎠

⎞⎜⎜⎝

⎛−−+−−==

n

i

ii

i

iii

ztStM

ytLySML1

2

2))(ln(

)(ln)()ln()π2ln(),,(ll

• Um den stetigen und glatten Verlauf der Funktionen L(t), M(t) und S(t) zu gewährleisten, wird für jede dieser Funktionen bei derMaximisierung der Likelihood-Funktion ein Strafterm eingefügt (penalized likelihood):

αλ, αµ und ασ sind GlättungsparameterGenauigkeit vs. Glattheit

• Zwischen den beobachteten Zeitpunkten wird mittels kubischerSplines interpoliert. Optimierung iterativ über Fisher-Scoring, erste und zweite Ableitungen können explizit angegeben werden.

• Für Glättungsparameter werden Empfehlungen gegeben, die aber nicht immer gut funktionieren. Auch Tests angeboten. Praxis häufig: nach optischem Eindruck, zuerst für µ, dann λ und σ.

{ } { } { } maxd)(α21d)(α

21d)(α

21 2

σ2

µ2

λ →′′−′′−′′− ∫∫∫ ttSttMttLl

Mögliches Problem: unterschiedliche Variabilität über Alter

• Glättung kann in unterschiedlichen Bereichen der unabhängigen Variablen (z.B. Alter) unterschiedlich wichtig sein.

• Vorschlag von Pan und Cole (2004): Zweiphasiges Vorgehen

– erst „normale“ Anpassung der Perzentilkurven,

– dann Alter umskalieren, so dass 50%-Perzentilkurve auf eine Gerade verzerrt wird (falls überhaupt monotoner Verlauf),

– neue Analyse mit umskaliertem Alter.

LMS-Programm nach Cole and Green (1992)Basisversion frei

Anwendungen im KN Angeborene Herzfehler• Erste passive Anwendungen der LMS-Methode in unserem Institut bei

Qualitätssicherung in Pädiatrischer Endokrinologie (Röhl, Mohnike)• Aktuell aktive Normwertermittlungen bzw. Alters- und Geschlechts-

standardisierung für kardiologische Diagnostik, z.B. Laufbandergometrie.

Protokoll zur Laufbandergometrie

Dubowy/Baden 2002

• 1,5 Minuten Stufen; Erholung 2km/h eben• Stufe I 2,5 km/h 0% Steigung• Stufe II 3,0 km/h 3% Steigung• Stufe III 3,5 km/h 6% Steigung• Stufe IV 4,0 km/h 9% Steigung• Stufe V 4,5 km/h 12% Steigung• Stufe VI 5,0 km/h 15% Steigung• Stufe VII 5,5 km/h 18% Steigung• Stufe VIII 6,0 Km/h 21% Steigung• Stufe IX 6,5km/h 21% Steigung etc.

n = 484 weiblich

n = 580 männlich

Alter 3 – 78 Jahre

• Ca. 50 Funktionsparameter betrachtet.

• Werden automatisch bestimmt und in Standardformat gespeichert.

• Kenntnis aller altersabhängigen Normwerte auch für Experten schwer, Umrechnung in Perzentile oder SDS-Werte hilfreich.

• Entsprechende Tabellen getrennt nach Geschlecht erarbeitet.

• Bereitstellung von Umrechnungen / Grafiken über Web vorbereitet,zusammen mit CIO in Göttingen, noch Copyrights zu klären.

• Einbau der Umrechnungen in Ergometrie-Software angedacht.

• In klinischen Studien Einfluss von Störgrößen weitgehend eliminiert ⇒ Biasvermeidung und Varianzreduktion;

im Kindes- und Jugendalter besonders wichtig.

Beispiele:

HR_peak (1/min)

120

140

160

180

200

220

3 8 13 18 23 28 33 38 43 48 53 58 63 68 73

Alter in Jahren

männlichHR_peak (1/min)

100

120

140

160

180

200

220

4 9 14 19 24 29 34 39 44 49 54 59 64 69 74Alter in Jahren

Perz.

979075502510

3

weiblich

Erweiterung auf 2-dimensionale Grundmenge

Frage: Lassen sich Referenzbereiche angeben in Abhängigkeit von zwei Einflussgrößen, z.B. Alter und Größe oder Alter und BMI?

Erste (zaghafte) Versuche!

Problem scheint interessant zu sein, aber sehr rechenintensiv, wenn man nur (ausgewählte) Standardsoftware benutzen will.

Grundansatz ähnlich wie bislang:

bleibt, nur dass die Parameter λ, µ und σ jetzt aus Funktionen zweier unabhängiger Variabler ausgelesen werden: L(t,u), M(t,u) und S(t,u).

)0λ(σ/)µ/ln(σ/bzw.)0λ(λσ

1)µ/(σ

λ

===≠−== yxzyxz

∑=

⎟⎟⎠

⎞⎜⎜⎝

⎛−−+−−=

n

i

iii

ii

iiii

zutSutM

yutLy1

2

2)),(ln(

),(ln),()ln()π2ln(l

Loglikelihood analog:

Glättung jetzt in zwei Dimensionen:

max

dd),(α),(αα),(α21

dd),(α),(αα),(α21

dd),(α),(αα),(α21

2

2

2

σ

22

σσ

2

2

2

σ

2

2

2

µ

22

µµ

2

2

2

µ

2

2

2

λ

22

λλ

2

2

2

λ

→⎪⎭

⎪⎬⎫

⎪⎩

⎪⎨⎧

⎥⎦

⎤⎢⎣

⎡∂∂+⎥

⎦

⎤⎢⎣

⎡∂∂

∂+⎥⎦

⎤⎢⎣

⎡∂∂−

⎪⎭

⎪⎬⎫

⎪⎩

⎪⎨⎧

⎥⎦

⎤⎢⎣

⎡∂∂+⎥

⎦

⎤⎢⎣

⎡∂∂

∂+⎥⎦

⎤⎢⎣

⎡∂∂−

⎪⎭

⎪⎬⎫

⎪⎩

⎪⎨⎧

⎥⎦

⎤⎢⎣

⎡∂∂+⎥

⎦

⎤⎢⎣

⎡∂∂

∂+⎥⎦

⎤⎢⎣

⎡∂∂−

∫

∫

∫

ututSu

utSut

utSt

ututMu

utMut

utMt

ututLu

utLut

utLt

uutt

uutt

uuttl

Glättungsparameter für L, M und S in jeweils zwei Richtungen

• Größere Variabilität zwingt bei gleichbleibenden Stichprobenumfängenzu stärkerer Glättung.

• Im Einzelfall abzuwägen, ob man sich durch zweite Einflussgrößewirklich verbessert oder ob man sich besser anders behilft (z.B. durch dieNutzung von Zielgröße, die bereits auf Körperoberfläche standardisiert ist).

• Man kann durch geeignete Wahl der t- und u- Komponenten der Glättungsparameter Prioritäten setzen hinsichtlich der Genauigkeit der Modellierung des Einflusses der beiden unabhängigen Variablen.

Probleme bei Realisierung mit (bei uns) vorhandener Software:

SAS/STAT bzw. SAS/IML:

• unterstützt zwar Splines, aber nicht zweidimensional (?)

MatLab:

• unterstützt prinzipiell zweidimensionale Splines

• keine Speicherung der Spline-Koeffizienten, nur Ausgabe der fertig gefitteten Daten auf Zielraster ⇒ rein numerische Bildung der Ableitungen, Integrale und Optimierung.

• Ausgangswerte für Interpolation müssen auf Rechteckgitter vorliegen

t (Alter)

u

(Körper-größe)

Beispiel: aus 12 Messwerten werden 70 Gitterpunkte gebildet

⇒ beträchtliche Erhöhung der Dimension des Optimierungsproblems⇒ sehr unattraktive Rechenzeiten (> 24 Stunden; Mehrfachläufe mit

verschiedenen Glättungen erforderlich)

Beispieldarstellung zweidimensionaler Referenzwerte: ein diagnostischer Parameter über Alter und Körpergröße

Median, Ober- und Untergrenze (3 %, 50 %, 97%)

Darstellung hier nur anhand eines Teils der Daten und mit noch zu schwacher Glättung über die Körpergröße

AlterKörpergröße

Diagnostischer Parameter

Für praktische Anwendung wäre Software-Unterstützung erforderlich oder mehrere zweidim. Tabellen (z.B. Parameter über Alter in separaten Tabellen nach Körpergröße

Förderhinweis:

Die Arbeit wurde unterstützt durch das Kompetenznetz Angeborene Herzfehler, gefördert vom Bundesministerium für Bildung und Forschung (Förderkennzeichen: 01G10210)

Literatur

• Cole, T.J. and Green, P.J. (1992). Smoothing reference centile curves: the LMS method and penalized likelihood. Statistics in Medicine 11, 1305-1319.

• Dubowy, K.-O., Bernitzki, ST., Peters, B. (2006). Objektive Belastbar-keit von Patienten mit angeborenen Herzfehlern. Die medizinische Welt 57, 158-163.

• Pan, H. and Cole, T.J. (2004). A comparison of goodness of fit tests forage-related reference ranges. Statistics in Medicine 23, 1749-1765.