VDI/VDE-GMA - Workshop Neue Entwicklungen in der Messunsicherheitsermittlung“ ” 6. und 7. April 2006, Berghotel Gabelbach, Ilmenau Th¨ uringen

Nichtlineare Regression: flexibler aber auch anspruchsvoller als die lineare Regression Andreas Ruckstuhl

Abstract - Deutsch Eine nahe liegende Erweiterung der linearen multiplen Regression wird dadurch erreicht, dass die in den unbekannten Parametern lineare Regressionsfunktion durch eine nichtlineare ersetzt wird. Diese Erweiterung schafft grenzenlose M¨oglichkeiten. Nur schon f¨ ur sigmoide Kurven gibt es unz¨ahlige Funktionsfamilien, die sich z.T. aber kaum unterscheiden. In der Praxis sind deshalb jene Ans¨atze, die sich aus dem Anwendungsgebiet ableiten lassen, am erfolgreichsten. Will man nichtlineare Funktionen an Daten anpassen, so geht man analog wie in der linearen Regressionsrechnung vor. Allerdings stellen sich einige zus¨atzliche Herausforderungen, weil die Sch¨atzgleichungen nicht mehr explizit gegeben sind: (1) Die Gleichungen m¨ ussen iterativ gel¨ost werden und es braucht deshalb geeignete Startwerte. (2) Wenn der Algorithmus nicht konvergiert hat, k¨onnen die L¨osungen grob falsch sein. (3) Die Vertrauensintervalle sind immer approximativ und je nach Datensatz k¨onnen sie sogar unbrauchbar sein, auch wenn der Algorithmus konvergiert hat. Der Vortrag wird anhand von praktischen Beispielen die erw¨ahnten Problemkreise aufzeigen und Hinweise geben, wie man damit umgehen kann. Abstract - Englisch A straightforward extension of the linear regression models can be obtained by replacing the linear regression function by a function which is non-linear in its unknown parameters. This extension opens a wide range of functions to fit. For example, there are many functions for describing sigmoid curves some of which are hardly to distinguish however. Successful fits of non-linear functions are often based on functions which are derived from the theory in the application area. Non-linear functions are fitted to data basically in the same way as linear regression

1

functions. However, there are some challenges to cope with: (1) The estimation equations must be solved iteratively and hence reasonable stating values are needed. (2) If the algorithm does not converge, the results may be completely invalid. (3) The standard confidence intervals are always approximations and, depending on the dataset, may be useless even the algorithm has converged. In the talk, the addressed challenges are introduced using examples from practice and some advices are provided to handle them.

1 Das nichtlineare Regressionsmodell Regression studiert den Zusammenhang zwischen einer Zielgro ¨sse Y und einer oder mehreren erkl¨ arenden (oder Ausgangs-) Variablen x(j) . Das allgemeine Modell lautet   (1) (2) (m) Yi = h xi , xi , . . . , xi ; θ1 , θ2 , . . . , θp + Ei . Dabei ist h eine geeignete Funktion, die von den erkl¨arenden Variablen und von Para(1) (2) (m) metern abh¨angt, die wir zu Vektoren zusammenfassen wollen, x = [xi , xi , . . . , xi ]T und θ = [θ1 , θ2 , . . . , θp ]T . Die unstrukturierten Abweichungen von der Funktion h werden durch zuf¨allige Fehler Ei beschrieben. F¨ ur die Verteilung dieser zuf¨alligen Fehler wird die Normalverteilung angenommen, also
Ei ∼ kN 0, σ 2 ,

unabh¨angig.

In der (multiplen) linearen Regression werden Funktionen h betrachtet, die linear sind in den Parametern θj ,   (1) (2) (m) (1) (2) (p) h xi , xi , . . . , xi ; θ1 , θ2 , . . . , θp = θ1 xei + θ2 xei + . . . + θp xei , wobei die xe(j) beliebige Funktionen der urspr¨ unglichen erkl¨arenden Variablen x(j) sein k¨onnen. (Die Parameter werden dort u ¨blicherweise mit βj statt θj bezeichnet.)

In der nichtlinearen Regression werden Funktionen h untersucht, die sich nicht als lineare Funktionen in den Parametern schreiben lassen. Oft wird eine solche Funktion aus der Theorie abgeleitet. Es bestehen im Prinzip unbeschr¨ankte M¨oglichkeiten, den deterministischen Teil des Modells anzusetzen. Wie wir sehen werden, wird diese Flexibilit¨at erkauft durch einen gr¨osseren Aufwand, statistische Aussagen zu gewinnen.

2

Geschwindigkeit

Geschwindigkeit

200

150

100

50 0.0

0.2

0.4

0.6

0.8

1.0 Konzentration

Konzentration

Abbildung 1: Beispiel Puromycin. (a) Daten. • behandeltes Enzym; △ unbehandelt) und (b)

typischer Verlauf der Regressionsfunktion.

⊲ Beispiel Puromycin. Die Geschwindigkeit, mit der eine enzymatische Reaktion abl¨auft, h¨angt von der Konzentration eines Substrates ab. Gem¨ass den Angaben in [1] wurde untersucht, wie eine Behandlung des Enzyms mit einer weiteren Substanz namens Puromycin diese Reaktionsgeschwindigkeit beeinflusst. Als Zielvariable wurde die Anfangsgeschwindigkeit der Reaktion gew¨ahlt, welche u ¨ ber Radioaktivit¨at gemessen wird. (Die Einheit der Zielvariablen ist Anzahl/min2 ; die Anzahl Anschl¨age in einem Geigerz¨ahler pro Zeiteinheit misst ja die Quantit¨at der vorhandenen Substanz, und die Reaktionsgeschwindigkeit ist proportional zu deren Ver¨anderung pro Zeiteinheit) Der Zusammenhang der Zielgr¨osse mit der Substrat-Konzentration x (in ppm) wird beschrieben durch die Michaelis-Menten-Funktion h (x; θ) =

θ1 x . θ2 + x

F¨ ur unendlich grosse Substratkonzentration (x → ∞) ergibt sich die “asymptotische” Geschwindigkeit θ1 . Es wurde vermutet, dass diese Gr¨osse durch das Hinzuf¨ ugen von Puromycin beeinflusst wird. Das Experiment wurde deshalb einmal mit dem mit Puromycin behandelten Enzym und einmal mit dem unbehandelten Enzym durchgef¨ uhrt. Abbildung 1 zeigt das Resultat. In diesem Abschnitt werden die Daten des behandelten Enzyms benutzt. ⊲ Beispiel Sauerstoffverbrauch. Um den biochemischen Sauerstoffverbrauch zu bestimmen, werden Flusswasserproben mit gel¨osten organischen N¨ahrstoffen, mit anorganischen Materialien und mit gel¨ostem Sauerstoff angereichert und in verschiedene Flaschen abgef¨ ullt (Marske, 1967, siehe [1]). Jede Flasche wird dann mit einer Mischkultur von Mikroorganismen geimpft und verschlossen in eine Klimakammer mit konstanter Temperatur gestellt. Die Flaschen werden periodisch ge¨offnet und nach gel¨ostem Sauerstoffgehalt analysiert. Daraus wird der biochemische Sauerstoffverbrauch [mg/l] berechnet. Das verwendete Modell, das den kumulierten biochemischen Sauerstoffverbrauch Y mit der Inkubinationszeit x in Verbindung bringt, basiert auf exponentiellem Abfall der Zuw¨achse,

3

20

• •

16

•

•

4

5

Sauerstoffverbrauch

Sauerstoffverbrauch

18

14 12 •

10 8

• 1

2

3

6

7 Tage

Tage

Abbildung 2: Beispiel Sauerstoffverbrauch. (a) Daten und (b) typischer Verlauf der Regressi-

onsfunktion

was zu h (x, θ) = θ1 1 − e−θ2 x




f¨ uhrt. Abbildung 2 zeigt die Daten und die zu verwendende Regressionsfunktion. ⊲ Beispiel aus der Membrantrenntechnologie [2]. Das Verh¨altnis von protonierten zu deprotonierten Carboxylgruppen in den Poren von Cellulosemembranen ist vom pH-Wert x der Aussenl¨osung abh¨angig. Die Protonierung des Carboxylkohlenstoffatoms kann mit 13 C-NMR erfasst werden. Wir nehmen an, dass der Zusammenhang mit der erweiterten Henderson-Hasselbach” Gleichung“ f¨ ur Polyelektrolyte beschrieben werden kann,   θ1 − y = θ3 + θ4 x , log10 y − θ2 wobei die unbekannten Parameter θ1 , θ2 und θ3 > 0 und θ4 < 0 sind. Aufl¨osung nach y f¨ uhrt zum Modell θ1 + θ2 10θ3 +θ4 xi Yi = h(xi ; θ) + Ei = + Ei . 1 + 10θ3 +θ4 xi ur ein willk¨ urlich gew¨ahltes θ in Abbildung 3(b) und die Daten Die Regressionsfunktion h(xi , θ) ist f¨ sind in Abbildung 3(a) dargestellt.

Linearisierbare Regressionsfunktionen. Einige nichtlineare Regressionsfunktionen lassen sich durch Transformationen der Zielgr¨osse und der erkl¨arenden Variablen linearisieren. Beispielsweise l¨asst sich eine Potenzfunktion, h (x; θ) = θ1 xθ2

4

162

delta

delta

163

161 160

(a) 2

4

6

8

10

(b)

12

pH

pH

Abbildung 3: Beispiel Membrantrenntechnologie. Links: Daten; rechts: ein typischer Verlauf der

Regressionsfunktion.

zu einer (in den Parametern) linearen Funktion transformieren, ln (h (x; θ)) = ln (θ1 ) + θ2 ln (x) = β0 + β1 xe ,

wobei β0 = ln (θ1 ), β1 = θ2 und xe = ln (x) ist. Die Regressionsfunktion h nennen wir linearisierbar, wenn sie durch Transformationen der Argumente und eine monotone Transformation des Resultats in eine in den (unbekannten) Parametern lineare Funktion verwandelt werden kann. Eine lineare Regression mit der linearisierten Regressionsfunktion beruht im genannten Beispiel auf dem statistisch vollst¨andigen Modell ln (Yi ) = β0 + β1 xei + Ei ,

wobei die Zufallsfehler Ei alle der gleichen Normalverteilung folgen. Transformieren wir dieses Modell zur¨ uck, so erhalten wir ei Yi = θ1 · xθ2 · E

ei = exp (Ei ). Die Fehler E ei , i = 1, . . . , n wirken nun multiplikativ und sind lognormit E mal verteilt! Die Annahme u ¨ber die Zufallsabweichung ver¨andert sich also recht drastisch gegen¨ uber einem Modell, das direkt auf h beruht, Yi = θ1 · xθ2 + Ei∗

5

mit Zufallsabweichungen Ei∗ , die wie u ¨blich additiv wirken und einer spezifischen Normalverteilung folgen.

Eine Linearisierung der Regressionsfunktion ist deshalb nur dann angebracht, wenn damit die Annahmen u ¨ber die Zufallsabweichungen besser befriedigt werden k¨onnen – im Beispiel, falls tats¨achlich die Fehler eher multiplikativ als additiv wirken und lognormal statt normal verteilt sind. Diese Annahmen m¨ ussen mit Residuen-Analyse gepr¨ uft werden.

2 Methodik zur Sch¨ atzung der Parameter Um Sch¨atzungen f¨ ur die Parameter θ = [θ1 , θ2 , . . ., θp ]T zu erhalten, wenden man, wie in der linearen Regresionsrechnung, das Prinzip der Kleinsten Quadrate an. Die Summe der quadrierten Abweichungen S(θ) :=

n X

(yi − ηi (θ))2

mit ηi (θ) := h (xi ; θ)

(1)

i=1

soll also minimiert werden. Die Schreibweise, die h (xi ; θ) durch ηi (θ) ersetzt, ist sinnvoll, weil nach der Messung oder Beobachtung die Daten [xi , yi] gegeben sind und nun die Parameter θ zu bestimmen bleiben. Leider lassen sich das Minimum der Quadratsumme und damit die Sch¨atzung nicht wie in der linearen Regression explizit angeben. Iterative numerische Verfahren helfen weiter. Die Grundidee, die hinter dem u ¨blichen Algorithmus steckt, soll hier skizziert werden. Sie bilden auch die Basis f¨ ur die einfachste Art, Tests und Vertrauensbereiche herzuleiten. Geometrischen Veranschaulichung. Die beobachteten Werte Y = [Y1 , Y2, . . . , Yn ]T legen einen Punkt im n-dimensionalen Raum fest. Das Gleiche gilt f¨ ur die “Modellwerte” T η (θ) = [η1 (θ) , η2 (θ) , . . . , ηn (θ)] f¨ ur gegebenes θ. Achtung! Die u ¨ bliche geometrische Veranschaulichung von Daten, die beispielsweise in der multivariaten Statistik grundlegend ist, betrachtet die Beobachtungen, die durch m Variable x(j) , j = 1, 2, . . . , m, festgelegt sind, als Punkte im m-dimensionalen Raum. Hier aber betrachten wir die Y - und η-Werte aller n Beobachtungen als Punkte im n-dimensionalen Raum.

6

Leider h¨ort unsere Anschauung bei drei Dimensionen, also bei drei Beobachtungen auf. Versuchen wir es also f¨ ur eine solche eingeschr¨ankte Situation, zuerst bei der einfachen linearen Regression. Wie gesagt, die beobachteten Werte Y = [Y1 , Y2, Y3 ]T legen einen Punkt im 3-dimensionalen Raum fest. F¨ ur gegebene Parameter β0 = 5 und β1 = 1
k¨onnen wir die Modellwerte ηi β = β0 + β1 xi ausrechnen und den entsprechenden Vek
tor η β = β0 1 + β1 x ebenfalls als Punkt darstellen. Nun fragen wir, wo alle Punkte liegen, die durch Variation der Parameter erreicht werden k¨onnen. Sie sind die m¨oglichen Linearkombinationen der beiden Vektoren 1 und x und bilden deshalb die Ebene die “durch 1 und x aufgespannt” wird. Beim Sch¨atzen der Parameter nach dem Prinzip der Kleinsten Quadrate wird, geometrisch ausgedr¨ uckt, der quadrierte Abstand zwischen
Y und η β minimiert. Gesucht ist also der Punkt auf der Ebene, der den kleinsten Abstand zu Y hat. Dieser wird auch die Projektion von Y auf die Ebene genannt. Die Parameterwerte, die diesem Punkt ηb entsprechen, sind dann die gesch¨atzten Parameterwerte βb = [βb0 , βb1 ]T .

Soll nun eine nichtlineare Funktion, z.B. h (x; θ) = θ1 exp (1 − θ2 x), an die gleichen drei Beobachtungen angepasst werden. Wir k¨onnen uns wieder fragen, wo alle Punkte η (θ) liegen, die durch Variation der Parameter θ1 und θ2 erreicht werden k¨onnen. Sie liegen auf einer zweidimensionalen, gekr¨ ummten Fl¨ache (im Folgenden Modellfl¨ ache genannt) im dreidimensionalen Raum. Das Sch¨atzproblem besteht wieder darin, den Punkt ηb auf der Modellfl¨ache zu finden, der Y am n¨achsten liegt. Die Parameterwerte, die diesem Punkt ηb entsprechen, sind dann die gesch¨atzten Parameterwerte θb = [θb1 , θb2 ]T .

Die Hauptidee des u ¨blichen Algorithums’ zur Minimierung der Summe der quadrierten Abweichungen (vgl. (1)) l¨auft wie folgt: Wenn ein vorl¨aufig bester Wert θ(ℓ) vorliegt,
approximiert man die Modellfl¨ache durch die Ebene, die die Fl¨ache im Punkt η θ(ℓ) =
h x; θ(ℓ) ber¨ uhrt. Nun sucht man den Punkt in dieser Ebene, der am n¨achsten bei Y liegt. Das l¨auft auf die Sch¨atzung in einem linearen Regressionsproblem hinaus. Dieser neue Punkt liegt auf der Ebene, aber nicht auf der Fl¨ache, die dem nicht-linearen Problem entspricht. Er legt aber einen Parametervektor θ(ℓ+1) fest, und mit diesem geht man in die n¨achste Iterations-Runde. Um die approximierende Ebene zu bestimmen, brauchen wir die partiellen Ableitungen ∂ηi (θ) (j) Ai (θ) := , ∂θj

7

die wir zu einer n×p-Matrix A zusammenfassen k¨onnen. Die Approximation der Modellfl¨ache η (θ) durch die Tangentialebene“ in einem Parameterwert θ∗ lautet ” (1)

(p)

ηi (θ) ≈ ηi (θ ∗ ) + Ai (θ∗ ) (θ1 − θ1∗ ) + ... + Ai (θ∗ ) (θp − θp∗ ) oder, in Matrixschreibweise, η (θ) ≈ η (θ ∗ ) + A (θ∗ ) (θ − θ∗ ) .

(2)

Wenn wir nun den Zufallsfehler wieder hinzuf¨ ugen, erhalten wir ein lineares Regressionsmodell Y − η (θ ∗ ) = A (θ ∗ ) β + E (3) mit den vorl¨aufigen Residuen“ Yi −ηi (θ∗ ) als Zielgr¨osse, den Spalten von A als Regresso” ren und den Koeffizienten βj = θj −θj∗ (ein Modell ohne Achsenabschnitt β0 ). Der Gauss-Newton-Algorithmus besteht darin, ausgehend von einem Startwert θ(0) uhrte lineare Regressionsproblem f¨ ur θ∗ = θ (0) zu l¨osen, um eif¨ ur θ das gerade eingef¨ ne Korrektur β und daraus einen verbesserten Wert θ(1) = θ(0) + β zu erhalten. F¨ ur diesen wird wieder das approximierende Modell ausgerechnet, also die vorl¨aufigen Resi”

duen“ Y − η θ(1) und die partiellen Ableitungen A θ(1) bestimmt, und daraus ergibt sich θ2 . Mit diesem Iterationsschritt wird so lange weitergefahren, bis die Korrektur β vernachl¨assigbar wird. (Weitere Details finden Sie in [1].) Es kann nicht garantiert werden, dass dieses Verfahren tats¨achlich das Minimum der Quadratsumme findet. Die Chancen daf¨ ur stehen besser, je besser sich die pdimensionale Modellfl¨ache im Minimum θb = (θb1 , . . . , θbp )T durch eine p-dimensinale Ebene“ lokal approximieren l¨asst, und je n¨aher der Startwert θ (0) zur gesuchten L¨osung ” ist. Komfortable Algorithmen bestimmen die Ableitungsmatrix A numerisch. In komplexeren Problemen kann die numerische N¨aherung ungen¨ ugend sein und Konvergenzprobleme verursachen. Dann ist es von Vorteil, wenn analytische Ausdr¨ ucke f¨ ur die ersten partiellen Ableitungen angegeben werden k¨onnen. Damit kann die Ableitungsmatrix numerisch zuverl¨assiger bestimmt werden und das Verfahren konvergiert eher. Ein iteratives Verfahren ben¨otigt Startwerte, damit es u ¨berhaupt angewandt werden kann. Gute Startwerte helfen, dass das iterative Verfahren schneller und sicherer die L¨osung findet. Einige M¨oglichkeiten, diese mehr oder weniger einfach zu gewinnen, werden hier kurz vorgestellt.

8

•

Geschwindigkeit

1/Geschwindigkeit

0.015 • • •

0.010

0.005

••• 0

••

• •

150 • • • •

100

• •

• 50

10

• •

• •

200

0.020

20

30

40

50

• 0.0

1/Konzentration

0.2

0.4

0.6

0.8

1.0

Konzentration

Abbildung 4: Beispiel Puromycin. Links: Regressionsgerade im linearisierten Problem. Rechts:

Regressionsfunktion h (x; θ) f¨ ur die Startwerte θ = θ(0) ( Sch¨ atzung θ = θb (——–).

) und f¨ ur die Kleinste-Quadrate-

Wie schon in der Einleitung bemerkt, stammen die nichtlinearen Modelle vielfach aus ¨ theoretischen Uberlegungen im jeweiligen Anwendungsgebietes. Bestehen schon Vorkenntnisse aus ¨ahnlichen Experimenten, so k¨onnen diese verwendet werden, um Startwerte zu gewinnen. Um sicher zu gehen, dass der gew¨ahlte Startwert passt, ist es ratsam, die Regressionsfunktion h (x; θ) f¨ ur verschiedene m¨ogliche Startwerte θ = θ 0 graphisch mit den Daten zusammen darzustellen (z.B. so wie in Abbildung 4, rechts). Manchmal ist man wegen der Verteilung der Fehler gezwungen, in Modellen mit linearisierbaren Regressionsfunktionen bei der nichtlinearen Form zu verbleiben. Das linearisierte Modell kann aber Startwerte liefern. ⊲ Beispiel Puromycin. In diesem Beispiel ist die Regressionsfunktion linearisierbar: Die Kehrwerte der beiden Variablen erf¨ ullen 1 1 θ2 1 1 ye = ≈ = + = β0 + β1 xe . y h (x; θ) θ1 θ1 x

Die Kleinste-Quadrate-L¨osung f¨ ur dieses modifizierte Problem ist βb = [βb0 , βb1 ]T = (0.00511, 0.000247)T (Abbildung 4 (a)). Das liefert die Startwerte (0)

(0)

θ1 = 1/βb0 = 196 ,

θ2 = βb1 /βb0 = 0.048 .

Oft ist es hilfreich, die analytischen Besonderheiten der Regressionsfunktion zu untersuchen. Wertet man z.B. die Funktion h (xi ; θ) = θ1 + θ2 exp (−θ3 x), θ3 > 0, an der Stelle x = 0 aus, so erh¨alt man h (0; θ) = θ1 +θ2 . Damit ist eine visuelle Sch¨atzung des yWertes f¨ ur x = 0 ein guter Startwert f¨ ur θ1 +θ2 . Da f¨ ur x → ∞ die Funktion h (x; θ) → θ1

9

2

(a) 163

1

δ

y

162 0

161 −1 160

(b)

−2 2

4

6

8

10

12

2

pH

4

6

8

10

12

pH

Abbildung 5: Beispiel aus der Membrantrenn-Technologie. (a) Regressionsgerade, die zur Be-

stimmung der Startwerte f¨ ur θ3 und θ4 gebraucht wird. (b) Regressionsfunktion h (x; θ) f¨ ur die (0) b ) und f¨ ur die der Kleinste-Quadrate-Sch¨ atzung θ = θ (——–). Startwerte θ = θ (

geht, kann der an der maximalen x-Stelle beobachtete y-Wert als Startwert f¨ ur θ1 ben¨ utzt werden. Beides zusammen ergibt dann einen Startwert f¨ ur θ2 . ⊲ Beispiel Im Beispiel Puromycin k¨onnen wir auch auf eine andere Art zu Startwerten gelangen: θ1 ist der y-Wert f¨ ur x = ∞. Da die Regressionsfunktion monoton steigend ist, k¨onnen wir den maximalen yi -Wert oder einen visuell bestimmten “asymptotischen Wert” θ10 = 207 als Startwert f¨ ur θ1 ben¨ utzen. Der Parameter θ2 ist der x-Wert, bei dem ydie H¨alfte des asymptotischen Wertes θ1 erreicht. Das ergibt θ20 = 0.06. ⊲ Beispiel aus der Membrantrenn-Technologie. Lassen wir im Beispiel aus der MembrantrennTechnologie x → ∞ gehen, so geht h (x; θ) → θ1 (da θ4 < 0); f¨ ur x → −∞ geht h (x; θ) → θ2 . Aus Abbildung 3(a) und den Daten geht hervor, dass θ1 ≈ 163.7 und θ2 ≈ 159.5 ist. Sind θ1 und θ2 bekannt, so kann man die Regressionsfunktion linearisieren durch ! (0) θ1 − y = θ3 + θ4 x . ye := log10 (0) y − θ2 Man spricht von einer bedingt linearisierbaren Funktion. Die lineare Regression f¨ uhrt zu den (0) (0) Startwerten θ3 = 1.83 und θ4 = −0.36. b b b Mit diesen Startwerten konvergiert der Algorithmus  zur  L¨osung θ1 = 163.7, θ2 = 159.8, θ3 = 2.67
(0) b b und θ 4 = −0.512. Die Funktionen h ·; θ und h ·; θ sind in Abbildung 5(b) eingezeichnet.

Die Eigenschaft der bedingten Linearit¨at von Funktionen kann auch dazu benutzt werden, einen dieser Situation speziell angepassten Algorithmus zu entwickeln (siehe z.B. [1]).

10

3 Gen¨ aherte Vertrauensintervalle Die Sch¨atzung θb liefert den Wert von θ, der optimal zu den Daten passt. Nun fragen wir, welche Parameterwerte θ mit den Beobachtungen vertr¨aglich sind. Der Vertrauensbereich ist die Menge all dieser Werte. F¨ ur einen einzelnen Parameter θj wird der Vertrauensbereich zum Vertrauensintervall oder Konfidenzintervall. Die Resultate, die nun folgen, beruhen darauf, dass die Sch¨atzung θb asymptotisch multivariat normalverteilt ist. F¨ ur einen einzelnen Parameter f¨ uhrt das zu einem “z-Test” und zum entsprechenden Vertrauensintervall; f¨ ur mehrere Parameter kommt der entsprechende Chiquadrat-Test zum Zug und liefert elliptische Vertrauensbereiche. Die asymptotischen Eigenschaften der Sch¨atzung k¨onnen aus der linearen Approximation hergeleitet werden. Das Problem der nichtlinearen Regression ist ja n¨aherungsweise gleich dem in (3) erw¨ahnten linearen Regressionsproblem, wenn der Parametervektor θ∗ , der f¨ ur die Linearisierung verwendet wird, nahe bei der L¨osung liegt. Im b ur β im linearen Problem exakt = 0 – sonst w¨are es L¨osungspunkt θ ist die L¨osung f¨ nicht die L¨osung. Die Standardfehler der Koeffizienten β – und allgemeiner die Kovarianzmatrix von βb – geben aber n¨aherungsweise die entsprechenden Gr¨ossen f¨ ur θb wieder. Etwas genauer: Die Standardfehler geben ja die Unsicherheiten wieder, die durch die Zufallsschwankungen der Daten erzeugt werden. Die vorliegenden Daten haben zum Sch¨atzwert θb gef¨ uhrt. W¨aren die Daten etwas anders ausgefallen, dann w¨are θb immer noch ungef¨ahr richtig, also – so nehmen wir an – gut genug f¨ ur die Linearisierung. Die Sch¨atzung von β f¨ ur den neuen Datensatz w¨ urde also so weit vom Sch¨atzwert f¨ ur den vorliegenden Daten weg liegen, wie es der Verteilung der Parameter im linearisierten Problem entspricht. ¨ Aus dieser Uberlegung folgt: Asymptotisch ist die Kleinste-Quadrate-Sch¨atzung θb konsistent und normalverteilt,   V (θ) a b θ ∼ N θ, , n

mit asymptotischer Kovarianzmatrix V (θ) = σ 2 (A (θ)T A (θ))−1 , wobei A (θ) die n × p Matrix der partiellen Ableitungen ist (vgl. (2)). Um die Kovarianzmatrix V (θ) explizit zu bestimmen, wird A (θ) an der Stelle θb statt an der unbekannten Stelle θ berechnet,

11

und f¨ ur die Fehlervarianz σ 2 wird die u ¨bliche Sch¨atzung eingesetzt,     −1 T 2 [ b V (θ) = σb A θ A θb mit

n     2 1 1 X b σb = S θ = yi − ηi θb . n−p n − p i=1 2

Damit ist die Verteilung der gesch¨atzten Parameter n¨aherungsweise bestimmt, und daraus lassen sich wie in der linearen Regression Standardfehler und Vertrauensintervalle herleiten, ebenso Vertrauens-Ellipsen (oder -Ellipsoide), wenn mehrere Parameter gemeinsam betrachtet werden. Der Nenner n − p in σb2 wurde in der linearen Regression eingef¨ uhrt, um die Sch¨atzung erwartungstreu zu machen. Tests und Vertrauensintervalle wurden nicht mit der Normalund Chiquadrat-Verteilung bestimmt, sondern mit der t- und F-Verteilung. Damit wurde ber¨ ucksichtigt, dass die Sch¨atzung von σ 2 eine zus¨atzliche Zufallsschwankung bewirkt. Auch wenn die Verteilungen nicht mehr exakt gelten, so werden die N¨aherungen doch genauer, wenn man dies bei der nichtlinearen Regression ebenfalls tut. Asymptotisch geht der Unterschied gegen null. ⊲ Beispiel aus der Membrantrenntechnologie. Eine Computer-Ausgabe f¨ ur das Beispiel aus der Membrantrenn-Technologie zeigt Tabelle 1. Die Sch¨atzungen der Parameter stehen in der Kolonne “Value”, gefolgt von den gesch¨atzten approximativen Standardfehler und den Teststatistiken (“t value”), die approximativ tn−p -verteilt sind. In der letzten Zeile wird die gesch¨atzte Standardabweichung σb der Zufallsfehler Ei angegeben. Aus diesen Angaben k¨onnen wie in der linearen Regression die Vertrauensintervalle f¨ ur die Parameter bestimmt werden: Das approximative 95%-Vertrauensintervall f¨ ur den Parameter θ1 ist t35 163.706 ± q0.975 · 0.1262 = 163.706 ± 0.256 .

¨ ⊲ Beispiel Puromycin. Zur Uberpr¨ ufung eines Einflusses der Behandlung des Enzyms mit Puromycin von der postulierten Form kann ein gemeinsames Modell f¨ ur die Daten mit und ohne Behandlung folgendermassen formuliert werden: Yi =

(θ1 + θ3 zi )xi + Ei . θ2 + θ4 zi + xi

Dabei ist z die Indikatorenvariable f¨ ur die Behandlung (zi = 1, wenn behandelt, sonst =0). Tabelle 2 zeigt, dass der Parameter θ4 nicht signifikant von 0 verschieden ist, denn der t-Wert von t19 1.44 ist kleiner als die kritische Grenze q0.975 = 2.09. Die Behandlung hat aber einen eindeutigen

12

Parameters: Value θ1 163.706 θ2 159.784 θ3 2.675 θ4 -0.512

Std. Error 0.1262 0.1595 0.3813 0.0703

t value 1297.21 1002.03 7.02 -7.28

Residual standard error: 0.293137 on 35 degrees of freedom

Tabelle 1: Computer-Ausgabe f¨ ur das Beispiel aus der Membrantrenn-Technologie Einfluss, der sich durch θ3 ausdr¨ uckt; das 95% Vertrauensintervall u ¨ berdeckt den Bereich 52.398 ± 9.5513 · 2.09 = [32.4, 72.4]. Parameters: Value θ1 160.286 θ2 0.048 θ3 52.398 θ4 0.016

Std. Error 6.8964 0.0083 9.5513 0.0114

t value 23.24 5.76 5.49 1.44

Residual standard error: 10.4 on 19 degrees of freedom

Tabelle 2: Computer-Ausgabe f¨ ur das Beispiel Puromycin

Neben den Parametern ist oft der Funktionswert h (x0 , θ) f¨ ur ein beliebiges x0 von In
teresse. In der linearen Regression wird der Funktionswert h x0 , β = xT0 β durch xT0 βb gesch¨atzt, und das (1 − α)-Vertrauensintervall daf¨ ur ist q tn−p xTo βb ± σb xTo (X T X)−1 xo q1−α/2 . ¨ Durch analoge Uberlegungen und asymptotische N¨aherung kann man Vertrauensintervalle f¨ ur den Funktionswerte h (x0 ; θ) f¨ ur nicht lineare h angeben. Wird die Funktion η0 (θ) := h (x0 , θ) an der Stelle θ linear approximiert, so erh¨alt man   ∂h (xo , θ) ηo θb ≈ ηo (θ) + aTo (θb − θ) mit ao = . ∂θ

(Wenn x0 gleich einem beobachteten xi ist, dann ist a0 gleich der entsprechenden Zeile der Matrix A aus (2).) Das Vertrauensintervall f¨ ur den Funktionswert h (x0 , θ) ist dann approximativ r        −1 T t n−p T b b h xo , θ ± q σbx mit σbx = σb ab A θ A θb abo . (4) 1−α/2

o

o

13

o

•

•• •

2 • • • 1

0

• •

• •

••

•

•

• •

• •

30

• • •

Sauerstoffverbrauch

log(PCB Konzentration)

3

•

• • •

20

•

• •

15

•

•

10

•

5

• • 1.0

25

0 1.2

1.4

1.6 1.8 Jahre^(1/3)

2.0

2.2

0

2

4 Tage

6

8

Abbildung 6: Links: Vertrauensband f¨ ur eine gesch¨ atzte Gerade bei einem linearen Problem.

Rechts: Vertrauensband f¨ ur die gesch¨atzte Kurve h (x, θ) im Beispiel Sauerstoffgehalt.

In dieser Formel wurden wieder die unbekannten Gr¨ossen durch ihre Sch¨atzungen ersetzt. ur beliebiges Der Ausdruck f¨ ur das (1 − α)-Vertrauensintervall f¨ ur ηo (θ) := h (xo , θ) gilt f¨ xo . Es ist wie in der linearen Regression nahe liegend, die Grenzen dieses Intervalls als Funktion von xo als “Vertrauensband” aufzuzeichnen, wie dies Abbildung 6 f¨ ur die beiden Beispiele Puromycin und Sauerstoffverbrauch zeigt. Vertrauensb¨ander f¨ ur lineare und nichtlineare Regressionsfunktionen verhalten sich verschieden: Bei linearen Funktionen ist das Vertrauensband beim Schwerpunkt der Ausgangs-Variablen am engsten und wird gegen aussen allm¨ahlich breiter (vgl. Abbildung 6, links). Im nichtlinearen Fall k¨onnen die B¨ander beliebiger sein. Weil die Funktionen in den Beispielen durch den Nullpunkt gehen m¨ ussen, schrumpft dort das Intervall zu einem Punkt. Beide Modelle haben eine horizontale Asymptote und deshalb wird das Band f¨ ur grosse x eine konstante Breite erreichen(vgl. Abbildung 6, rechts) .

14

4 Genauere Vertrauensintervalle Die Qualit¨at der approximativen Vertrauensbereiche ist stark von der Qualit¨at der linearen Approximation abh¨angig. Ebenfalls werden die Konvergenzeigenschaften der Optimierungsalgorithmen durch die Qualit¨at der linearen Approximation beeinflusst. Mit etwas gr¨osserem Rechenaufwand l¨asst sich die Linearit¨at grafisch u ufen, und gleich¨berpr¨ zeitig erh¨alt man genauere Vertrauensintervalle. Um eine Nullhypothese θ = θ∗ f¨ ur den ganzen Parametervektor oder auch θj = θj∗ f¨ ur eine einzelne Komponente zu testen, kann man, wie in der linearen Regression, den F-Test zum Vergleich von Modellen verwenden. Man vergleicht dabei die Quadratsumme   ∗ S (θ ), die sich unter der Nullhypothese ergibt, mit der Quadratsumme S θb . (F¨ ur n → ∞ stimmt der F-Test mit dem so genannten Likelihood-Quotienten-Test u ¨berein, und die Quadratsumme ist, bis auf eine Konstante, gleich der Log-Likelihood.) Zun¨achst wollen wir eine Nullhypothese θ = θ∗ f¨ ur den ganzen Parametervektor betrachten. Die Teststatistik ist   ∗ ) − S θb a S (θ n−p   T = ∼ Fp,n−p . p S θb Daraus erh¨alt man als Vertrauensbereich   n θ S (θ) ≤ S θb 1 + F

p n−p

q

o

p,n−p wobei q = q1−α das (1 − α)-Quantil der F-Verteilung mit p und n − p Freiheitsgraden ist.

In der linearen Regression erh¨alt man genau den gleichen Vertrauensbereich, wenn man die (multivariate) Normalverteilung der Sch¨atzung βb ben¨ utzt. Im nichtlinearen Fall sind die Ergebnisse verschieden. Der Bereich, der auf dem F-Test beruht, st¨ utzt sich nicht auf die lineare Approximation in (2) ab und ist deshalb (viel) exakter. Falls p = 2 ist, k¨onnen wir den exakten Vertrauensbereich finden, indem wir S (θ) auf einem Gitter von θ-Werten berechnen und durch Interpolation die Grenzen dieses Bereichs bestimmen wie das f¨ ur Kontur-Plots gel¨aufig ist. In Abbildung 7 sind die Konturen zusammen mit den elliptischen Bereichen, die sich aus der linearen Approximation ergeben, f¨ ur die Beispiele Puromycin (links) und Sauerstoffverbrauch (rechts) wiedergegeben.

15

0.10

10

0.09

8

0.08 6

θ2

θ2

0.07

4

0.06 2

0.05 0.04

0 190

200

210

220

230

240

0

10

θ1

20

30

40

50

60

θ1

Abbildung 7: Nominale 80 und 95% Likelihood-Konturen (——) und die Vertrauensellipsen aus

der asymptotischen Approximation (– – – –). + bezeichnet die Kleinste-Quadrate L¨ osung. Im ¨ Beispiel Pyromycin (links) ist die Ubereinstimmung gut und im Beispiel Sauerstoffverbrauch (rechts) schlecht.

F¨ ur p > 2 gibt es keine Kontur-Plots. Grafische Hilfsmittel, die auch f¨ ur h¨ohere Dimen¨ sionen funktionieren, beruhen auf den folgenden Uberlegungen. Es soll gepr¨ uft werden, ob ein einzelner Parameter θk gleich einem bestimmten Wert ∗ ¨ θk sein kann. Uber die u ¨brigen Parameter macht eine solche Nullhypothese keine Aussage. Das Modell, das der Nullhypothese entspricht und am besten zu den Daten passt, ist durch eine Kleinste-Quadrate-Sch¨atzung der u ¨brigen Parameter bei festem θk = θk∗ bestimmt. Es wird also S (θ1 , . . . , θk∗ , . . . , θp ) minimiert in Bezug auf alle θj , j 6= k. Das Minimum bezeichnen wir mit Sek und die Werte θj , die zu ihm f¨ uhren, mit θej . Beide Gr¨ossen h¨angen von θk∗ ab. Wir schreiben deshalb Sek (θk∗ ) und θej (θk∗ ).

Die F-Teststatistik f¨ ur den Test “θk = θk∗ ” ist

  Sek (θk∗ ) − S θb   Tek = (n − p) . S θb

Sie hat (gen¨ahert) eine F1,n−p -Verteilung.

F1,n−p Ein Vertrauensintervall erh¨alt man daraus, indem man die Gleichung Tek = q0.95 nu∗ merisch nach θk aufl¨ost.Sie hat eine L¨osung, die kleiner als θbk ist, und eine, die gr¨osser ist.

16

In der linearen Regression und im vorhergehenden Kapitel haben wir Tests und Vertrauensintervalle aus einer Testgr¨osse ausgerechnet, die einer t-Verteilung folgt (t-Test f¨ ur die Koeffizienten). Ist das ein anderer Test? Es stellt sich heraus, dass die Teststatistik des t-Tests in der linearen Regression in die Teststatistik des F-Tests u ¨bergeht, wenn man sie quadriert, und die beiden Tests sind ¨aquivalent. In der nichtlinearen Regression ist der F-Test nicht ¨aquivalent mit dem im letzten Kapitel besprochenen t-Test. Aber wir k¨onnen den F-Test in einen t-Test verwandeln, der genauer ist als der des letzten Kapitels: Aus der Teststatistik des F-Tests ziehen wir die Wurzel und versehen diese mit dem Vorzeichen von θbk − θk∗ , r   ∗ e S (θ ) − S θb k   k ∗ ∗ b Tk (θk ) := sign θ k − θk . (5) σb   (sign (a) bezeichnet das Vorzeichen von a, und es ist σb2 = S θb /(n − p).) Diese Teststatistik ist (gen¨ahert) tn−p -verteilt. Im linearen Regressionsmodell ist Tk , wie erw¨ahnt, gleich der Teststatistik des u ¨ blichen t-Tests, θbk − θ∗ Tk (θk∗ ) =  k . se θbk

Wir k¨onnen mit dieser Technik auch Vertrauensintervalle fu ¨r einen Funktionswert an einer Stelle xo bestimmen. Dazu reparametrisieren wir das urspr¨ ungliche Problem so, dass ein Parameter, sagen wir φ1 , den Funktionswert h (xo ) repr¨asentiert und gehen dann wie besprochen vor.

¨ Die grafischen Hilfsmittel zur Uberpr¨ ufung der linearen Approximation beruhen auf dem gerade besprochenen t-Test, der ja eben diese N¨aherung nicht ben¨ utzt. Wir betrachten die Teststatistik Tk (5) als Funktion ihres Argumentes θk und nennen sie Profil-tFunktion (im letzten Kapitel war das Argument jeweils mit θk∗ bezeichnet, nun lassen wir das ∗ einfachheitshalber weg). F¨ ur die lineare Regression erh¨alt man, wie aus (4) ersichtlich, eine Gerade, w¨ahrend f¨ ur die nichtlineare Regression eine monoton steigende Funktion herauskommt. Den grafischen Vergleich von Tk (θk ) mit einer Geraden erm¨oglicht der so genannte Profil-t-Plot. Es ist u ¨ blich, auf der horizontalen Achse

17

θ1 190

210

θ1 230

20

40

60

80

100 0.99

4

4 0.99

0.0 0.80

0.80 0.0

0

0.80

-2 -4

-2 0.99

Niveau

0

T1 (θ1 )

0.80

Niveau

T1 (θ1 )

2 2

0.99

-6

-4 -4

-2

0

δ (θ1 )

2

4

0

10

δ (θ1 )

20

30

ur jeweils den ersten Parameter in den Beispielen Puromycin und Abbildung 8: Profil-t-Plot f¨ Sauerstoffverbrauch. Die gestrichelte Linien zeigen die verwendete lineare Approximation und die gepunktete Linie die Konstruktion des 99% Vertrauensintervalls mit Hilfe von T1 (θ1 ).

nicht θk , sondern die auf Grund der linearen Approximation bestimmte standardisierte Version θk − θbk  δk (θk ) :=  se (θbk )

zu verwenden. Die Vergleichsgerade wird dann die Diagonale“, also die Gerade mit ” Steigung 1 und Achsenabschnitt 0. Je st¨arker die Profil-t-Funktion gekr¨ ummt ist, desto st¨arker ist die Nichtlinearit¨at in einer Umgebung von θk . Folglich zeigt diese Darstellung, wie gut die lineare Approximation in einer Umgebung von θbk ist. (Die Umgebung, die f¨ ur die Statistik wichtig ist, ist etwa durch |δk (θk ) | ≤ 2.5 bestimmt.) In Abbildung 8 zeigt sich, dass im Beispiel Puromycin die Nichtlinearit¨at minim, im Beispiel Sauerstoffverbrauch dagegen gross ist. Aus den Darstellungen kann man die Vertrauensintervalle gem¨ass (5) ablesen. Der Bequemlichkeit halber sind auf der rechten vertikalen Achse die Wahrscheinlichkeiten P (Tk ≤ t) gem¨ass der t-Verteilung markiert. Im Beispiel des Sauerstoff-Verbrauchs ergibt sich ein Vertrauensintervall ohne obere Grenze!

18

• •• • •••

Gewicht

800

•• ••• •••• •• •

600

Gewicht

1000

•• •• •

400

••

••• •

200 0

1

2

3

4 Konzentration

Konzentration

Abbildung 9: Beispiel Kresse. Links: Darstellung der Daten. Rechts : Ein typischer Verlauf

der zu verwendenden Regressionsfunktion.

5 Prognose und Eichung Neben der Frage, welche Parameter f¨ ur ein Modell bez¨ uglich gegebener Daten plausibel sind, ist oft die Frage, in welchem Bereich k¨ unftige Beobachtungen zu liegen kommen, von zentralem Interesse. Wir gehen davon aus, dass der unbekannte Parameter θ aus vorgegebenen Daten mit der Kleinste Quadrate Methode gesch¨atzt wird. (Wenn Sie wollen, k¨onnen Sie diesen Datensatz Trainingsdatensatz nennen.) Was k¨onnen wir nun mit Hilfe unseres Modells u unftige Beobachtung Yo an einer vorgegebener ¨ ber eine zuk¨ Stelle xo aussagen? ⊲ Beispiel Kresse. Die Konzentration eines agrochemischen Stoffes in Erdproben kann durch das Wachstumsverhalten einer bestimmten Kressesorte (Nasturtium) studiert werden. Es wurden je 6 Messungen der Zielgr¨osse Y an 7 Erdproben mit vorgegebenen (oder mit gr¨osst m¨oglicher Pr¨azision gemessenen) Konzentrationen x gemacht. In jedem Fall tun wir so, als ob die x-Werte keinen Messfehler aufweisen. Die Zielgr¨osse ist das Gewicht der Kresse pro Fl¨acheneinheit nach 3 Wochen. Ein “logit-log” Modell wurde verwendet, um den Zusammenhang zwischen Konzentration und Gewicht zu beschreiben:  θ if x = 0 1 h (x, θ) = θ1  if x > 0. 1+exp(θ2 +θ3 ln(x))

(Die Daten und die Funktion h (·) sind in Abbildung 9 graphisch dargestellt.) Wir k¨onnen uns jetzt fragen, welche Gewichtswerte sind bei einer Konzentration von z.B. x0 = 3 zu erwarten?

19

Wir k¨onnen  den  Erwartungswert E (Yo ) = h (xo , θ) der Zielvariablen Y an der Stelle xo b durch h xo , θ sch¨atzen. Wollen wir zus¨atzlich einen Bereich angeben, wo eine zuk¨ unftige Beobachtung h¨ochst wahrscheinlich zu liegen kommt, dann m¨ ussen wir nicht nur die   Streuung der Sch¨atzung h xo , θb ber¨ ucksichtigen sondern auch den zuf¨alligen Fehler Eo . Analog zur linearen Regression l¨asst sich ein zumindest gen¨ ahertes (1 − α/2)Prognoseintervall durch   tn−p b h xo , θ ± σbv qα/2 (6) angeben. Dabei ist σb2v die Summe der Fehlervarianz und der Varianz der Funktionssch¨atzung bei xo : σb2v = σb2 + σb2xo . Die Berechnung von σb2xo finden Sie in (4).

Falls die Stichprobengr¨osse n im Trainingsdatensatz sehr gross ist, wird die gesch¨atzte Varianz σb2v durch die Fehlervarianz σb2 dominiert. Das bedeutet, dass die Unsicherheit in der Prognose dann vor allem durch die Beobachtungsfehler bestimmt wird. Der zweite Summand im Ausdruck von σb2v wiedergibt die Unsicherheit, die durch die Sch¨atzung von θ verursacht wird.

Es ist daher klar, dass das Prognose-Intervall breiter ist als das Vertrauensintervall f¨ ur den Erwartungswert, da ja noch die zuf¨allige Streuung der Beobachtung ber¨ ucksichtigt werden musste. Die Endpunkte solcher Intervalle f¨ ur xo -Werte in einem ausgew¨ahlten Bereich wurden in Abbildung 10 links miteinander verbunden und als Band dargestellt. Die Bestimmung des Prognoseintervalls in (6) beruht auf den gleichen N¨aherungen wie sie in Kapitel 3 verwendet wurden. Die Qualit¨at der N¨aherung kann wieder mit dem Profil-t-Plot untersucht werden. Das eigentliche Ziel des Experimentes im Beispiel Kresse ist, aus dem Gewicht der Kresse die Konzentration des agrochemischen Stoffes zu sch¨atzen. Folglich m¨ochte man die Regressionsbeziehung in der “falschen” Richtung verwenden. Daraus ergeben sich Probleme bei der Inferenz. Trotzdem ist ein solches Vorgehen h¨aufig erw¨ unscht, um eine Mess-Methode zu eichen oder um aus dem Resultat einer (billigen) Mess-Methode das Resultat einer anderen (teuren) zu sch¨atzen. Die Regressions-Kurve wird in diesem Zusammenhang oft Eich-Kurve genannt. – Die Behandlung dieses Problems findet man in Englisch unter den Stichw¨ortern inverse regression und calibration.

20

Wir wollen hier eine einfache Behandlung vorstellen, die ein brauchbares Resultat ergibt, wenn vereinfachende Annahmen zutreffen. Wir nehmen an, dass die x-Werte keinen Messfehler aufweisen. Das erreicht man, indem man im Beispiel die Konzentration des agrochemischen Stoffes sehr sorgf¨altig bestimmt. F¨ ur mehrere solche Erdproben mit m¨oglichst unterschiedlichen Konzentrationen f¨ uhrt man jeweils mehrere (m¨oglichst) unabh¨angige Messungen der Gr¨osse Y durch. Das ergibt dann einen Trainings-Datensatz, mit dem die unbekannten Parameter gesch¨atzt und die Standardfehler der Sch¨atzungen bestimmt werden k¨onnen. Wenn nun f¨ ur eine zu messende Probe der Wert yo abgelesen wird, ist es nahe liegend, den zugeh¨origer xo -Wert folgendermassen zu bestimmen:   xbo = h−1 yo , θb .

Dabei bezeichnet h−1 die Umkehrfunktion von h. Dieses Vorgehen ist allerdings nur dann korrekt, wenn h (·) monoton wachsend oder fallend ist. Diese Voraussetzung ist jedoch bei Eichproblemen meistens erf¨ ullt. Anschliessend an diese Rechnung stellt sich die Frage, wie genau dieser Wert eigentlich ist. Diese Problemstellung sieht ¨ahnlich aus wie das Prognoseproblem am Anfang dieses Kapitels. Im Gegensatz zum Prognose-Problem wird jedoch yo beobachtet und der entsprechende Wert xo muss gesch¨atzt werden. Die Antwort l¨asst sich formulieren, indem wir xo als Parameter ansehen, f¨ ur den ein Vertrauensintervall gesucht ist. Ein solches Intervall ergibt sich (wie immer) aus einem Test. Nehmen wir als Nullhypothese xo = xH an! Wie wir in (6) gesehen haben, liegt Y mit Wahrscheinlichkeit 0.95 im Prognose-Intervall q q   tn−p tn−p 2 2 b b T A) b −1 ao , h xo , θ ± bxo mit bxo = q0.975 σb + σbxo = q0.975 σb 1 + aTo (A

das in Abbildung 10 – gleich f¨ ur alle m¨oglichen xo in einem vorgegebenen Bereich – f¨ ur das Beispiel Kresse dargestellt ist. Das Intervall bildet deshalb ein Annahmeintervall f¨ ur die Gr¨osse Yo (die hier die Rolle einer Test-Statistik spielt) unter der Nullhypothese xo = xH . Die Abbildung 10 rechts veranschaulicht nun den weiteren Gedankengang am Beispiel Kresse: Messwerte yo sind mit Parameterwerten xo vereinbar im Sinne des Tests, wenn der Punkt [xo , yo ] zwischen den eingezeichneten Kurven liegt. In der Abbildung kann man deshalb ohne Schwierigkeiten die Menge der xo -Werte bestimmen, die mit der Beobachtung yo vertr¨aglich sind. Sie bilden das eingezeichnete Intervall, welches auch durch

21

Gewicht

800

•• •• ••• • •• • • ••

600

1000 800 •• • ••

Gewicht

1000

•• •

400

••

200

••• •

•• •• ••• • •• • • ••

600

•• • •• •• •

400

••

200

0

••• •

0 0

1

2

3

4

5

0

1

Konzentration

2

3

4

5

Konzentration

Abbildung 10: Beispiel Kresse. Links: Vertrauensband f¨ ur die gesch¨ atzte Regressionskurve

(gestrichelt) und Prognoseband (ausgezogen). Rechts: Schematische Darstellung, wie ein Kalibrationsintervall bestimmt wird, an den Stellen y0 = 650 und y0 = 350. Die resultierenden Intervalle sind [0.4, 1.22] respektive [1.73, 4.34].

die Menge

n   o p tn−p x : |y0 − h x, θb | ≤ σb2 + σb2x q1−α/2 .

beschrieben werden kann. Dieses Intervall ist nun das gesuchte Vertrauensintervall (oder Kalibrationsintervall) f¨ ur xo . Falls wir m Werte haben, um yo zu bestimmen, so wenden wir das obige Verfahren auf P y¯o = m j=0 yoj /m an: (



 b x : |¯ yo − h x, θ | ≤

r

σb2 m

t

n−p + σb2x q1−α/2

)

.

6 Schlussbemerkungen ⊲ Beispiel Sauerstoffverbrauch. Warum haben wir so viele Schwierigkeiten mit dem Beispiel Sauerstoffverbrauch? Betrachten wir die Abbildung 2 und erinnern uns, dass der Parameter θ1 den erwarteten Sauerstoffverbrauch bei unendlicher Inkubinationszeit repr¨asentiert, so ist klar, dass θ1 schwierig zu sch¨atzen ist, weil die horizontale Asymptote durch die Daten schlecht bestimmt ist. H¨atten wir noch weitere Beobachtungen mit l¨angeren Inkubinationszeiten, so h¨atten wir die Schwierigkeiten mit der Qualit¨at des Vertrauensintervalles von θ vermeiden k¨onnen.

22

Gerade auch bei nichtlinearen Modellen ist eine gute (statistische) Versuchsplanung (experimental design) unerl¨asslich. Der Informationsgehalt der Daten wird durch die Wahl der Versuchsbedingungen festgelegt, und kein (statistisches) Verfahren ist in der Lage, Informationen u ¨ ber das Modell zu liefern, welche in den Daten nicht enthalten sind. Neben der Profil-t-Funktion sind die Likelihood-Profilspuren (likelihood profile tra¨ ces) weitere n¨ utzliche Hilfsmittel zurUberpr¨ ufung der Approximation. Eine geeignete Darstellung der Likelihood-Profilspuren zeigt jedoch nicht nur Nichtlinearit¨aten, sie enth¨alt auch n¨ utzliche Hinweise, wie sich die Parameter gegenseitig beeinflussen, und welche Effekte diesbez¨ uglich Parametertransformationen haben. Mehr dar¨ uber finden Sie in meinen Unterlagen f¨ ur den Nachdiplomkurs in angewandter Statistik [3] oder im Buch von Bates und Watts [1]. Diese Unterlagen beruhen auf Unterlagen f¨ ur den Nachdiplomkurs in angewandter Statistik [3] und dem Buch von Bates und Watts [1]. Eine mathematischere Diskussion u ¨ber die statistischen und numerischen Methoden in der nichtlinearen Regression findet sich in Seber und Wild [4]. Das Buch von Ratkowsky [5] z¨ahlt zahlreiche m¨ogliche nichtlineare Funktionen h (·) auf, die vor allem in biologischen Bereichen ihre Anwendung finden. Seit einiger Zeit wird zur Bestimmung von Vertrauensbereichen auch der Bootstrap benutzt (siehe [6]). In diesem Buch wird auch der Fall mit nichtkonstanter Varianz (heteroskedastic models) besprochen. Dazu lohnt auch ein Blick in das Buch von Carrol and Ruppert [7]. Heutzutage enthalten die meisten Statistik-Pakete eine Prozedur, welche nichtlineare Modelle anpassen und asymptotische Vertrauensintervalle f¨ ur die Parameter berechnen kann. Prinzipiell ist es dann auch m¨oglich, t-Profile“ (und Likelihood-Profilspuren) zu ” berechnen, weil sie auch auf dem Anpassen von nichtlinearen Modellen, jedoch mit einem Parameter weniger, beruhen. In den beiden Implementationen S-Plus und R der Statistik-Sprache S sind die Funktionen nls und profile verf¨ ugbar, die auf der Arbeit von Bates und Watts [1] basieren. Zusammenfassungen der Methode finden sich deshalb in Kapitel 10 von [8], in den SPlus-Manuals oder in [9]. Die Library“ nlme enth¨alt S-Funktionen, die nichtlineare ” Regressionsmodelle mit korrelierten Fehlern (gnls) und zuf¨alligen Effekten (nlme) an Daten anpassen k¨onnen. Diese Implementationen basieren auf dem Buch “Mixed Effects Models in S and S-Plus” von [10].

23

Literatur [1] D. M. Bates and D. G. Watts, Nonlinear Regression Analysis & Its Applications. John Wiley & Sons, 1988. [2] I. Rapold-Nydegger, Untersuchungen zum Diffusionsverhalten von Anionen in carboxylierten Cellulosemembranen. PhD thesis, ETH Zurich, 1994. [3] A. Ruckstuhl, “Einf¨ uhrung in die nichtlineare Regression.” Unterlagen f¨ ur den NDK in angewandter Statistik der ETH Z¨ urich, 2006. [4] G. Seber and C. Wild, Nonlinear regression. Wiley, New York, 1989. [5] D. A. Ratkowsky, Handbook of Nonlinear Regression Models. Marcel Dekker, New York, 1989. [6] S. Huet, A. Bouvier, M.-A. Gruet, and E. Jolivet, Statistical Tools for Nonlinear Regression: A Practical Guide with S-Plus Examples. Springer-Verlag, New York, 1996. [7] R. Carroll and D. Ruppert, Transformation and Weighting in Regression. Wiley, New York, 1988. [8] J. M. Chambers and T. J. Hastie, Statistical Models in S. Brooks/Cole, Pacific Grove, California, 1992.

Wadsworth &

[9] W. N. Venables and B. Ripley, Modern Applied Statistics with S. Statistics and Computing, Springer-Verlag, New York, fourth ed., 2002. [10] J. C. Pinheiro and D. M. Bates, Mixed-Effects Models in S and S-PLUS. Statistics and Computing, Springer, 2000. Autorenangabe Prof. Dr. Andreas Ruckstuhl Institut f¨ ur Datenanalyse und Prozessdesign Z¨ urcher Hochschule Winterthur J¨agerstrasse 2 Postfach 805 CH-8401 Winterthur Tel.: +41 (0)52 267 78 12 Fax: +41 (0)52 268 78 12 E-mail: [email protected] 24