08-Nichtlineare Regression

Applied Statistics – A Practical Course

Thomas Petzoldt

2025-10-15

Achtung: deutsche Übersetzung muss noch geprüft werden

Nichtlineare Regression

Viele Phänomene in der Natur sind nichtlinear!

Linear oder nichtlinear?

  • Einige nichtlineare Probleme können mit linearen Methoden gelöst werden
  • z.B. Polynome oder periodische (Sinus- und Cosinus-) Funktionen
  • andere können durch Transformationen angepasst werden

Beispiel

\[y = a \cdot x^b\] kann umgewandelt werden in:

\[\ln(y) = \ln(a) + b \cdot \ln(x)\]

  • aber: Linearisierung transformiert auch die Residuen
  • die Transformation ist manchmal richtig und manchmal falsch
  • Homogenität der Varianzen

Linearisierung oder direkte nichtlineare Anpassung?


Linearisierende Transformationen

  • Logarithmus, Quadratwurzel, Kehrwerte …
  • können angewandt werden, wenn die Restvarianz homogen bleibt (oder wird).
  • aber in einigen Fällen führen die Transformationen zu Heteroskedastizität
    \(\Rightarrow\) verfälschte Regressionsparameter.

Nichtlineare Regression

  • sehr flexible, benutzerdefinierte Funktionen
  • keine Transformation erforderlich
  • aber: erfordert iterative Optimierungsmethoden
  • theoretisch können nur lokale Optima gefunden werden
  • Parameter sind nicht in allen Fällen identifizierbar

Nichtlineare Regression


Iterative Suche nach dem Minimum der Summe der Quadrate

Wie lässt sich das globale Minimum finden?


Anpassungsgüte: Residuensumme der quadrierten Differenzen \(RSS\):

\[ RSS = \sum_{i=1}^n\left(y_i- f(\mathbf x_i, \mathbf p)\right)^2 = \min ! \]

mit \(y\): abhängige Variable, \(\mathbf x\): Matrix der unabhängigen Variablen,\(\mathbf p\): Parametervektor, \(n\): Stichprobenumfang

Verwendung eines iterativen Lösers \(\rightarrow\) siehe nächste Folien

Nichtlinearer Bestimmungskoeffizient \(R^2\)

  • steht in keinem Zusammenhang mit dem (linearen) Korrelationskoeffizienten
  • kann aus den Rest- und Gesamtvarianzen berechnet werden

\[ R^2 = 1 - \frac{\text{Varianz der Residuen}}{\text{Varianz der y-Daten}} = 1- \frac{s^2_\varepsilon}{s^2_y} \]

  • nichtlinear \(r^2\) misst den Anteil der erklärten Varianz
  • … andere Indizes können zusätzlich verwendet werden, z.B. MSE, RMSE, NSE, …

Allgemeines Prinzip von Optimierungsalgorithmen


Das Minimum der quadrierten Residuen (RSS) wird durch Iteration gesucht.:


  1. Beginne mit einer Anfangsschätzung für den Parametersatz.
  2. Versuche, einen Parametersatz mit niedrigerer RSS zu finden.
  3. Ist der neue Parametersatz besser?
    • Nein: Verwerfe die neuen Parameter und gehe zu 2
    • Ja: Akzeptiere die neuen Parameter und springe zu 4
  4. Ist der neue Parametersatz gut genug?
    • Nein: Springe zu 2
    • Ja: Springe zu 5
  5. Parametersatz gefunden.

Deterministische Methoden


Gradient Descent

  • einen Schritt in die Richtung des steilsten Abstiegs gehen
  • \(\rightarrow\) relativ einfach
  • \(\rightarrow\) relativ robust für „kompliziertere“ Probleme

Newton-Methode

  • quadratische Approximation der RSS-Funktion
  • Versuche, das Minimum zu schätzen
  • \(\rightarrow\) berücksichtigt die Krümmung
  • \(\rightarrow\) ist schneller für „gut funktionierende“ Probleme
  • mehrere Versionen:
    quasi-Newton, Gauss-Newton, L-BFGS, …

Levenberg-Marquardt

  • interpoliert zwischen Gauß-Newton und Gradient Descent.

Die Newton-Methode (rot) nutzt Krümmungsinformationen, um schneller zu konvergieren als der Gradient Descent (grün).
Quelle: Wikipedia.

Stochastische Methoden


Anwendung statistischer Prinzipien (Zufallssuche)

Klassische Methoden

  • Simuliertes Glühen: simuliert die Erwärmung und Abkühlung von Materie \(\rightarrow\) „Kristallisationsprozess“.
  • Kontrollierte Zufallssuche Price (1983)

Evolutionäre Algorithmen

  • Analogie zur Genetik: Mutation und Selektion
  • verfolgt eine „Population“ von mehreren Parametersätzen parallel
  • Informationsaustausch („Crossover“) zwischen Parametersätzen
  • \(\rightarrow\) für komplizierte Probleme mit einer großen Anzahl von Parametern
  • heutzutage in Microsoft Excel und LibreOffice Calc eingebaut

… und vieles mehr.

Beispiele

  • Enzymkinetik
  • Wachstum von Organismen
  • Kalibrierung von komplexen Modellen
    • in Chemie, Biologie, Ingenieurwesen, Finanzwirtschaft und Sozialwissenschaften
    • Wasser: Hydrologie, Hydrophysik, Grundwasser, Abwasser, Wasserqualität …

Enzymkinetik


… kann mit der bekannten Michaelis-Menten-Funktion beschrieben werden:

\[ v = v_{max} \frac{S}{k_m + S} \]

Linearisierung vs. (echte) nichtlineare Regression


Linearisierende Transformation

[>] Angemessen, wenn die Transformation die Homogenität der Varianzen verbessert.
[+] Schnell, einfach und leicht.
[+] Analytische Lösung liefert das globale Optimum.
[-] Nur eine begrenzte Anzahl von Funktionen kann angepasst werden.
[-] Kann zu einer falsch transformierten Fehlerstruktur und verzerrten Ergebnissen führen.

Nichtlineare Regression

[>] Geeignet, wenn die Fehlerstruktur bereits homogen ist und/oder keine analytische Lösung existiert.
[+] Kann zur Anpassung beliebiger Funktionen verwendet werden, sofern die Parameter identifizierbar sind.
[-] Benötigt Startwerte und beträchtliche Berechnungszeit.
[-] Beste Lösung (globales Optimum) ist nicht garantiert.

Nichtlineare Regression in R: einfache exponentielle Regression

Modell anpassen

# Beispieldaten
x <- 1:10
y <- c(1.6, 1.8, 2.1, 2.8, 3.5, 
       4.1, 5.1, 5.8, 7.1, 9.0)

# Anfangsparameter für den Optimierer
pstart <- c(a = 1, b = 1)

# nichtlineare kleinste Quadrate
fit <- nls(y ~ a * exp(b * x), start = pstart)
summary(fit)

Formula: y ~ a * exp(b * x)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a 1.263586   0.049902   25.32 6.34e-09 ***
b 0.194659   0.004716   41.27 1.31e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1525 on 8 degrees of freedom

Number of iterations to convergence: 13 
Achieved convergence tolerance: 5.956e-08

Plotte Ergebnisse

# zusätzliche x-Werte, für eine geglättete Kurve
x1 <- seq(1, 10, 0.1)
y1 <- predict(fit, data.frame(x = x1))
plot(x, y)
lines(x1, y1, col = "red")

Angepasste Parameter



Formula: y ~ a * exp(b * x)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a 1.263586   0.049902   25.32 6.34e-09 ***
b 0.194659   0.004716   41.27 1.31e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1525 on 8 degrees of freedom

Number of iterations to convergence: 13 
Achieved convergence tolerance: 5.956e-08
  • Estimate: die angepassten Parameter
  • Std. Error: \(s_{\bar{x}}\): zeigt die Zuverlässigkeit der Schätzung an
  • t- und p-Werte: keine Überinterpretation!
  • In der nichtlinearen Welt können „nicht-signifikante“ Parameter strukturell notwendig sein.

Bestimmungskoeffizient \(r^2 = 1-\frac{s^2_\varepsilon}{s^2_y}\)

(Rsquared <- 1 - var(residuals(fit))/var(y))
[1] 0.9965644

Michaelis-Menten-Kinetik


Code

S <-c(25, 25, 10, 10, 5, 5, 2.5, 2.5, 1.25, 1.25)
V <-c(0.0998, 0.0948, 0.076, 0.0724, 0.0557,
      0.0575, 0.0399, 0.0381, 0.017, 0.0253)

## Michaelis-Menten-Kinetik
f <- function(S, Vm, K) { Vm * S/(K + S) }

pstart <- c(Vm = 0.1, K = 5)
model_fit <- nls(V ~ f(S, Vm, K), start = pstart)
summary(model_fit)

plot(S, V, xlim = c(0, max(S)), ylim = c(0, max(V)))
x1 <-seq(0, 25, length = 100)
lines(x1, predict(model_fit, data.frame(S = x1)), col = "red")


Bestimmungskoeffizient

(1 - var(residuals(model_fit))/var(V))
[1] 0.9895147

Ergebnisse


Formula: V ~ f(S, Vm, K)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
Vm  0.11713    0.00381   30.74 1.36e-09 ***
K   5.38277    0.46780   11.51 2.95e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.003053 on 8 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 6.678e-06

Michaelis-Menten-Kinetik mit Selbststart

  • Funktion SSmicmen bestimmt Startparameter automatisch.
  • Nur wenige Selbststartfunktionen in R verfügbar.

Code

S <- c(25, 25, 10, 10, 5, 5, 2.5, 2.5, 1.25, 1.25)
V <- c(0.0998, 0.0948, 0.076, 0.0724, 0.0557,
       0.0575, 0.0399, 0.0381, 0.017, 0.0253)

model_fit <- nls(V ~ SSmicmen(S, Vm, K))

plot(S, V, xlim = c(0, max(S)), ylim = c(0, max(V)))
x1 <-seq(0, 25, length = 100)
lines(x1, predict(model_fit, data.frame(S = x1)), col="red")

Ergebnisse

summary(model_fit, correlation=TRUE)

Formula: V ~ f(S, Vm, K)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
Vm  0.11713    0.00381   30.74 1.36e-09 ***
K   5.38277    0.46780   11.51 2.95e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.003053 on 8 degrees of freedom

Correlation of Parameter Estimates:
  Vm  
K 0.88

Number of iterations to convergence: 3 
Achieved convergence tolerance: 6.678e-06

Plot

Anmerkung: Korrelation der Parameter

  • Hohe absolute Korrelationswerte deuten auf die Nichtidentifizierbarkeit von Parametern hin.
  • kritischer Wert hängt von den Daten ab
  • manchmal können bessere Startwerte oder ein anderer Optimierungsalgorithmus helfen

Praktische Hinweise


  1. Daten plotten
  2. Finde gute Ausgangswerte durch Nachdenken oder durch Trial and Error
  3. Vermeide sehr kleine und/oder sehr große Zahlen \(\longrightarrow\) skaliere das Problem auf Werte zwischen etwa 0,001 und 1000 um
  4. Beginne mit einer einfachen Funktion und füge nach und nach Terme und Parameter hinzu
  5. Nimm die Signifikanz von Parametern nicht zu ernst. Ein nicht signifikanter Parameter kann für die Struktur des Modells notwendig sein, sein Wegfall macht das gesamte Modell ungültig.

Weiterführende Literatur


  • Paket growthrates für Wachstumskurven: https://cran.r-project.org/package=growthrates

  • Paket FME für komplexere Modellanpassungsaufgaben (Identifizierbarkeitsanalyse, eingeschränkte Optimierung, mehrere abhängige Variablen und MCMC): (Soetaert & Petzoldt, 2010), https://cran.r-project.org/package=FME

  • Mehr über Optimierung in R: https://cran.r-project.org/web/views/Optimization.html

Anhang

Lineweaver-Burk-Transformation vs. nichtlineare Anpassung


Code der Methodenvergleichsfolie

S <-c(25, 25, 10, 10, 5, 5, 2.5, 2.5, 1.25, 1.25)
V <-c(0.0998, 0.0948, 0.076, 0.0724, 0.0557, 
      0.0575, 0.0399, 0.0381, 0.017, 0.0253)
model_fit<-nls(V ~ SSmicmen(S, Vm, K))

par(mfrow=c(1,2), las=1, lwd=2)
## Lineweaver-Burk
x <- 1/S; y <- 1/V

plot(x, y, xlab="1/S", ylab="1/V", xlim=c(-0.2,1), ylim=c(0, 80), pch=16,
  main="Linearisation\n(Lineweaver-Burk)")
abline(h=0, lwd=1, col="grey")
abline(v=0, lwd=1, col="grey")
m <- lm(y ~ x)
abline(m, col = "red")
Vm_l <- 1/coef(m)[1]
Km_l <- coef(m)[2] * Vm_l
#legend("topleft", c("vmax = 1/intercept", "km = slope * vmax"), 
#        box.lwd=1, bg="white")
text(0.35, 75, "1/V = 1/vmax + km / vmax * S")

## retransformed, both together
plot(S, V, xlim = c(0, max(S)),ylim=c(0, max(V)), pch=16, main="No Transformation")
x1 <-seq(0, 25, length=100)
lines(x1, Vm_l * x1 / (Km_l + x1), col="red")
lines(x1, predict(model_fit, data.frame(S=x1)), col="blue")
legend("bottomright", legend=c("linearisation", "nonlinear"), 
       box.lwd=1, lwd=2, col=c("red", "blue"))
text(15, 0.04, "V = S * vmax / (km + S)")

Siehe: https://en.wikipedia.org/wiki/Lineweaver-Burk_plot

Referenzen


Price, W. L. (1977). A controlled random search procedure for global optimization. The Computer Journal, 20(4), 367–370.
Price, W. L. (1983). Global optimization by controlled random search. Journal of Optimization Theory and Applications, 40(3), 333–348.
Soetaert, K., & Petzoldt, T. (2010). Inverse modelling, sensitivity and monte carlo analysis in R using package FME. Journal of Statistical Software, 33(3), 1–28. https://doi.org/10.18637/jss.v033.i03