[Statistischer Test 2. Klasse / quasi 1. Klasse] Regressives Analysetraining mit Python (2)

Einführung

[Statistischer Test 2. Klasse / quasi 1. Klasse] Das Training zur Regressionsanalyse in Python (1) - Qiita ist eine Fortsetzung. Dieses Mal werden wir tatsächlich eine multiple Regressionsanalyse mit zwei oder mehr unabhängigen Variablen (erklärenden Variablen) unter Verwendung von Python versuchen und eine statistische Interpretation in Betracht ziehen. Auch dieses Mal werde ich Bibliotheken wie Numpy und Pandas verwenden, um den Berechnungsteil anhand der mathematischen Formeln selbst zu schreiben.

Nachschlagewerk

Katsuya Hasegawa "Multivariate Analyse, die von Grund auf gut verstanden werden kann" Technical Review Company (2004) Kapitel 2

Als Referenz für den theoretischen Teil verweise ich auch auf Folgendes. Abteilung für Statistik der Fakultät für Geisteswissenschaften der Universität Tokio "Einführung in die Statistik (Grundlegende Statistik I)" Tokyo University Press (1991) Kapitel 13

Gegenstand

Ich wünschte, ich könnte es mit den JR-Routendaten tun, die ich zuletzt verwendet habe, aber es schien schwierig, die unabhängigen Variablen zu erfassen, und so entschied ich mich, dieses Mal andere Daten zu verwenden. Homepage des Statistikbüros / Japanische Statistik 2019 - Kapitel 2 Bevölkerung / Haushalt

Die folgenden Daten für jede Präfektur werden kombiniert und zur Verwendung im CSV-Format verarbeitet.

―― 2-2 Bevölkerung nach Präfektur und Bevölkerungszunahme / -abnahme ―― 2-13 Tagesbevölkerung nach Präfektur und Anzahl der Personen, die außerhalb des Hauses arbeiten / zur Schule gehen ―― 2-14 Anzahl der Personen, die nach Präfektur ein- und ausziehen ―― 2-16 Anzahl der Geburten / Todesfälle und Anzahl der Ehen / Scheidungen nach Präfektur

** (2020/01/07) Nur die Spalten Bevölkerung und Bevölkerungsverhältnis waren um eine Zeile versetzt, daher habe ich sie korrigiert. ** ** **

Präfekturen,Population,Population性比,Population増減率,昼夜間Population比率,Anzahl der Übermigranten,Geburtenrate,Mortalität,Natürliche Zu- / Abnahmerate,Heiratsrate,Scheidungsrate
Hokkaido,5320,89.1,-0.455,99.945,-6569,6.4,11.8,-5.4,4.5,1.92
Aomori,1278,88.6,-0.965,99.849,-6075,6.3,13.8,-7.5,4.0,1.64
:
:
Okinawa,1626,88.5,0.582,99.968,-1112,11.3,8.4,3.0,5.7,2.44

――Bevölkerung (Einheit: 1.000 Personen) ist die geschätzte Bevölkerung im Jahr 2017.

Nach wie vor speichern wir dies mit dem ** Zeichencode als UTF-8 und dem Dateinamen ** populations.csv.

Trainieren

Im Folgenden haben wir den Vorgang mit Python 3.6.8 bestätigt.

Dateneingabe

import numpy as np
import pandas as pd
df = pd.read_csv("population.csv")
print(df) #Überprüfen Sie den Inhalt der importierten Daten

Multivariate Regressionsanalyse (Abschnitte 2-3)

Lassen Sie uns anhand dieser Daten über die Regressionsanalyse nachdenken, wenn zwei oder mehr unabhängige Variablen vorhanden sind, dh "multiple Regressionsanalyse". Es gibt verschiedene mögliche Kombinationen von unabhängigen Variablen und abhängigen Variablen, aber zunächst scheint sich eine leicht verständliche Beziehung hinsichtlich der Bestätigung der Methode zu ergeben.

Lass es uns versuchen mit.

[^ Geburtenrate]: Anzahl der Geburten und Todesfälle pro 1.000 Einwohner

Stellen Sie sich ein Regressionsmodell von $ y = \ beta_0 + \ beta_1 x_1 + \ beta_2 x_2 $ vor, bei dem die Geburtenrate $ x_1 $, die Sterblichkeitsrate $ x_2 $ und die Bevölkerungswachstumsrate $ y $ beträgt. Wie für eine Variable, partieller Regressionskoeffizient $ \ beta_0, um die Summe der Quadrate der Residuen jedes Datenpunkts (Differenz zwischen dem beobachteten Wert (realisierter Wert) und dem theoretischen Wert (vorhergesagter Wert) der abhängigen Variablen) zu minimieren. , \ beta_1, \ beta_2 $. Jeder Datenpunkt (beobachteter Wert) $ (x_1 ^ {(1)}, x_2 ^ {(1)}, y ^ {(1)}), ..., (x_1 ^ {(n)}, x_2 ^ { Für (n)}, y ^ {(n)}) $ ist der theoretische Wert der abhängigen Variablen $ \ hat {y} ^ {(i)} = \ beta_0 + \ beta_1 x_1 ^ {(i)} + \ beta_2 Wenn x_2 ^ {(i)} $, ist die Summe der Quadrate der Residuen $ L $

\begin{align}
L &= \sum_{i=1}^n (y^{(i)} - \hat{y}^{(i)})^2 \\
&= \sum_{i=1}^n (y^{(i)} - \beta_0 - \beta_1 x_1^{(i)} - \beta_2 x_2^{(i)})^2
\end{align}

Es wird sein. Nehmen wir für $ \ beta_0, \ beta_1, \ beta_2 $ an, dass die Lösung, wenn der partielle Differentialkoeffizient 0 ist, der geschätzte Wert $ \ hat \ beta_0, \ hat \ beta_1, \ hat \ beta_2 $ ist.

\begin{align}
\hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y}
\end{align}

(Es wird angenommen, dass $ ^ \ Top $ die Transposition einer Matrix / eines Vektors darstellt). jedoch

\begin{align}
\hat{\boldsymbol{\beta}} &= (\hat\beta_0, \hat\beta_1, \hat\beta_2)^\top \\
\boldsymbol{X} &= \begin{pmatrix}
1 & x_1^{(1)} & x_2^{(1)} \\
\vdots & \vdots & \vdots \\
1 & x_1^{(n)} & x_2^{(n)} \\
\end{pmatrix}
\\
\boldsymbol{y} &= (y^{(1)}, ..., y^{(n)})^\top
\end{align}

Nehmen wir also an, dass $ \ boldsymbol {X} ^ \ top \ boldsymbol {X} $ eine inverse Matrix hat. Versuchen Sie anhand dieses Ergebnisses, den partiellen Regressionskoeffizienten zu berechnen

n = len(df)
X = np.c_[np.ones(n), df[["Geburtenrate", "Mortalität"]].values]
y = df["Bevölkerungswachstums- / -abnahmerate"]
b = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print(b) # Output: [ 1.00978055  0.10175178 -0.18246347]

Daher wird $ \ hat {y} = 1,010 + 0,102 x_1 - 0,182 x_2 $ als Regressionsgleichung des theoretischen Wertes der Bevölkerungszunahme- / -abnahmerate von $ \ hat {y} $ erhalten.

――Wenn die Geburtenrate pro 1.000 Personen um 1 steigt, erhöht sich die jährliche Bevölkerungs- / Abnahmerate um 0,102 Punkte. ――Wenn die Sterblichkeitsrate pro 1.000 Menschen um 1 steigt, sinkt die jährliche Zunahme / Abnahme der Bevölkerung um 0,182 Punkte.

Wenn Sie versuchen, den Entscheidungskoeffizienten wie bei einer Variablen zu ermitteln

TSS = ((y - y.mean()) ** 2).sum()        #Gesamtschwankung
ESS = ((X.dot(b) - y.mean()) ** 2).sum() #Regressive Fluktuation
RSS = ((y - X.dot(b)) ** 2).sum()        #Restschwankung
R_2 = ESS / TSS                          #Entscheidungskoeffizient
print(R_2) # 0.9014307808651625

Es wurde gebeten, ungefähr 0,901 zu sein. Ist dieser Wert groß oder klein? Ich werde es später versuchen, wenn ich eine andere unabhängige Variable nehme, also vergleiche die Werte.

Ableitung der Formel zur Berechnung des partiellen Regressionskoeffizienten

Die Details werden im obigen Nachschlagewerk übersprungen, aber lassen Sie uns die Formel für die partiellen Regressionskoeffizienten $ \ beta_0, \ beta_1, \ beta_2 $ ableiten. Es sollte möglich sein, es zu erhalten, indem man es teilweise normal differenziert und als simultane Gleichung löst. Da jedoch die Formeltransformation, die selbst mit einer unabhängigen Variablen problematisch war, schwieriger wird, sollten wir sie mit einer klaren Sicht unter Verwendung der "Differenzierung des Vektors" organisieren.

"Differenzierung nach Vektor / Differenzierung nach Matrix" Offizielle Zusammenfassung --Qiita

Verallgemeinern Sie hier die Anzahl der unabhängigen Variablen auf $ p $ und finden Sie die Formel für die partiellen Regressionskoeffizienten $ \ beta_0, \ beta_1, ..., \ beta_p $.

\begin{align}
\boldsymbol{\beta} &= (\beta_0, \beta_1, ...,  \beta_p)^\top \\
\boldsymbol{x}^{(i)} &= (1, x_1^{(i)}, ..., x_p^{(i)})^\top \\
\boldsymbol{X} &= (\boldsymbol{x}^{(1)}, ..., \boldsymbol{x}^{(n)})^\top = \begin{pmatrix}
1 & x_1^{(1)} & \cdots & x_p^{(1)} \\
\vdots & \vdots & & \vdots \\
1 & x_1^{(n)} & \cdots & x_p^{(n)} \\
\end{pmatrix}
\\
\boldsymbol{y} &= (y^{(1)}, ..., y^{(n)})^\top \\
\epsilon^{(i)}
&= y^{(i)} - \hat{y}^{(i)} = y^{(i)} - (\beta_0 + \beta_1 x_1^{(i)} + ... + \beta_p x_p^{(i)} ) \\
&= y^{(i)} - \boldsymbol{\beta}^\top \boldsymbol{x}^{(i)} \\
\boldsymbol{\epsilon} &= (\epsilon^{(1)}, ..., \epsilon^{(n)})^\top \\ 
&= \boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}\\
\end{align}

Wird besorgt. Beachten Sie, dass wir der Einfachheit halber $ \ boldsymbol {x} ^ {(i)} $ eine Dimension hinzugefügt haben, die dem konstanten Term entspricht. Wir werden auch den Nullvektor (einen Vektor mit allen Nullkomponenten) als $ \ boldsymbol {0} $ schreiben. Das Setzen des partiellen Differentialkoeffizienten auf 0 für $ \ beta_0, \ beta_1, ..., \ beta_p $ bedeutet $ \ frac {\ partielles L} {\ partielles \ Boldsymbol {\ beta}} = \ Boldsymbol {0} Es ist dasselbe wie $ setzen.

Suchen Sie $ \ boldsymbol {\ beta} $ anhand der entsprechenden Formel auf der obigen Seite. Es ist eine leicht sorgfältig geschriebene Version der Formel mit den folgenden Inhalten. Einführung in die Statistik - Kapitel 7 | Garaku Taton Chenkan> Trivia Room> Trivia Corner

\begin{align}
L &= \sum_{i=1}^n (\epsilon^{(i)})^2  = \boldsymbol{\epsilon}^\top \boldsymbol{\epsilon} \\
&= (\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta})^\top (\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}) \\
&= \boldsymbol{y}^\top \boldsymbol{y} - \boldsymbol{y}^\top (\boldsymbol{X} \boldsymbol{\beta}) - (\boldsymbol{X} \boldsymbol{\beta})^\top \boldsymbol{y} + \boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} \\
&= \boldsymbol{y}^\top \boldsymbol{y} - 2 \boldsymbol{y}^\top \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} \\
\frac{\partial L}{\partial \boldsymbol{\beta}} &= \frac{\partial}{\partial \boldsymbol{\beta}} \boldsymbol{y}^\top \boldsymbol{y} - 2 \frac{\partial}{\partial \boldsymbol{\beta}}(\boldsymbol{y}^\top \boldsymbol{X}) \boldsymbol{\beta} + \frac{\partial}{\partial \boldsymbol{\beta}}\boldsymbol{\beta}^\top (\boldsymbol{X}^\top \boldsymbol{X}) \boldsymbol{\beta} \\
&= -2 (\boldsymbol{y}^\top \boldsymbol{X})^\top + (\boldsymbol{X}^\top \boldsymbol{X} + (\boldsymbol{X}^\top \boldsymbol{X})^\top) \boldsymbol{\beta} \\
&= -2 \boldsymbol{X}^\top \boldsymbol{y} + 2\boldsymbol{X}^\top \boldsymbol{X}\boldsymbol{\beta}
\end{align}

Wenn es also eine inverse Matrix in $ \ boldsymbol {X} ^ \ top \ boldsymbol {X} $ gibt, ist $ \ frac {\ partielles L} {\ partielles \ boldsymbol {\ beta}} = \ boldsymbol {0} $ Die Lösung von sei der geschätzte Wert $ \ hat {\ boldsymbol {\ beta}} $

\begin{align}
\hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y}
\end{align}

Es kann erhalten werden.

Statistische Theorie der Regressionsanalyse (Abschnitte 2-4)

Nachdem wir nun eine Regressionsgleichung haben, betrachten wir statistisch, "wie zuverlässig die resultierende Regressionsgleichung ist". Bisher habe ich nur über lineare Algebra und nicht über Statistik gesprochen.

Bisher bestand die Idee darin, aus den beobachteten Werten der abhängigen Variablen eine optimale Regressionsgleichung zu finden. Was ist nun, wenn ** die abhängige Variable als stochastische Variable betrachtet, die einer bestimmten Wahrscheinlichkeitsverteilung folgt **? Mit anderen Worten, selbst wenn der Wert der unabhängigen Variablen $ \ boldsymbol {X} $ gleich ist, erhalten Sie bei vielen Beobachtungen jedes Mal einen anderen Beobachtungswert $ \ boldsymbol {y} $. Zu diesem Zeitpunkt ist der partielle Regressionskoeffizient $ \ hat {\ boldsymbol {\ beta}} $, der durch die obige Formel erhalten wird, auch eine stochastische Variable, die sich mit jeder Beobachtung ändert. ** Wir kennen den Wert des wahren partiellen Regressionskoeffizienten nicht, können jedoch die Wahrscheinlichkeitsverteilung berücksichtigen, der der partielle Regressionskoeffizient folgt. ** Die Idee dazu ist dieselbe wie die Intervallschätzung des Bevölkerungsdurchschnitts des statistischen Testgrad-2-Bereichs (der wahre Bevölkerungsdurchschnitt kann nicht bekannt sein, und die Wahrscheinlichkeitsverteilung des Bevölkerungsdurchschnitts wird berechnet).

Wenn Sie beispielsweise tatsächlich eine Regressionsgleichung finden, gibt es möglicherweise unabhängige Variablen, die eine wichtige Rolle in der Regressionsgleichung spielen, und unabhängige Variablen, die nur geringe Auswirkungen haben. Wenn $ \ beta_j = 0 $ ist, hat die entsprechende unabhängige Variable $ x_j $ keine Auswirkung auf den Wert der abhängigen Variablen. Wenn daher die Regressionsgleichung aus dem beobachteten Wert erhalten wird, werden die Nullhypothese $ \ beta_j = 0 $ und die Alternativhypothese $ \ beta_j \ neq 0 $ gesetzt, und wenn die Nullhypothese verworfen wird, ist die unabhängige Variable die Regression. Sie können sich das als einen Ausdruck vorstellen, der die abhängige Variable erklärt. Zu diesem Zweck ist es notwendig, die Wahrscheinlichkeitsverteilung von $ \ hat {\ beta_j} $ zu berücksichtigen, die aus den beobachteten Werten erhalten wird.

Hier nehmen wir an, dass der Fehler $ \ epsilon ^ {(i)} $ jeder Daten einer Normalverteilung mit Mittelwert 0, Varianz $ \ sigma ^ 2 $ unabhängig voneinander folgt. Dann

Das kann man sagen. Der Durchschnitt dieses $ \ hat {\ boldsymbol {\ beta}} $ entspricht dem wahren partiellen Regressionskoeffizienten $ \ boldsymbol {\ beta} $ (Beweis unten). Außerdem ist $ S_ {ij} $ als die $ (i, j) $ -Komponente von $ (\ boldsymbol {X} ^ \ top \ boldsymbol {X}) ^ {-1} $ definiert (jedoch $ \ beta_j $ Oben links ist die (0,0) -Komponente gemäß der Nummerierung), und die Verteilung jedes $ \ hat {\ beta_j} $ ist $ \ sigma ^ 2 S_ {jj} $ (der Beweis wird später beschrieben). Daher wurde $ \ hat {\ beta_j} $ im Durchschnitt auf 0 und bei der Verteilung auf 1 standardisiert.

Z_j = \frac{\hat{\beta_j} - \beta_j}{\sqrt{\sigma^2 S_{jj}}}

Folgt der Standardnormalverteilung, aber wie üblich ist es nicht möglich, den Wert der Fehlervarianz $ \ sigma ^ 2 $ zu kennen. Tatsächlich bleibt keine andere Wahl, als aus dem beobachteten Wert zu schätzen, aber um eine unvoreingenommene Schätzung vorzunehmen, ist der Nenner die Anzahl der Daten minus (unabhängige Variable + 1). Die Anzahl der unabhängigen Variablen sei $ p $

s^2 = \frac{\sum_{i=1}^n (y^{(i)} -  \hat{y}^{(i)})^2}{N-p-1}

Ist eine unvoreingenommene Schätzung von $ \ sigma ^ 2 $, was die geschätzte Fehlervarianz ist. $ \ Sigma ^ 2 $ wurde durch $ s ^ 2 $ ersetzt

t_j = \frac{\hat{\beta_j} - \beta_j}{\sqrt{s^2 S_{jj}}}

Folgt einer $ t $ -Verteilung mit $ n-p-1 $ Freiheitsgraden. Der Grund, warum dies der $ t $ -Verteilung folgt, ähnelt der folgenden "t-Statistik (10.4)". [Statistische Teststufe 2] Zusammenfassung der Wahrscheinlichkeitsverteilungen, die in hypothetischen Tests für normale Populationen auftreten - Qiita

Daher kann beim Testen, ob $ \ beta_j = 0 $ gesagt werden kann,

t_j = \frac{\hat{\beta_j} - 0}{\sqrt{s^2 S_{jj}}} = \frac{\hat{\beta_j}}{\sqrt{(y^{(i)} - \sum_{i=1}^n y^{(i)})^2 S_{jj} / (N-p-1)}}

Berechnen Sie den Wert von und ermitteln Sie den berechneten Wert aus der Verteilungstabelle $ t $ mit $ n-p-1 $ Freiheit, um den Wert $ P $ zu ermitteln. Je größer der absolute Wert des $ t $ -Werts ist, desto kleiner ist der $ P $ -Wert, was bedeutet, dass die unabhängige Variable eine Bedeutung hat.

Auch die zuvor erhaltene geschätzte Fehlervarianz $ s ^ 2 $ stellt den Teil dar, der nicht durch die Regressionsgleichung ausgedrückt werden kann. Man kann also sagen, dass die Regressionsgleichung (ungefähr) umso besser (ungefähr) ist, je kleiner diese ist. Machen wir das.

Die Einführung ist lang geworden, wird aber endlich geübt. Berechnen wir mit den obigen Daten die $ t $ -Statistik und den $ P $ -Wert für die Nullhypothese $ \ beta_j = 0 $ für die partiellen Regressionskoeffizienten, die den beiden unabhängigen Variablen entsprechen. Bei der Berechnung des $ P $ -Werts im Test wird die $ t $ -Verteilungstabelle verwendet, die hier jedoch vom Modul scipy.stats berechnet wird.

from scipy import stats
p = X.shape[1] - 1                               #Anzahl unabhängiger Variablen
S = np.linalg.inv(X.T.dot(X))                    # (X^T X)Inverse Matrix von
s_2 = ((y - X.dot(b)) ** 2).sum() / (n - p - 1)  #Geschätzter Wert der Fehlerstreuung
t = b / np.sqrt(s_2 * S.diagonal())              #t Statistiken
P = (1 - stats.t.cdf(np.abs(t), n - p - 1)) * 2  #P-Wert
print(s_2) # 0.015396963025934803
print(t)   # [  3.273055     3.94856867 -13.93114362]
print(P)   # [0.00207535 0.00028009 0.        ]

** In Statistical Test Level 2 gibt es immer eine Frage, um das Ausgabeergebnis von R zu lesen, aber diese Ausgabe wird durch diese Art der Berechnung erhalten, nicht wahr? ** ** **

Bei Betrachtung dieses Ergebnisses werden sowohl $ \ beta_1 = 0 als auch \ beta_2 = 0 $ mit einem Signifikanzniveau von 5% abgelehnt. Und es kann gelesen werden, dass die Sterblichkeitsrate besser funktioniert als die Geburtenrate. Andererseits wird ab $ s ^ 2 = 0,0154 $ der Standardfehler der Bevölkerungszunahme- / -abnahmerate [%] auf $ s = 0,124 $ geschätzt, aber ich kann nicht sagen, ob dieses Ergebnis gut ist oder nicht, da es kein Vergleichsziel gibt. ..

Ableitung des Mittelwerts des partiellen Regressionskoeffizienten

Lesen Sie unten den erwarteten Wert des Vektors als einen Vektor, in dem die erwarteten Werte jeder Komponente angeordnet sind. Angenommen, der Durchschnitt des Fehlers $ \ epsilon ^ {(i)} $ jeder Daten ist 0

\begin{align}
{\rm E}[\hat{\boldsymbol{\beta}}] &= {\rm E}[(\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y}] \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top {\rm E}[\boldsymbol{y}] \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top {\rm E}[\boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}] \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} \\
&= \boldsymbol{\beta}
\end{align}

Somit entspricht der Durchschnitt der partiellen Regressionskoeffizienten dem wahren partiellen Regressionskoeffizienten von $ \ boldsymbol {\ beta} $ (unverzerrte Schätzung). Es gibt verschiedene Methoden zum Schätzen des partiellen Regressionskoeffizienten, die sich von der Methode des minimalen Quadrats unterscheiden. Das Erhalten einer unverzerrten Schätzung ist jedoch ein großer Vorteil der Methode des minimalen Quadrats.

Ableitung der Varianz des partiellen Regressionskoeffizienten

In Bezug auf diese Seite habe ich die Notation hinzugefügt, um die Zwischenberechnung zu ergänzen. Abschnitt 2.4 Partieller Regressionskoeffizient und sein Test | Datenanalyse / Grundlagen und Anwendungen

Wenn Sie im Folgenden $ {\ rm V} [\ hat {\ boldsymbol {\ beta}}] $ schreiben, wird eine verteilte, gemeinsam verteilte Matrix dargestellt. Das heißt, die $ (i, j) $ -Komponente ist die Kovarianz von $ \ hat {\ beta_i}, \ hat {\ beta_j} $ (die $ (i, j) $ -Komponente ist die Verteilung von $ \ hat {\ beta_i} $ Es ist jedoch eine Matrix, so dass oben links die (0,0) -Komponente gemäß der Nummerierung von $ \ beta_j $) ist.

\begin{align}
\hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y} \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top (\boldsymbol{\epsilon} + \boldsymbol{X} \boldsymbol{\beta}) \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon} + (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon} + \boldsymbol{\beta}
\end{align}

Mit sein

\begin{align}
{\rm V}[\hat{\boldsymbol{\beta}}] &= {\rm E}[(\hat{\boldsymbol{\beta}} - {\rm E}[\hat{\boldsymbol{\beta}}])(\hat{\boldsymbol{\beta}} - {\rm E}[\hat{\boldsymbol{\beta}}])^\top] \\
&= {\rm E}[(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top] \\
&= {\rm E}[((\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon})((\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon})^\top] \\
&= {\rm E}[(\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon} \boldsymbol{\epsilon}^\top \boldsymbol{X} ((\boldsymbol{X}^\top \boldsymbol{X})^{-1})^\top] \\
&= {\rm E}[(\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon} \boldsymbol{\epsilon}^\top \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1}] \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top {\rm E}[\boldsymbol{\epsilon} \boldsymbol{\epsilon}^\top] \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top (\sigma^2 \boldsymbol{I}) \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \\
&= \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \\

\end{align}

($ \ Boldsymbol {I} $ ist eine Einheitsmatrix). Wenn $ S_ {ij} $ als die $ (i, j) $ -Komponente von $ (\ boldsymbol {X} ^ \ top \ boldsymbol {X}) ^ {-1} $ definiert ist, wird jedes $ \ hat { Sie können sehen, dass die Verteilung von \ beta_j} $ $ \ sigma ^ 2 S_ {jj} $ ist.

Untersuchung der Regressionsgleichung (Abschnitte 2-6)

Bisher haben wir zwei unabhängige Variablen verwendet, die Geburtenrate und die Sterblichkeitsrate, um die Bevölkerungsänderungsrate abzuschätzen. Hier werden wir prüfen, ob diese unabhängige Variable wirklich gut ist oder ob es besser ist, andere unabhängige Variablen zu kombinieren. Vielleicht gibt es neben der Geburtenrate und der Sterblichkeitsrate noch andere Faktoren, die die Bevölkerungsänderungsrate beeinflussen.

Wenn Sie willkürlich zwei unabhängige Variablen aus anderen Variablen als der Bevölkerungszunahme- / -abnahmerate (und dem Namen der Präfektur, die die nominale Skala darstellt) auswählen, welche Art von Kombination ist besser?

import itertools
y = df["Bevölkerungswachstums- / -abnahmerate"]
vars = [x for x in df.columns if x not in ["Präfekturen", "Bevölkerungswachstums- / -abnahmerate"]]
for exp_vars in itertools.combinations(vars, 2):     #Wählen Sie zwei unabhängige Variablen
    X = np.c_[np.ones(n), df[list(exp_vars)].values]
    b = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)    #Berechnen Sie den partiellen Regressionskoeffizienten
    TSS = ((y - y.mean()) ** 2).sum()                #Gesamtschwankung
    ESS = ((X.dot(b) - y.mean()) ** 2).sum()         #Regressive Fluktuation
    R_2 = ESS / TSS                                  #Entscheidungskoeffizient
    p = X.shape[1] - 1                               #Anzahl unabhängiger Variablen
    S = np.linalg.inv(X.T.dot(X))                    # (X^T X)Inverse Matrix von
    s_2 = ((y - X.dot(b)) ** 2).sum() / (n - p - 1)  #Geschätzter Wert der Fehlerstreuung
    t = b / np.sqrt(s_2 * S.diagonal())              #t Statistiken
    P = (1 - stats.t.cdf(np.abs(t), n - p - 1)) * 2  #P-Wert
    print("Unabhängige Variable:        ", exp_vars)
    print("Partieller Regressionskoeffizient:      ", b)
    print("Entscheidungskoeffizient:        ", R_2)
    print("t Statistiken:         ", t)
    print("P-Wert:             ", P)
    print("Geschätzter Wert der Fehlerstreuung:", s_2)
    print()

Da es zu lang ist, alle Ergebnisse auszugeben, ist das Folgende ein Auszug der Ausgabe für eine Kombination mit einem Entscheidungskoeffizienten von 0,9 oder mehr. ** (2020/01/07) Die Ergebnisse wurden aufgrund von Datenrevisionen ersetzt. ** ** **

Unabhängige Variable:         ('Population', 'Natürliche Zu- / Abnahmerate')
Partieller Regressionskoeffizient:       [1.35396495e-01 2.92466698e-05 1.38354513e-01]
Entscheidungskoeffizient:         0.9202813093155153
t Statistiken:          [ 2.50498414  4.06232118 16.12219163]
P-Wert:              [0.01602233 0.00019698 0.        ]
Geschätzter Wert der Fehlerstreuung: 0.012452424232620031

Unabhängige Variable:         ('Anzahl der Übermigranten', 'Natürliche Zu- / Abnahmerate')
Partieller Regressionskoeffizient:       [2.36387964e-01 6.30652813e-06 1.43339245e-01]
Entscheidungskoeffizient:         0.9235317925311585
t Statistiken:          [ 6.24974791  4.36740999 18.56543006]
P-Wert:              [1.44848286e-07 7.53394644e-05 0.00000000e+00]
Geschätzter Wert der Fehlerstreuung: 0.011944683881961054

Unabhängige Variable:         ('Geburtenrate', 'Mortalität')
Partieller Regressionskoeffizient:       [ 1.00978055  0.10175178 -0.18246347]
Entscheidungskoeffizient:         0.9014307808651792
t Statistiken:          [  3.273055     3.94856867 -13.93114362]
P-Wert:              [0.00207535 0.00028009 0.        ]
Geschätzter Wert der Fehlerstreuung: 0.015396963025934803

Unabhängige Variable:         ('Geburtenrate', 'Natürliche Zu- / Abnahmerate')
Partieller Regressionskoeffizient:       [ 1.02915847 -0.08316346  0.18254551]
Entscheidungskoeffizient:         0.9030150527461578
t Statistiken:          [ 3.35453348 -2.39398938 14.07003206]
P-Wert:              [0.00164462 0.02099475 0.        ]
Geschätzter Wert der Fehlerstreuung: 0.01514949250939055

Unabhängige Variable:         ('Mortalität', 'Natürliche Zu- / Abnahmerate')
Partieller Regressionskoeffizient:       [ 1.00487249 -0.08024897  0.10178654]
Entscheidungskoeffizient:         0.902445538157104
t Statistiken:          [ 3.30081126 -2.33256484  4.02629601]
P-Wert:              [0.00191782 0.02430887 0.0002203 ]
Geschätzter Wert der Fehlerstreuung: 0.01523845329397937

Unabhängige Variable:         ('Natürliche Zu- / Abnahmerate', 'Heiratsrate')
Partieller Regressionskoeffizient:       [-0.56327258  0.12603311  0.16092148]
Entscheidungskoeffizient:         0.9001397075218234
t Statistiken:          [-1.3456451   7.29460315  2.0734612 ]
P-Wert:              [1.85310560e-01 4.23937485e-09 4.40154160e-02]
Geschätzter Wert der Fehlerstreuung: 0.015598634589388555

Ist es in diesem Fall in Ordnung, weil die Kombination aus der Anzahl der Übermigranten und der natürlichen Zu- / Abnahmerate einen großen Bestimmungskoeffizienten und einen kleinen geschätzten Wert für die Fehlerstreuung aufweist? Die Anzahl der Übermigranten ist jedoch real und nicht prozentual. Daher sollten die Auswirkungen auf die natürliche Zu- oder Abnahmerate je nach ursprünglicher Bevölkerung variieren. Hier sollte stattdessen das, was durch die Bevölkerung geteilt wird, als Übermigrationsrate bezeichnet werden und eine unabhängige Variable sein.

df["Übermigrationsrate"] = df["Anzahl der Übermigranten"] / (df["Population"] * 1000) * 100
n = len(df)
exp_vars = ("Übermigrationsrate", "Natürliche Zu- / Abnahmerate")
y = df["Bevölkerungszunahme / -abnahme"]
X = np.c_[np.ones(n), df[exp_vars].values]
b = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
TSS = ((y - y.mean()) ** 2).sum()                #Gesamtschwankung
ESS = ((X.dot(b) - y.mean()) ** 2).sum()         #Regressive Fluktuation
R_2 = ESS / TSS                                  #Entscheidungskoeffizient
p = X.shape[1] - 1                               #Anzahl unabhängiger Variablen
S = np.linalg.inv(X.T.dot(X))                    # (X^T X)Inverse Matrix von
s_2 = ((y - X.dot(b)) ** 2).sum() / (n - p - 1)  #Geschätzter Wert der Fehlerstreuung
t = b / np.sqrt(s_2 * S.diagonal())              #t Statistiken
P = (1 - stats.t.cdf(np.abs(t), n - p - 1)) * 2  #P-Wert
print("Unabhängige Variable:        ", exp_vars)
print("Partieller Regressionskoeffizient:      ", b)
print("Entscheidungskoeffizient:        ", R_2)
print("t Statistiken:         ", t)
print("P-Wert:             ", P)
print("Geschätzter Wert der Fehlerstreuung:", s_2)
print()
Unabhängige Variable:         ('Übermigrationsrate', 'Natürliche Zu- / Abnahmerate')
Partieller Regressionskoeffizient:       [0.20976966 0.75632976 0.10990394]
Entscheidungskoeffizient:         0.9580723972439072
t Statistiken:          [ 7.50186534  8.42827532 14.31096232]
P-Wert:              [2.11598050e-09 9.87150361e-11 0.00000000e+00]
Geschätzter Wert der Fehlerstreuung: 0.006549283387531298

Im Vergleich zur Verwendung der überschüssigen Anzahl von Zuwanderern sind der Bestimmungskoeffizient und der geschätzte Wert der Fehlerstreuung besser. Es ist zu beachten, dass sowohl die Übermigrationsrate als auch die natürliche Zu- / Abnahmerate die Bevölkerung erhöhen, so dass der entsprechende partielle Regressionskoeffizient sicherlich ein positiver Wert ist.

Was ist mit drei unabhängigen Variablen?

for exp_vars in itertools.combinations(vars, 3):     #Wählen Sie 3 unabhängige Variablen
    X = np.c_[np.ones(n), df[list(exp_vars)].values]
    b = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)    #Berechnen Sie den partiellen Regressionskoeffizienten
    TSS = ((y - y.mean()) ** 2).sum()                #Gesamtschwankung
    ESS = ((X.dot(b) - y.mean()) ** 2).sum()         #Regressive Fluktuation
    R_2 = ESS / TSS                                  #Entscheidungskoeffizient
    p = X.shape[1] - 1                               #Anzahl unabhängiger Variablen
    S = np.linalg.inv(X.T.dot(X))                    # (X^T X)Inverse Matrix von
    s_2 = ((y - X.dot(b)) ** 2).sum() / (n - p - 1)  #Geschätzter Wert der Fehlerstreuung
    t = b / np.sqrt(s_2 * S.diagonal())              #t Statistiken
    P = (1 - stats.t.cdf(np.abs(t), n - p - 1)) * 2  #P-Wert
    print("Unabhängige Variable:        ", exp_vars)
    print("Partieller Regressionskoeffizient:      ", b)
    print("Entscheidungskoeffizient:        ", R_2)
    print("t Statistiken:         ", t)
    print("P-Wert:             ", P)
    print("Geschätzter Wert der Fehlerstreuung:", s_2)
    print()

Zum Beispiel können Sie eine solche Kombination finden. (Dies ist nicht das Beste. Es ist ein Beispiel, das die Interpretation der Ergebnisse erleichtert.)

Unabhängige Variable:         ('Natürliche Zu- / Abnahmerate', 'Heiratsrate', 'Scheidungsrate')
Partieller Regressionskoeffizient:       [-0.30054533  0.12984192  0.19858489 -0.25130329]
Entscheidungskoeffizient:         0.9109243128862279
t Statistiken:          [-0.72219064  7.82613187  2.61426129 -2.28169067]
P-Wert:              [4.74086687e-01 8.36975600e-10 1.22779282e-02 2.75178580e-02]
Geschätzter Wert der Fehlerstreuung: 0.014237611977607154

Je höher die natürliche Rate der Zunahme / Abnahme und die Rate der Ehe ist, desto höher ist die Bevölkerung und je höher die Scheidungsrate, desto niedriger die Bevölkerung. Es ist ein Ergebnis, das zu überzeugen scheint. Alle drei partiellen Regressionskoeffizienten (mit Ausnahme des konstanten Terms) sind 5% signifikant. Schauen wir uns auf der anderen Seite einige nicht so gute Beispiele an.

Unabhängige Variable:         ('Geburtenrate', 'Mortalität', 'Natürliche Zu- / Abnahmerate')
Partieller Regressionskoeffizient:       [ 1.03883742 -0.2559058   0.17133451  0.35353742]
Entscheidungskoeffizient:         0.9034996234109349
t Statistiken:          [ 3.3482056  -0.68534239  0.4646739   0.96013752]
P-Wert:              [0.00169918 0.49680614 0.6445094  0.34235404]
Geschätzter Wert der Fehlerstreuung: 0.015424353851135318

Das seltsame Ergebnis ist, dass mit steigender Geburtenrate die Bevölkerung abnimmt und mit steigender Sterblichkeitsrate die Bevölkerung zunimmt. Dies liegt daran, dass zwischen den unabhängigen Variablen "natürliche Anstiegs- / Abnahmerate = Geburtenrate-Sterblichkeitsrate" eine Beziehung besteht, die das problematische Problem der multiplen Co-Linearität verursacht. Was ist multiple Co-Linearität? Weblio-Wörterbuch

Wenn keine Variablenauswahl durchgeführt wird, ist es besser, diejenigen mit hoher Korrelation zwischen unabhängigen Variablen nicht einzubeziehen. Wenn einer von ihnen nicht unabhängig ist (z. B. sind die Variablen A und B und ihre Summe C = A + B enthalten), schlägt die Analyse fehl. In einigen Fällen stimmt das Vorzeichen des Korrelationskoeffizienten zwischen jeder unabhängigen Variablen und der abhängigen Variablen möglicherweise nicht mit dem Vorzeichen des partiellen Regressionskoeffizienten überein. (Weggelassen) Dies geschieht, weil es eine Mischung aus stark korrelierten unabhängigen Variablen gibt (wenn eine Variable den überprognostizierten Teil mit einer anderen Variablen aufhebt). Es gibt). Da dies jedoch bei der Betrachtung des Kausalzusammenhangs unpraktisch ist, ist es ratsam, nach einer multiplen Regressionsgleichung zu suchen, die unabhängige Variablen ausschließt, deren Vorzeichen nicht übereinstimmen.

Es kann möglich sein, festzustellen, ob bis zu einem gewissen Grad eine multiple Co-Linearität auftritt oder nicht, aber es ist für Maschinen schwierig, die Geschichte zu verstehen, dass "je höher die Geburtenrate, desto mehr wird die Bevölkerung nicht zunehmen". Es scheint, dass "multiple Regressionsgleichungen ohne unabhängige Variablen, deren Vorzeichen nicht übereinstimmen", tatsächlich ziemlich schwierig sind. Die Interpretation des Endergebnisses ist dort, wo Menschen denken müssen.

Damit

Die lange Regressionsanalyse ist vorerst beendet. Nächstes Mal planen wir Kapitel 3 "Hauptkomponentenanalyse".

Recommended Posts

[Statistischer Test 2. Klasse / quasi 1. Klasse] Regressives Analysetraining mit Python (2)
[Statistischer Test 2. Klasse / quasi 1. Klasse] Regressives Analysetraining mit Python (1)
Regressionsanalyse mit Python
In Python ② erlernte statistische Wahrscheinlichkeitsverteilung für Testgrad 2
In Python ① erlernte statistische Wahrscheinlichkeitsverteilung für Testgrad 2
Einfache Regressionsanalyse mit Python
Erste einfache Regressionsanalyse in Python
Statistischer Test (Mehrfachtest) in Python: scikit_posthocs
[Statistische Teststufe 2] Diskrete Wahrscheinlichkeitsverteilung
2. Multivariate Analyse in Python 1-1. Einfache Regressionsanalyse (Scikit-Learn)
2. Multivariate Analyse in Python 7-3. Entscheidungsbaum [Rückgabebaum]
2. Multivariate Analyse in Python 2-1. Multiple Regressionsanalyse (Scikit-Learn)
2. Multivariate Analyse in Python 1-2. Einfache Regressionsanalyse (Algorithmus)
Assoziationsanalyse in Python
2. Multivariate Analyse in Python 5-3. Logistische Regressionsanalyse (Statistikmodelle)
Mehrfacher Regressionsausdruck in Python
Algorithmus in Python (Haupturteil)
Python Statistical Techniques-Statistische Analyse gegen Python-
Axialsymmetrische Spannungsanalyse mit Python
Online lineare Regression in Python
Stellen Sie den Python-Test in Jenkins ein
2. Multivariate Analyse in Python 6-2. Ridge-Regression / Lasso-Regression (Scikit-Learn) [Ridge-Regression vs. Lasso-Regression]
2. Multivariate Analyse in Python 2-3. Multiple Regressionsanalyse [COVID-19-Infektionsrate]
Gehirnwellenanalyse mit Python: Python MNE-Tutorial
Schreiben Sie Selentestcode in Python
Planare Skelettanalyse in Python (2) Hotfix
Einfache Implementierung einer Regressionsanalyse mit Keras
Logistische Regressionsanalyse Selbst erstellt mit Python
2. Multivariate Analyse in Python 6-1. Ridge-Regression / Lasso-Regression (Scikit-Learn) [multiple Regression vs. Ridge-Regression]
2. Multivariate Analyse in Python 8-2. K Nachbarschaftsmethode [Gewichtungsmethode] [Rückgabemodell]