[PYTHON] Einführung in die statistische Modellierung für die Datenanalyse Generalized Linear Model (GLM)
Fortsetzung der vorherigen Zeit (Kapitel 2), Kapitel 3. Untersuchung eines statistischen Modells, das ** erklärende Variablen ** enthält.
3.1 Beispiel: Wenn die durchschnittliche Anzahl der Samen von Individuum zu Individuum variiert
Als Unterschied zur vorherigen Zeit werden die Körpergröße $ x_i $, die eines der Attribute des Individuums ist, und das Vorhandensein oder Fehlen einer Kontrolle (nichts tun: ** C **, Dünger hinzufügen: ** T **) als erklärende Variablen verwendet.
Daten lesen
[Hier] Laden Sie data3a.csv in Kapitel 3 von () herunter und verwenden Sie es.
Die diesmal verwendeten Daten werden als DataFrame-Typ bezeichnet.
# d <- read.csv("data3a.csv")
>>> import pandas
>>> d = pandas.read_csv("data3a.csv")
>>> d
y x f
0 6 8.31 C
1 6 9.44 C
2 6 9.50 C
3 12 9.07 C
...(Weggelassen)...
98 7 10.86 T
99 9 9.97 T
[100 rows x 3 columns]
# d$x
>>> d.x
0 8.31
1 9.44
2 9.50
...(Weggelassen)...
98 10.86
99 9.97
Name: x, Length: 100, dtype: float64
# d$y
>>> d.y
0 6
1 6
2 6
...(Weggelassen)...
98 7
99 9
Name: y, Length: 100, dtype: int64
# d$f
>>> d.f = d.f.astype('category')
>>> d.f
0 C
1 C
2 C
3 C
...(Weggelassen)...
98 T
99 T
Name: f, Length: 100, dtype: category
Categories (2, object): [C < T]
In R heißt Spalte f ** Faktor ** (Faktor), In Python (Pandas) ist es eine Kategorie.
So überprüfen Sie den Datentyp
# class(d)Kann für einen Moment nicht bestätigt werden.
# class(d$y)Yara Klasse(d$x)Yara Klasse(d$f)
>>> d.y.dtype
dtype('int64')
>>> d.x.dtype
dtype('float64')
>>> d.f.dtype
dtype('O')
Eine Übersicht über den Datenrahmen finden Sie unter Serie (wie beim letzten Mal).
# summary(d)
>>> d.describe()
y x
count 100.000000 100.000000
mean 7.830000 10.089100
std 2.624881 1.008049
min 2.000000 7.190000
25% 6.000000 9.427500
50% 8.000000 10.155000
75% 10.000000 10.685000
max 15.000000 12.400000
>>> d.f.describe()
count 100
unique 2
top T
freq 50
Name: f, dtype: object
Visualisierung
Ich konnte keinen Weg finden, das Streudiagramm in einer Zeile wie R farblich zu kennzeichnen. .. ..
# plot(d$x, d$y, pch=c(21, 19)[d$f])
>>> plt.plot(d.x[d.f == 'C'], d.y[d.f == 'C'], 'bo')
[<matplotlib.lines.Line2D at 0x1184d0490>]
>>> plt.plot(d.x[d.f == 'T'], d.y[d.f == 'T'], 'ro')
[<matplotlib.lines.Line2D at 0x111b58dd0>]
>>> plt.show()
>>> d.boxplot(column='y', by='f')
<matplotlib.axes._subplots.AxesSubplot at 0x11891a9d0>
Was Sie aus diesen Zahlen sehen können
- Es scheint, dass die Anzahl der Samen $ y $ mit zunehmender Körpergröße $ x_i $ zunimmt.
- Die Wirkung von Dünger scheint nicht $ f $ zu sein
Statistisches Modell der Poisson-Regression
Erstellen Sie ein statistisches Modell, das in der Lage zu sein scheint, die Seed-Number-Daten auszudrücken, bei denen es sich um die Count-Daten handelt.
Statistisches Modell, das von der Körpergröße $ x_i $ des einzelnen $ i $ abhängt
** Erklärende Variable **: $ x_i $ ** Antwortvariable **: $ y_i $
Die Wahrscheinlichkeit $ p (y_i | \ lambda_i) $, dass die Anzahl der Samen $ y_i $ in einem einzelnen $ i $ ist, folgt der Poisson-Verteilung.
p(y_i|\lambda_i) = \frac{\lambda _i ^{y_i} \exp (-\lambda_i)}{y_i !}
Annehmen. Hier ist die durchschnittliche Anzahl von Samen $ \ lambda_i $ eines einzelnen $ i $
\begin{eqnarray}
\lambda_i &=& \exp(\beta_1 + \beta_2 x_i )
\log \lambda_i &=& \beta_1 + \beta_2 x_i
\end{eqnarray}
Nehme an, dass. $ \ Beta_1 $ heißt ** section **, $ \ beta_2 $ heißt ** tilt **, Die rechte Seite $ \ beta_1 + \ beta_2 x_i $ heißt ** linearer Prädiktor **.
Wenn ($ \ lambda_i $ function) = (linearer Prädiktor), Die Funktion auf der linken Seite heißt ** Link-Funktion **.
Poisson Rückkehr
Die Poisson-Regression findet $ \ beta_1, \ beta_2 $, die die folgende Gleichung maximiert.
\log L(\beta_1, \beta_2) = \sum_i \log \frac{\lambda_i^{y_i}\exp(-\lambda_i)}{y_i !}
In R scheint es ein Schuss zu sein, wenn Sie so etwas wie "glm" verwenden. Mit Python muss ich endlos mit der Bibliothek "statsmodels" arbeiten ... (Referenz)
# fit <- glm(y ~ x, data=d, familiy=poisson)
>>> import statsmodels.api as sm
>>> import statsmodel.formula.api as smf
>>> model = smf.glm('y ~ x', data=d family=sm.families.Poisson())
>>> results = model.fit()
>>> print(results.summary())
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 98
Model Family: Poisson Df Model: 1
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -235.39
Date:Mond, 01 6 2015 Deviance: 84.993
Time: 23:06:45 Pearson chi2: 83.8
No. Iterations: 7
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
const 1.2917 0.364 3.552 0.000 0.579 2.005
x 0.0757 0.036 2.125 0.034 0.006 0.145
==============================================================================
>>> results.llf
-235.38625076986077
Es scheint, dass es Intercept in R heißt, In Python scheint die untere Konstante $ \ beta_1 $ (Abschnitt) und x $ \ beta_2 $ (Neigung) darzustellen.
Mit anderen Worten, die wahrscheinlichsten Schätzungen sind $ \ beta_1 = 1,2917 $ und $ \ beta_2 = 0,0757 $.
Der Standardfehler gibt den ** Standardfehler ** an, der die Variation von $ \ beta_1 und \ beta_2 $ angibt.
Wenn dieselbe Anzahl unterschiedlicher Daten mit derselben Erhebungsmethode erneut abgetastet wird, ändert sich auch der wahrscheinlichste Schätzbetrag und der Variationsgrad.
z ist eine Statistik namens ** z-Wert **, die die wahrscheinlichste Schätzung geteilt durch den Standardfehler ist.
Es scheint, dass dieses z verwendet werden kann, um das ** Wald-Konfidenzintervall ** zu schätzen. Anscheinend enthält dieses Konfidenzintervall nicht 0 in dem Wert, den die wahrscheinlichste Schätzung annehmen kann siehe. html # 08)
Vorhersage durch Poisson-Regressionsmodell
Das Schätzergebnis der Poisson-Regression wird verwendet, um die durchschnittliche Anzahl von Samen $ \ lambda $ bei Körpergröße $ x $ vorherzusagen.
\lambda = \exp(1.29 + 0.0757x)
Wird mit dargestellt.
# plot(d$x, d$y, pch=c(21, 19)[d$f])
# xx <- seq(min(d$x), max(d$x), length=100)
# lines(xx, exp(1.29 + 0.0757 * xx), lwd = 2)
>>> plt.plot(d.x[d.f == 'C'], d.y[d.f == 'C'], 'bo')
>>> plt.plot(d.x[d.f == 'T'], d.y[d.f == 'T'], 'ro')
>>> xx = [d.x.min() + i * ((d.x.max() - d.x.min()) / 100.) for i in range(100)]
>>> yy = [numpy.exp(1.29 + 0.0757 * i )for i in xx]
>>> plt.plot(xx, yy)
>>> plt.show()
3.5 Statistisches Modell mit erklärenden Variablen des Faktortyps
Betrachten Sie ein statistisches Modell, das die Düngemittelanwendung $ f_i $ als erklärende Variable enthält.
In GLM wird ** Dummy-Variable ** als erklärende Variable für den Faktortyp verwendet. Das heißt, der Durchschnittswert des Modells,
\begin{eqnarray}
\lambda_i &=& \exp (\beta_1 + \beta_3 d_i) \\
\therefore d_i &=& \left\{ \begin{array}{ll}
0 & (f_i =Für C.) \\
1 & (f_i =Für T.)
\end{array} \right.
\end{eqnarray}
Wird geschrieben werden.
# fit.f <- glm(y ~ f, data=d, family=poisson)
>>> model = smf.glm('y ~ f', data=d, familiy=sm.families.Poisson())
>>> model.fit().summary()
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 98
Model Family: Poisson Df Model: 1
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -237.63
Date:Holz, 04 6 2015 Deviance: 89.475
Time: 17:20:11 Pearson chi2: 87.1
No. Iterations: 7
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept 2.0516 0.051 40.463 0.000 1.952 2.151
f[T.T] 0.0128 0.071 0.179 0.858 -0.127 0.153
==============================================================================
>>> model.llf
-237.62725696068682
Wenn aus diesem Ergebnis $ f_i $ des Individuums $ i $ C ist,
\lambda_i = \exp(2.05 + 0.0128 * 0) = \exp(2.05) = 7.77
Wenn T.
\lambda_i = \exp(2.05 + 0.0128 * 1) = \exp(2.0628) = 7.87
Es wird vorausgesagt, dass die durchschnittliche Anzahl der Samen nur geringfügig zunimmt, wenn Dünger verabreicht wird.
Es kann gesagt werden, dass die logarithmische Wahrscheinlichkeit (llf) nicht anwendbar ist, weil sie kleiner geworden ist.
3.6 Statistisches Modell mit erklärenden Variablen für den quantitativen Typ + Faktortyp
Stellen Sie sich ein statistisches Modell vor, das sowohl die individuelle Körpergröße $ x_i $ als auch die Befruchtungsbehandlung $ f_i $ umfasst.
\log \lambda_i = \beta_1 + \beta_2 x_i + \beta_3 d_i
>>> model = smf.glm('y ~ x + f', data=d, family=sm.families.Poisson())
>>> result = model.fit()
>>> result.summary()
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 97
Model Family: Poisson Df Model: 2
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -235.29
Date:Holz, 04 6 2015 Deviance: 84.808
Time: 17:31:34 Pearson chi2: 83.8
No. Iterations: 7
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept 1.2631 0.370 3.417 0.001 0.539 1.988
f[T.T] -0.0320 0.074 -0.430 0.667 -0.178 0.114
x 0.0801 0.037 2.162 0.031 0.007 0.153
==============================================================================
>>> result.llf
-235.29371924249369
Da die maximale Log-Wahrscheinlichkeit am größten ist, kann gesagt werden, dass sie besser passt als die beiden oben genannten statistischen Modelle.
Einfaches Verständnis der logarithmischen Verknüpfungsfunktion: zu multiplizierender Effekt
Durch Transformation des obigen statistischen Modells vom quantitativen Modell + Faktortyp
\begin{eqnarray}
\log \lambda_i &=& \beta_1 + \beta_2 x_i + \beta_3 d_i \\
\lambda_i &=& \exp(\beta_1 + \beta_2 x_i + \beta_3 d_i) \\
\lambda_i &=& \exp(\beta_1) * \exp(\beta_2 x_i) * \exp(\beta_3 d_i) \\
\lambda_i &=& \exp(Konstante) * \exp(Einfluss der Körpergröße) * \exp(Wirkung der Befruchtungsbehandlung)
\end{eqnarray}
Es ist ersichtlich, dass sich die Effekte vervielfachen.
Wenn keine Verknüpfungsfunktion vorhanden ist, wird diese als ** konstante Verknüpfungsfunktion ** bezeichnet.
In GLM
Der Fall der Normalverteilung und der oralen Verknüpfungsfunktion wird als ** lineares Modell ** (** lineares Modell, LM ) oder ** allgemeines lineares Modell ** ( allgemeines lineares Modell **) bezeichnet.
Recommended Posts