48  Lösungen Nicht-Lineare Regression

Hier finden Sie die Lösungen zu den Übungsaufgaben von Abschnitt 44.4.

Die hier vorgestellten Lösungen stellen immer nur eine mögliche Vorgehensweisen dar und sind sicherlich nicht der Weisheit letzter Schluss. In R führen viele Wege nach Rom, und wenn Sie mit anderem Code zu den richtigen Ergebnissen kommen, dann ist das völlig in Ordnung.

48.1 Lösung zur Aufgabe 44.4.1

a) Erstellen Sie ein Datenframe mit den Variablen Stunden und Bakterien.
# Daten erzeugen
df <- data.frame(Stunden=c(1:9),
                 Bakterien=c(25, 28, 47, 65, 86, 121, 190, 290, 362))
b) Erzeugen Sie ein Scatterplot. Welche Regression würden Sie auf Grundlage des Plots vorschlagen?
# plot()
plot(df$Stunden, df$Bakterien)
# ggplot()
ggplot(df, aes(x=Stunden, y=Bakterien)) +
  geom_point()

Die Punktwolken sprechen für einen exponentiellen Anstieg.

c) Berechnen Sie die quadratischen und exponentiellen Modelle für die Bakterienvermehrung über die Zeit.
# quadratisch
q <- lm(Bakterien ~ Stunden + I(Stunden^2), data=df)
summary(q)

Call:
lm(formula = Bakterien ~ Stunden + I(Stunden^2), data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-16.617  -8.297  -1.430   9.916  15.442 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   53.2381    15.9862   3.330   0.0158 *  
Stunden      -26.7420     7.3403  -3.643   0.0108 *  
I(Stunden^2)   6.8009     0.7159   9.500 7.75e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.56 on 6 degrees of freedom
Multiple R-squared:  0.9919,    Adjusted R-squared:  0.9892 
F-statistic: 368.8 on 2 and 6 DF,  p-value: 5.254e-07
# exponentiell
e <- lm(log(Bakterien) ~ Stunden, data=df)
summary(e)

Call:
lm(formula = log(Bakterien) ~ Stunden, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.12676 -0.06057  0.01145  0.03920  0.11190 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.75498    0.06174   44.62 7.41e-10 ***
Stunden      0.35199    0.01097   32.08 7.39e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.08498 on 7 degrees of freedom
Multiple R-squared:  0.9932,    Adjusted R-squared:  0.9923 
F-statistic:  1029 on 1 and 7 DF,  p-value: 7.389e-09
d) Plotten Sie das bessere Modell in die Punktwolke.
# Vorbereitung quadratisch
vorhersageQ <- predict(q, list(Stunden=df$Stunden))
# plot() quadratisch
plot(df$Stunden, df$Bakterien, 
     main="quadratischer Zusammenhang")
lines(df$Stunden, vorhersageQ, col="tan4")
# ggplot() quadratisch
ggplot(df, aes(x=Stunden, y=Bakterien)) + 
  geom_point() +
  geom_line(aes(y=vorhersageQ), col="tan4")
# Vorbereitung exponentiell
vorhersageE <- exp(predict(e, list(Stunden=df$Stunden)))

# plot() exponentiell
plot(df$Stunden, df$Bakterien, 
     main="exponentieller Zusammenhang")
lines(df$Stunden, vorhersageE, col="wheat")
# ggplot() exponentiell
ggplot(df, aes(x=Stunden, y=Bakterien)) + 
  geom_point() +
  geom_line(aes(y=vorhersageE), col="wheat")

e) Wie viele Bakterien werden nach dem besten Modell 3 Stunden nach Anlegen der Kultur vorhanden sein? Und nach 10 Stunden? Sind diese Vorhersagen zuverlässig?
# Vorhersage
exp(predict(e, list(Stunden=c(3, 10))))
        1         2 
 45.19322 531.05241 

Nach 3 Stunden können wir 46 Bakterien erwarten, nach 10 Stunden 532.

f) Machen Sie eine möglichst zuverlässige Vorhersage über die Zeit, die benötigt wird, um 100 Bakterien in der Kultur zu haben.
# neues Modell
df$BakterienLog <- log(df$Bakterien)
fit <- lm(BakterienLog ~ Stunden, data=df)
a <- exp(fit$coefficients[1])
b <- fit$coefficients[2]
# Vorhersage für 100 Bakterien
(log(100) - log(a)) / b
(Intercept) 
   5.256395 

Nach ca 5.3 Stunden sind 100 Bakterien zu erwarten.

48.2 Lösung zur Aufgabe 44.4.2

a) Laden Sie den Datensatz diet in Ihre R-Session.
# lade Datensatz
load(url("https://www.produnis.de/R/data/diet.RData"))
b) Erstellen Sie eine Punktwolke. Welche Art von Modell erklärt auf Grundlage der Punktwolke den Gewichtsverlust pro Diättag besser?
# plot()
plot(diet$days, diet$weight.loss)
# ggplot()
ggplot(diet, aes(x=days, y=weight.loss)) +
  geom_point()

c) Berechnen Sie das Regressionsmodell, welches den Gewichtsverlust mit der Anzahl an Diättagen am besten (im Vergleich zu anderen) erklären kann. Wird das Modell zuverlässige Vorhersagen machen?
# Vergleiche Bestimmheitsmaße verschiedener Modelle

# quadratisches Modell
q <- lm(weight.loss ~ days + I(days^2), data=diet)

# exponentielles Modell
e <- lm (log(weight.loss) ~ days, data=diet)

# logarithmisches Modell
l <- lm(weight.loss ~ log(days), data=diet)

# sigmoidales Modell
s <- lm(log(weight.loss) ~ I(1/days), data=diet)

result <- data.frame(Modell = c("quadratisch", "exponentiell",
                                "logarithmisch", "sigmoidal"),
                     R.square = c(summary(q)$r.square,
                                 summary(e)$r.square,
                                 summary(l)$r.square,
                                 summary(s)$r.square))
# Anzeigen
result[order(result$R.square, decreasing = TRUE),]
         Modell  R.square
4     sigmoidal 0.6662170
1   quadratisch 0.5397848
3 logarithmisch 0.5254856
2  exponentiell 0.4308936
## oder mit der compare.lm()-Funktion
jgsbook::compare.lm(diet$weight.loss, diet$days)
         Modell  R.square
6     sigmoidal 0.6662170
7        potenz 0.5684490
3       kubisch 0.5584355
2   quadratisch 0.5397848
5 logarithmisch 0.5254856
1        linear 0.4356390
4  exponentiell 0.4308936

Das sigmoidale Modell kann die Daten am besten erklären, da in diesem Modell das Bestimmtheitsmaß am größten ist.

d) Plotten Sie Ihr Modell.
# sigmoidales Modell
s <- lm(log(weight.loss) ~ I(1/days), data=diet)

# vorhersage vorbereiten
tage <- seq(min(diet$days), max(diet$days))
# alle "tage" vorhersagen
vorhersage <- predict(s, list(days=tage))
vorhersage[-1]=exp(vorhersage[-1])

# plot()
plot(diet$days, diet$weight.loss)
lines(tage, vorhersage, col="purple")

# ggplot()
# in Datenframe für ggplot speichern
helper <- data.frame(tage, vorhersage)

ggplot(diet, aes(x=days, y=weight.loss)) + 
  geom_point(color="#FF9999", size=3) +
    # Legend
    scale_linetype("Regression model") +
    # Sigmoidal model
    geom_line(data=helper, aes(x=tage, y=vorhersage, linetype="Sigmoidal"))

Auch die Vorhersagewerte für die Idealkurve können mittels compare.lm() und dem Parameter predict=TRUE erzeugt werden.

help <- jgsbook::compare.lm(diet$weight.loss, diet$days, 
                            predict=TRUE)

# anschauen                            
head(help)
  pred.x     line     quad     cube     expo     loga     sigm    power
1  14.00 7.645786 4.766305 3.526750 6.714088 5.166988 1.262199 4.966500
2  14.01 7.647113 4.770032 3.532726 6.715040 5.171506 3.537806 4.969006
3  14.02 7.648440 4.773757 3.538700 6.715992 5.176020 3.542430 4.971512
4  14.03 7.649767 4.777483 3.544671 6.716944 5.180531 3.547052 4.974018
5  14.04 7.651094 4.781207 3.550641 6.717896 5.185040 3.551675 4.976523
6  14.05 7.652422 4.784931 3.556608 6.718848 5.189544 3.556296 4.979027
  logistic
1 4.261378
2 4.264669
3 4.267962
4 4.271256
5 4.274552
6 4.277850
# plot()
plot(diet$days, diet$weight.loss)
lines(help$pred.x, help$sigm, col="purple")

# ggplot()
ggplot(diet, aes(x=days, y=weight.loss)) + 
  geom_point(color="#FF9999", size=3) +
    # Legend
    scale_linetype("Regression model") +
    # Sigmoidal model
    geom_line(data=help, aes(x=pred.x, y=sigm, linetype="Sigmoidal"))

e) Berechnen Sie das Regressionsmodell, das den Gewichtsverlust anhand der Tage der Diät für die Gruppe der Personen, die sich nicht regelmäßig körperlich betätigen, am besten erklärt.
# Subset bilden
df1 <- subset(diet, exercise=="no")

# Vergleiche Bestimmheitsmaße verschiedener Modelle

# quadratisches Modell
q1 <- lm(weight.loss ~ days + I(days^2), data=df1)

# exponentielles Modell
e1 <- lm (log(weight.loss) ~ days, data=df1)

# logarithmisches Modell
l1 <- lm(weight.loss ~ log(days), data=df1)

# sigmoidales Modell
s1 <- lm(log(weight.loss) ~ I(1/days), data=df1)

result1 <- data.frame(Modell = c("quadratisch", "exponentiell",
                                 "logarithmisch", "sigmoidal"),
                     R.square = c(summary(q1)$r.square,
                                  summary(e1)$r.square,
                                  summary(l1)$r.square,
                                  summary(s1)$r.square))
result1[order(result1$R.square, decreasing = TRUE),]
         Modell  R.square
4     sigmoidal 0.7401212
1   quadratisch 0.7100610
3 logarithmisch 0.6494521
2  exponentiell 0.5222832
# oder mittels compare.lm()
jgsbook::compare.lm(df1$weight.loss, df1$days)
         Modell  R.square
6     sigmoidal 0.7401212
3       kubisch 0.7151929
2   quadratisch 0.7100610
7        potenz 0.6700051
5 logarithmisch 0.6494521
1        linear 0.5286338
4  exponentiell 0.5222832

Das sigmoidale Modell liefert wieder die beste Erklärung der Daten.

f) Wiederholen Sie die Analyse für die Gruppe, die sich regelmäßig körperlich betätigt.
# Subset bilden
df2 <- subset(diet, exercise=="yes")

# Vergleiche Bestimmheitsmaße verschiedener Modelle

# quadratisches Modell
q2 <- lm(weight.loss ~ days + I(days^2), data=df2)

# exponentielles Modell
e2 <- lm (log(weight.loss) ~ days, data=df2)

# logarithmisches Modell
l2 <- lm(weight.loss ~ log(days), data=df2)

# sigmoidales Modell
s2 <- lm(log(weight.loss) ~ I(1/days), data=df2)

result2 <- data.frame(Modell = c("quadratisch", "exponentiell",
                                 "logarithmisch", "sigmoidal"),
                     R.square = c(summary(q2)$r.square,
                                  summary(e2)$r.square,
                                  summary(l2)$r.square,
                                  summary(s2)$r.square))
result2[order(result1$R.square, decreasing = TRUE),]
         Modell  R.square
4     sigmoidal 0.8305013
1   quadratisch 0.7791671
3 logarithmisch 0.7885173
2  exponentiell 0.4945564
# oder mittels compare.lm()
jgsbook::compare.lm(df2$weight.loss, df2$days)
         Modell  R.square
3       kubisch 0.8326179
6     sigmoidal 0.8305013
5 logarithmisch 0.7885173
2   quadratisch 0.7791671
7        potenz 0.6704843
1        linear 0.6623502
4  exponentiell 0.4945564

Das kubische Modell liefert hier die beste Erklärung der Daten.

g) Benutzen Sie die erstellen Modelle, um den Gewichtsverlust nach 30 und nach 100 Tagen Diät für Personen, die sich körperlich betätigen, und für solche, die dies nicht tun, vorherzusagen. Sind diese Vorhersagen zuverlässig?
# Vorhersage Kein Sport
exp(predict(s1, list(days=c(30,100))))
        1         2 
 7.808926 13.806339 
# Vorhersage Sport
exp(predict(s2, list(days=c(30,100))))
       1        2 
11.28578 21.13143 

48.3 Lösung zur Aufgabe 44.4.3

df <- data.frame(Stunden = c(2:8),
                 Konzentration = c(25, 36, 48, 64, 86, 114, 168))
a) Benutzen Sie ein exponentielles Modell, um die Konzentration nach 10 Stunden vorherzusagen. Ist die Vorhersage zuverlässig?
# exponentielles Modell
fit <- lm(log(Konzentration) ~ Stunden, data=df)
summary(fit)

Call:
lm(formula = log(Konzentration) ~ Stunden, data = df)

Residuals:
        1         2         3         4         5         6         7 
-0.023147  0.034218  0.014623 -0.004972 -0.016786 -0.042212  0.038276 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.627468   0.033673   78.03 6.55e-09 ***
Stunden     0.307277   0.006253   49.14 6.59e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03309 on 5 degrees of freedom
Multiple R-squared:  0.9979,    Adjusted R-squared:  0.9975 
F-statistic:  2415 on 1 and 5 DF,  p-value: 6.594e-08
# Vorhersage 10 Stunden
exp(predict(fit, list(Stunden=10)))
     1 
298.94 

Nach 10 Stunden beträgt die Konzentration 298,94 mg/dl.

Das Bestimmtheitsmaß R2 des Modells ist mit 0,9979 sehr groß. Die Daten werden sehr gut durch das Modell erklärt.

b) Benutzen Sie ein logarithmisches Modell um zu bestimmen, nach wie vielen Stunden eine Konzentration von 100 mg/dl erreicht sein wird.
# exponentielles Modell
fit <- lm(log(Konzentration) ~ Stunden, data=df)

# coefficienten
a <- fit$coefficient[1]
b <- fit$coefficient[2]

# Vorhersage 100 mg/dl
( log(100) - a ) / b
(Intercept) 
   6.436209 

Nach 6,436209 Stunden sind 100 mg/dl erreicht.