Mathematik für Biologiestudierende II¶

Sommersemester 2024

18.06.2024

© 2024 Prof. Dr. Rüdiger W. Braun

In [1]:
import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.formula.api as smf
import seaborn as sns
sns.set_theme()

Lineare Modelle¶

Themen heute¶

  • Vorhersagen
  • Konfidenzintervalle für den Mittelwert und für die Einzelbeobachtung
  • mehr als eine erklärende Variable
  • irrelevante erklärende Variable erkennen

Vorhersagen (prediction)¶

In [2]:
df = pd.read_csv('galton.csv')
df.head()
Out[2]:
family father mother midparentHeight children childNum gender childHeight
0 001 78.5 67.0 75.43 4 1 male 73.2
1 002 75.5 66.5 73.66 4 1 male 73.5
2 002 75.5 66.5 73.66 4 2 male 72.5
3 003 75.0 64.0 72.06 2 1 male 71.0
4 004 75.0 64.0 72.06 5 1 male 70.5

Aufbereitung eines Datensatzes von Galton. Die Aufbereitung stammt aus den Begleitdaten zum Buch "Linear Models with Python" von Faraway

In [3]:
formel = 'childHeight ~ father'
  • In dieser Formel ist father die erklärende und Childheight die erkärte Variable
  • Die erklärende Variable heißt auch exogen, die erklärte endogen
In [4]:
modell = smf.ols(formel, df)
res = modell.fit()
res.summary()
Out[4]:
OLS Regression Results
Dep. Variable: childHeight R-squared: 0.154
Model: OLS Adj. R-squared: 0.152
Method: Least Squares F-statistic: 87.17
Date: Mon, 24 Jun 2024 Prob (F-statistic): 3.74e-19
Time: 21:06:29 Log-Likelihood: -1105.8
No. Observations: 481 AIC: 2216.
Df Residuals: 479 BIC: 2224.
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 38.3626 3.308 11.596 0.000 31.862 44.863
father 0.4465 0.048 9.337 0.000 0.353 0.540
Omnibus: 8.610 Durbin-Watson: 1.468
Prob(Omnibus): 0.014 Jarque-Bera (JB): 12.731
Skew: -0.110 Prob(JB): 0.00172
Kurtosis: 3.766 Cond. No. 2.08e+03


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.08e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
In [5]:
sns.regplot(df, x='father', y='childHeight');
No description has been provided for this image
  • Aufgabe: Wie groß ist im Mittel der Sohn, wenn der Vater 69.8 Zoll groß ist?

Erste Lösung:

  • Wir lesen ab
In [6]:
m = 0.4465
b = 38.36
In [7]:
x = 69.8
y = m*x + b
y
Out[7]:
69.5257

Zweite Lösung:

  • Wir verwenden die Methode get_prediction
  • Dazu müssen die Daten der erklärenden Variablen in einen DataFrame geschrieben werden
In [8]:
anfrage = pd.DataFrame()
anfrage['father'] = [0, 69.8]   
#  rechte Seite ist auch dann ein array, wenn nur ein Wert berechnet werden soll
anfrage
Out[8]:
father
0 0.0
1 69.8
In [9]:
res.get_prediction(anfrage).summary_frame()
Out[9]:
mean mean_se mean_ci_lower mean_ci_upper obs_ci_lower obs_ci_upper
0 38.362581 3.308374 31.861862 44.863300 30.313002 46.412160
1 69.529859 0.114624 69.304631 69.755087 64.777270 74.282447
  • mean: Wert, der für den Sohn im Mittel zu erwarten ist
  • mean bei 0 ist das Intercept
  • mean_se: Standardabweichung für mean
  • die vier anderen Werte sind untere bzw. obere Grenzen für Konfidenzintervalle
  • mean_ci ist das Konfidenzintervall für den mittleren zu erwartenden Wert
  • obs_ci ist das Konfidentintervall für den individuelle beobachteten Wert (engl: "observed")
  • mean_ci_lower und mean_ci_upper begrenzen die blaue Kurve in dem regplot
  • obs_ci_lower und obs_ci_upper begrenzen einen Bereich, der 95% der Beobachtungen enthält
  • wir malen den mal hin
In [10]:
anfrage = pd.DataFrame()
anfrage['father'] = np.arange(625, 775) / 10
vorhersage = res.get_prediction(anfrage).summary_frame()
vorhersage
Out[10]:
mean mean_se mean_ci_lower mean_ci_upper obs_ci_lower obs_ci_upper
0 66.270244 0.336018 65.609992 66.930496 61.477301 71.063187
1 66.314896 0.331504 65.663515 66.966277 61.523167 71.106625
2 66.359548 0.326997 65.717023 67.002074 61.569015 71.150082
3 66.404201 0.322498 65.770515 67.037886 61.614845 71.193556
4 66.448853 0.318007 65.823992 67.073714 61.660657 71.237049
... ... ... ... ... ... ...
145 72.744822 0.391826 71.974912 73.514731 67.935546 77.554097
146 72.789474 0.396418 72.010542 73.568406 67.978746 77.600202
147 72.834126 0.401014 72.046162 73.622090 68.021927 77.646325
148 72.878778 0.405615 72.081775 73.675782 68.065091 77.692466
149 72.923431 0.410219 72.117379 73.729483 68.108237 77.738624

150 rows × 6 columns

In [11]:
ax = sns.scatterplot(df, x='father', y='childHeight')
sns.lineplot(x=anfrage.father, y=vorhersage['mean'], ax=ax)
sns.lineplot(x=anfrage.father, y=vorhersage.obs_ci_lower, ax = ax, color='orange')
sns.lineplot(x=anfrage.father, y=vorhersage.obs_ci_upper, ax=ax, color='orange');
No description has been provided for this image

Die orangen Linien sind die untere bzw. obere Vertrauensgrenze für die Einzelbeobachtungen

Anderes Konfidenzniveau $1-\alpha$

In [12]:
res.get_prediction(anfrage).summary_frame(alpha=0.02)
Out[12]:
mean mean_se mean_ci_lower mean_ci_upper obs_ci_lower obs_ci_upper
0 66.270244 0.336018 65.485924 67.054563 60.576660 71.963827
1 66.314896 0.331504 65.541114 67.088678 60.622755 72.007037
2 66.359548 0.326997 65.596286 67.122810 60.668827 72.050269
3 66.404201 0.322498 65.651440 67.156961 60.714879 72.093523
4 66.448853 0.318007 65.706574 67.191132 60.760908 72.136797
... ... ... ... ... ... ...
145 72.744822 0.391826 71.830239 73.659404 67.031837 78.457807
146 72.789474 0.396418 71.864173 73.714775 67.074763 78.504184
147 72.834126 0.401014 71.898096 73.770156 67.117669 78.550584
148 72.878778 0.405615 71.932010 73.825547 67.160553 78.597004
149 72.923431 0.410219 71.965914 73.880948 67.203415 78.643446

150 rows × 6 columns

Beispiel: Fische¶

  • Fische werden gezüchtet. In den ersten 24 Monaten wurden die folgenden Daten erhoben
  • Diesen Daten werden benutzt, um das Wachstum der nächsten Generation zu prognostizieren
In [13]:
rng = np.random.default_rng(123)
N = 70
In [14]:
df = pd.DataFrame()
df['Monat'] = rng.choice(np.arange(4, 25), size=N)
df['Höhe'] = 4.5*df.Monat + stats.norm(0, 2.2).rvs(size=N, random_state=rng)
df['Gewicht'] = 65*df.Monat + 4*df.Höhe + stats.norm(0, 50).rvs(size=N, random_state=rng)
#df.to_csv('fische.csv', index=False)
In [15]:
df = pd.read_csv('fische.csv')
df.describe()
Out[15]:
Monat Höhe Gewicht
count 70.000000 70.000000 70.000000
mean 13.585714 61.335852 1131.866720
std 5.619564 25.491528 481.798450
min 4.000000 15.600673 238.754392
25% 10.000000 42.654557 798.817581
50% 13.000000 58.437946 1130.206032
75% 18.000000 79.635416 1457.015505
max 24.000000 109.853556 2089.839600
  • Gewicht in g
  • Höhe in mm
  • Ein Züchter hat 1200 Fische in seinen Teichen, die alle gleichzeitig geschlüpft sind
  • Frage: Konfidenzintervall für das Gesamtgewicht dieser Fische nach 18 Monaten zum Konfidenzniveau 95%?
  • Frage: Wie muss das Netz gewählt werden, um nach 18 Monaten 97.5% der Fische zu fangen?
In [16]:
formel1 = 'Gewicht ~ Monat'
modell1 = smf.ols(formel1, df)
In [17]:
res = modell1.fit()
res.summary()
Out[17]:
OLS Regression Results
Dep. Variable: Gewicht R-squared: 0.989
Model: OLS Adj. R-squared: 0.989
Method: Least Squares F-statistic: 6398.
Date: Mon, 24 Jun 2024 Prob (F-statistic): 5.38e-69
Time: 21:06:30 Log-Likelihood: -371.83
No. Observations: 70 AIC: 747.7
Df Residuals: 68 BIC: 752.2
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -26.7757 15.660 -1.710 0.092 -58.024 4.472
Monat 85.2839 1.066 79.986 0.000 83.156 87.412
Omnibus: 2.313 Durbin-Watson: 2.229
Prob(Omnibus): 0.315 Jarque-Bera (JB): 1.580
Skew: -0.319 Prob(JB): 0.454
Kurtosis: 3.366 Cond. No. 38.8


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
In [18]:
anfrage = pd.DataFrame()
anfrage['Monat'] = [18]
res.get_prediction(anfrage).summary_frame()
Out[18]:
mean mean_se mean_ci_lower mean_ci_upper obs_ci_lower obs_ci_upper
0 1508.33412 7.585591 1493.19731 1523.470931 1407.869929 1608.798312
  • untere Vertrauensgrenze für das Gesamtgewicht von 1200 Fischen in kg:
In [19]:
1200 * 1493 / 1000
Out[19]:
1791.6
  • obere Vertrauensgrenze für das Gesamtgewicht von 1200 Fischen in kg:
In [20]:
1200 *  1523 / 1000
Out[20]:
1827.6

Mit 97.5% Sicherheit werden mindestens 1791 kg Fisch geerntet

In [21]:
formel2 = 'Höhe ~ Monat'
modell2 = smf.ols(formel2, df)
res = modell2.fit()
res.summary()
Out[21]:
OLS Regression Results
Dep. Variable: Höhe R-squared: 0.993
Model: OLS Adj. R-squared: 0.993
Method: Least Squares F-statistic: 9404.
Date: Mon, 24 Jun 2024 Prob (F-statistic): 1.24e-74
Time: 21:06:30 Log-Likelihood: -152.73
No. Observations: 70 AIC: 309.5
Df Residuals: 68 BIC: 313.9
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -0.0702 0.685 -0.103 0.919 -1.436 1.296
Monat 4.5199 0.047 96.974 0.000 4.427 4.613
Omnibus: 0.052 Durbin-Watson: 1.213
Prob(Omnibus): 0.974 Jarque-Bera (JB): 0.128
Skew: -0.061 Prob(JB): 0.938
Kurtosis: 2.830 Cond. No. 38.8


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
In [22]:
res.get_prediction(anfrage).summary_frame()
Out[22]:
mean mean_se mean_ci_lower mean_ci_upper obs_ci_lower obs_ci_upper
0 81.287975 0.331597 80.626284 81.949667 76.896279 85.679671

Um 97.5% der Fische zu fangen, muss das Netz so beschaffen sein, dass ein Fisch der Höhe 76.9mm nicht hindurch schlüpft

Mehrere erklärende Variablen¶

Lineares Modell mit einer erklärten und mehreren erklärenden Variablen¶

$$ y = m_1 \cdot x_1 + m_2 \cdot x_2 + \dots + m_n \cdot x_n + b $$

$y$ ist die erklärte und die $x_i$ sind die erklärenden Variablen

Beispiel: Körpergröße der Söhne hängt von der Körpergröße von Vater und Mutter ab

In [23]:
df = pd.read_csv('galton.csv')
In [24]:
formel = 'childHeight ~ father + mother'

Diese Formel hat 3 Unbekannte:

  • den Koeffizienten von father
  • den Koeffizienten von mother
  • den Ordinatenabschnitt
In [25]:
modell = smf.ols(formel, df)
In [26]:
res = modell.fit()
In [27]:
res.summary()
Out[27]:
OLS Regression Results
Dep. Variable: childHeight R-squared: 0.238
Model: OLS Adj. R-squared: 0.235
Method: Least Squares F-statistic: 74.62
Date: Mon, 24 Jun 2024 Prob (F-statistic): 6.25e-29
Time: 21:06:30 Log-Likelihood: -1080.7
No. Observations: 481 AIC: 2167.
Df Residuals: 478 BIC: 2180.
Df Model: 2
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 19.3128 4.095 4.716 0.000 11.266 27.359
father 0.4176 0.046 9.154 0.000 0.328 0.507
mother 0.3288 0.045 7.258 0.000 0.240 0.418
Omnibus: 10.653 Durbin-Watson: 1.592
Prob(Omnibus): 0.005 Jarque-Bera (JB): 14.542
Skew: -0.200 Prob(JB): 0.000695
Kurtosis: 3.752 Cond. No. 3.69e+03


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 3.69e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
  • Regressiongleichung: childHeight = 0.4176 * father + 0.3288 * mother + 19.3128
  • alle drei Koeffizienten haben statistisch signifikanten Einfluss

Zum Vergleich

In [28]:
formel2 = 'childHeight ~ father'
modell2 = smf.ols(formel2, df)
res = modell2.fit()
In [29]:
res.summary()
Out[29]:
OLS Regression Results
Dep. Variable: childHeight R-squared: 0.154
Model: OLS Adj. R-squared: 0.152
Method: Least Squares F-statistic: 87.17
Date: Mon, 24 Jun 2024 Prob (F-statistic): 3.74e-19
Time: 21:06:30 Log-Likelihood: -1105.8
No. Observations: 481 AIC: 2216.
Df Residuals: 479 BIC: 2224.
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 38.3626 3.308 11.596 0.000 31.862 44.863
father 0.4465 0.048 9.337 0.000 0.353 0.540
Omnibus: 8.610 Durbin-Watson: 1.468
Prob(Omnibus): 0.014 Jarque-Bera (JB): 12.731
Skew: -0.110 Prob(JB): 0.00172
Kurtosis: 3.766 Cond. No. 2.08e+03


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.08e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
  • Regressionsgleichung childHeight = 0.4465 * father + 38.3626
  • beide Koeffizienten haben statistisch signifikanten Einfluss
  • Das erste Modell ist genauer, denn dort ist der Wert von R-squared höher

Füge einen irrelevanten Term in den Datensatz ein

In [30]:
rng = np.random.default_rng()
df['Kontonummer'] = rng.integers(1000, 99999, size=481)
df.head()
Out[30]:
family father mother midparentHeight children childNum gender childHeight Kontonummer
0 001 78.5 67.0 75.43 4 1 male 73.2 85251
1 002 75.5 66.5 73.66 4 1 male 73.5 37646
2 002 75.5 66.5 73.66 4 2 male 72.5 66951
3 003 75.0 64.0 72.06 2 1 male 71.0 75302
4 004 75.0 64.0 72.06 5 1 male 70.5 65773
In [31]:
formel3 = 'childHeight ~ father + mother + Kontonummer'
modell3 = smf.ols(formel3, df)
res = modell3.fit()
In [32]:
res.summary()
Out[32]:
OLS Regression Results
Dep. Variable: childHeight R-squared: 0.239
Model: OLS Adj. R-squared: 0.234
Method: Least Squares F-statistic: 49.81
Date: Mon, 24 Jun 2024 Prob (F-statistic): 5.07e-28
Time: 21:06:30 Log-Likelihood: -1080.5
No. Observations: 481 AIC: 2169.
Df Residuals: 477 BIC: 2186.
Df Model: 3
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 19.0833 4.114 4.638 0.000 10.999 27.168
father 0.4188 0.046 9.167 0.000 0.329 0.509
mother 0.3292 0.045 7.261 0.000 0.240 0.418
Kontonummer 2.334e-06 3.77e-06 0.619 0.536 -5.08e-06 9.75e-06
Omnibus: 10.731 Durbin-Watson: 1.598
Prob(Omnibus): 0.005 Jarque-Bera (JB): 14.555
Skew: -0.204 Prob(JB): 0.000691
Kurtosis: 3.748 Cond. No. 2.25e+06


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.25e+06. This might indicate that there are
strong multicollinearity or other numerical problems.
  • Die Koeffizienten sind etwas verändert gegenüber dem ersten Modell
  • Der Unterschied ist nicht signifikant
  • Wert von R-squared unverändert gegenüber dem ersten Modell
  • Von den hier vorgestellten Modellen ist das erste das beste