[PYTHON] Einführung in verallgemeinerte Schätzungsgleichungen durch Statistikmodelle
Statsmodels ist eine statistische Analysesoftware, die auf einer Programmiersprache namens Python ausgeführt wird. Python muss auf Ihrem PC installiert sein, damit das Beispiel für die Statistikmodelle ausgeführt werden kann. Wenn Sie es noch nicht installiert haben, lesen Sie bitte Jupyter-Notebook installieren. Jupyter Notebook ist sehr nützlich für die Ausführung von Statistikmodellen.
Das meiste davon basiert auf der statsmodels-Referenz.
Verallgemeinerte Schätzgleichungen
Die verallgemeinerte Schätzungsgleichung ist ein Modell, das darauf abzielt, den durchschnittlichen Effekt von beobachteten Werten wie Panels, Clustern und iterativen Messdaten in einer Datenpopulation zu schätzen, die innerhalb des Clusters korreliert und außerhalb des Clusters nicht korreliert ist. , Verallgemeinertes lineares Modell. Obwohl gesagt wird, dass es eine Korrelation innerhalb des Clusters gibt, kann gesagt werden, dass die Ursache der Korrelation nicht identifiziert wurde oder nicht identifiziert werden kann. Die Durchschnittsstruktur und die Korrelationsstruktur können getrennt behandelt werden. GLM
- Die abhängige Variable $ y_i $ folgt der Exponentialverteilungsfamilie.
- Jede abhängige Variable ist unabhängig.
Unter der Annahme
- Annahme der Wahrscheinlichkeitsverteilung der Wahrscheinlichkeitsvariablen $ y_i $. -Stellen Sie ein lineares Regressionsmodell mit dem erwarteten Wert von $ y_i $ $ \ mu $ ein. --Erstellen Sie eine Wahrscheinlichkeitsfunktion für die angenommene Wahrscheinlichkeitsverteilung von $ y_i $ und suchen Sie $ \ beta (\ hat {\ beta}) $, um sie zu maximieren.
- Die wahrscheinlichste Schätzung $ \ hat {\ beta} $ hat Konsistenz und nahezu Normalität.
Schätzen Sie das Modell als.
In GEE
- Entspannen Sie die Annahme, dass die abhängige Variable $ y_i $ der Exponentialverteilungsfamilie folgt.
- Entspannen Sie die Annahme, dass die abhängige Variable unabhängig ist.
- Entspricht einer unbekannten Verteilung unter Verwendung der Pseudowahrscheinlichkeit.
- Führen Sie die funktionierende Korrelationsmatrix $ R_i (\ alpha) $ ein, um die Korrelationsstruktur darzustellen.
- Die verallgemeinerte Schätzungsgleichung ist in zwei Stufen unterteilt und wird durch numerische Berechnung gelöst.
Schätzen Sie das Modell in dieser Reihenfolge. Daher kann es mit der Vielfalt der Daten umgehen, benötigt aber auch die Menge an Informationen.
Beispiel: Anzahl der Anfälle bei Epilepsiepatienten
Thall und Vail (1990) analysierten den 8-wöchigen Basiszeitraum und 4 aufeinanderfolgende zweiwöchentliche Anfälle bei 59 Patienten mit Epilepsie. Während des Basiszeitraums wird dem Patienten keine Behandlung gegeben. Nach dem Basiszeitraum wird den Patienten zufällig Placebo oder Progavid verschrieben. Beobachtungen nach der Verschreibung werden für 4x2 = 8 Wochen durchgeführt.
Datendetails: Verwendung: Epilepsie Format: 236 Zeilen und 9 Spalten
y: Anzahl der Angriffe in 2 Wochen trt: Behandlung, Placebo oder Progavid Basis: Anzahl der Angriffe während des 8-wöchigen Basiszeitraums Alter: Alter des Probanden, Alter V4: 0/1 Indikatorvariable für 4 Perioden Betreff: Betreffnummer 1 bis 59. Zeitraum: Zeitraum 1 bis 4. lbase: Standardisieren Sie den Logarithmus und den Durchschnitt der Anzahl der Angriffe während des Basiszeitraums auf Null lage: Altersprotokoll, standardisiert auf Nulldurchschnitt
Epileptische Anfälle entsprechen wiederkehrenden Ereignissen in der Überlebenszeitanalyse. Wiederholungsereignisdaten werden durch wiederholtes Messen des Auftretens von Angriffen in demselben Subjekt erhalten und können als eine Art von Längsschnittdaten und wiederholten Messdaten angesehen werden. Es wird angenommen, dass es keine Korrelation in der Anzahl der Epilepsieanfälle bei verschiedenen Probanden gibt, aber es ist natürlich anzunehmen, dass es eine Korrelation in den Ergebnissen gibt, wenn die Messungen für dasselbe Subjekt durchgeführt werden.
Insgesamt 4 Beobachtungen wurden in Intervallen von 8 Wochen vor der Arzneimittelverabreichung und 2 Wochen nach der Arzneimittelverabreichung gemacht, und die Anzahl der Anfälle wurde für jeden Abschnitt untersucht. Nur das Alter des Subjekts ist eine Kovariate.
Initialisieren Sie zunächst, um die Daten abzurufen.
import statsmodels.formula.api as smf
data = sm.datasets.get_rdataset('epil', package='MASS').data
Beschreibende Statistik
Untersuchen Sie den Status der Daten.
data.head(10)
y trt base age V4 subject period lbase lage
0 5 placebo 11 31 0 1 1 -0.756354 0.114204
1 3 placebo 11 31 0 1 2 -0.756354 0.114204
2 3 placebo 11 31 0 1 3 -0.756354 0.114204
3 3 placebo 11 31 1 1 4 -0.756354 0.114204
4 3 placebo 11 30 0 2 1 -0.756354 0.081414
5 5 placebo 11 30 0 2 2 -0.756354 0.081414
6 3 placebo 11 30 0 2 3 -0.756354 0.081414
7 3 placebo 11 30 1 2 4 -0.756354 0.081414
8 2 placebo 6 25 0 3 1 -1.362490 -0.100908
9 4 placebo 6 25 0 3 2 -1.362490 -0.100908
Vergleichen wir die Verteilung der Anzahl epileptischer Anfälle, nachdem ein Medikament verschrieben wurde, das nicht den Wirkstoff Placebo enthält, und ein GABA-ähnliches Antiepileptikum namens Progavid. Die Anzahl der Daten beträgt 112 für die erstere und 124 für die letztere.
data[data.trt!="placebo"].y.hist()
data[data.trt=="placebo"].y.hist()
data[data.trt!="placebo"].y.count(),data[data.trt=="placebo"].y.count()
(124, 112)
Lassen Sie uns ein Frequenzdiagramm der Anzahl der Angriffe und des Alters des standardisierten Basiszeitraums erstellen.
data.lbase.hist()
data.lage.hist()
Als nächstes erhalten wir beschreibende Statistiken. Durchschnittliche und Standardabweichung der Anzahl der Angriffe während des Basiszeitraums
data.base.mean(),data.base.std()
(31.220338983050848, 26.705051109822254)
Durchschnitt und Standardabweichung der Anzahl der Angriffe nach Progavid
data[data.trt!="placebo"].y.mean(),data[data.trt!="placebo"].y.std()
(7.959677419354839, 13.92978863002629)
Durchschnittliche und Standardabweichung der Anzahl der Anfälle während des Basiszeitraums von Patienten, die Progavid erhielten
data[data.trt!="placebo"].base.mean(),data[data.trt!="placebo"].base.std()
(31.612903225806452, 27.638405492528324)
Durchschnitt und Standardabweichung der Anzahl der Anfälle nach Anwendung von Placebo
data[data.trt=="placebo"].y.mean(),data[data.trt=="placebo"].y.std()
(8.580357142857142, 10.369426989352696)
Durchschnittliche und Standardabweichung der Anzahl der Anfälle während des Basiszeitraums von Patienten, die Placebo erhielten
data[data.trt=="placebo"].base.mean(),data[data.trt=="placebo"].base.std()
(30.785714285714285, 25.749111266541444)
Das Protokoll des Alters des Probanden und der Anzahl der Anfälle während des Basiszeitraums wird aufgezeichnet.
plt.scatter(data.lage,data.lbase)
rg=sm.OLS(data.lbase,data.lage).fit()
plt.plot(data.lage,rg.fittedvalues)
Ratet mal Statistiken
Schätzen Sie das Modell. Es wird versucht, die Poisson-Verteilung für die Anzahl der Angriffe zu verfolgen.
Unterschied in der Arbeitskorrelationsmatrix
independence()
Für die Arbeitskorrelationsmatrix wird angenommen, dass die gepaarten Daten der wiederholten Messungen jedes Patienten unabhängig sind. "Betreff" kennzeichnet Gruppen.
fam = sm.families.Poisson()
ind = sm.cov_struct.Independence()
mod = smf.gee("y ~ age + trt + base", "subject", data,
cov_struct=ind, family=fam)
res = mod.fit()
print(res.summary())
GEE Regression Results
===================================================================================
Dep. Variable: y No. Observations: 236
Model: GEE No. clusters: 59
Method: Generalized Min. cluster size: 4
Estimating Equations Max. cluster size: 4
Family: Poisson Mean cluster size: 4.0
Dependence structure: Independence Num. iterations: 2
Date: Sat, 02 Nov 2019 Scale: 1.000
Covariance type: robust Time: 13:06:54
====================================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------------
Intercept 0.5730 0.361 1.589 0.112 -0.134 1.280
trt[T.progabide] -0.1519 0.171 -0.888 0.375 -0.487 0.183
age 0.0223 0.011 1.960 0.050 2.11e-06 0.045
base 0.0226 0.001 18.451 0.000 0.020 0.025
==============================================================================
Skew: 3.7823 Kurtosis: 28.6672
Centered skew: 2.7597 Centered kurtosis: 21.9865
==============================================================================
Nehmen Sie als nächstes den Durchschnittswert der Anzahl der Angriffe für jedes Subjekt und den Logarithmus der Varianz und tupfen Sie ihn ab.
yg = mod.cluster_list(np.asarray(data["y"]))
ymn = [x.mean() for x in yg]
yva = [x.var() for x in yg]
plt.grid(True)
plt.plot(np.log(ymn), np.log(yva), 'o')
plt.plot([-2, 6], [-2, 6], '-', color='grey')
plt.xlabel("Log variance", size=17)
plt.ylabel("Log mean", size=17)
exchangeable() In diesem Fall wird angenommen, dass die gepaarten Daten der wiederholten Messungen jedes Patienten für die Arbeitskorrelationsmatrix dieselbe Korrelation aufweisen.
fam = sm.families.Poisson()
ex = sm.cov_struct.Exchangeable()
mod = smf.gee("y ~ age + trt + base", "subject", data,
cov_struct=ex, family=fam)
res = mod.fit()
print(res.summary())
GEE Regression Results
===================================================================================
Dep. Variable: y No. Observations: 236
Model: GEE No. clusters: 59
Method: Generalized Min. cluster size: 4
Estimating Equations Max. cluster size: 4
Family: Poisson Mean cluster size: 4.0
Dependence structure: Exchangeable Num. iterations: 2
Date: Thu, 31 Oct 2019 Scale: 1.000
Covariance type: robust Time: 07:23:48
====================================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------------
Intercept 0.5730 0.361 1.589 0.112 -0.134 1.280
trt[T.progabide] -0.1519 0.171 -0.888 0.375 -0.487 0.183
age 0.0223 0.011 1.960 0.050 2.11e-06 0.045
base 0.0226 0.001 18.451 0.000 0.020 0.025
==============================================================================
Skew: 3.7823 Kurtosis: 28.6672
Centered skew: 2.7597 Centered kurtosis: 21.9865
==============================================================================
fig = res.plot_isotropic_dependence()
fig.get_axes()[0].set_ylim(0, 1)
rslts = res.params_sensitivity(0., 0.9, 5)
dp = [x.cov_struct.dep_params for x in rslts]
age_params = np.asarray([x.params[2] for x in rslts])
age_se = np.asarray([x.bse[2] for x in rslts])
plt.grid(True)
plt.plot(dp, age_params, '-', color='orange', lw=4)
plt.plot(dp, age_params + age_se, '-', color='grey')
plt.plot(dp, age_params - age_se, '-', color='grey')
plt.axvline(res.cov_struct.dep_params, ls='--', color='black')
plt.xlabel("Autocorrelation parameter", size=14)
plt.ylabel("Age effect", size=14)
plt.plot(res.fittedvalues, res.resid, 'o', alpha=0.5)
plt.xlabel("Fitted values", size=17)
plt.ylabel("Residuals", size=17)
Autoregressive() In diesem Fall wird angenommen, dass die gepaarten Daten der wiederholten Messungen jedes Patienten für die Arbeitskorrelationsmatrix autokorreliert sind.
fam = sm.families.Poisson()
au = sm.cov_struct.Autoregressive()
mod = smf.gee("y ~ age + trt + base", "subject", data,
cov_struct=au, family=fam)
res = mod.fit()
print(res.summary())
GEE Regression Results
===================================================================================
Dep. Variable: y No. Observations: 236
Model: GEE No. clusters: 59
Method: Generalized Min. cluster size: 4
Estimating Equations Max. cluster size: 4
Family: Poisson Mean cluster size: 4.0
Dependence structure: Autoregressive Num. iterations: 11
Date: Sat, 02 Nov 2019 Scale: 1.000
Covariance type: robust Time: 13:08:52
====================================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------------
Intercept 0.4123 0.422 0.977 0.329 -0.415 1.240
trt[T.progabide] -0.0704 0.164 -0.430 0.667 -0.391 0.250
age 0.0274 0.013 2.108 0.035 0.002 0.053
base 0.0223 0.001 18.252 0.000 0.020 0.025
==============================================================================
Skew: 3.9903 Kurtosis: 30.1753
Centered skew: 2.7597 Centered kurtosis: 21.9865
==============================================================================
fig = res.plot_isotropic_dependence()
fig.get_axes()[0].set_ylim(0, 1)
rslts = res.params_sensitivity(0., 0.9, 5)
dp = [x.cov_struct.dep_params for x in rslts]
age_params = np.asarray([x.params[2] for x in rslts])
age_se = np.asarray([x.bse[2] for x in rslts])
plt.grid(True)
plt.plot(dp, age_params, '-', color='orange', lw=4)
plt.plot(dp, age_params + age_se, '-', color='grey')
plt.plot(dp, age_params - age_se, '-', color='grey')
plt.axvline(res.cov_struct.dep_params, ls='--', color='black')
plt.xlabel("Autocorrelation parameter", size=14)
plt.ylabel("Age effect", size=14)
Ob die Verabreichung des Studienmedikaments die Anzahl der Anfälle im Vergleich zu Placebo verringerte oder nicht, wurde durch Ändern der Arbeitskorrelationsmatrix getestet. Das Ergebnis ist jedoch signifikant, wenn die verallgemeinerte Schätzungsgleichung aus dieser klinischen Studie verwendet wird Ich kann es nicht verstehen
Beispiel: Epilepsie (Leppik et al. 1987)
Schauen wir uns den Fall an, in dem die Ergebnisse klinischer Studien anhand der allgemeinen Schätzungsgleichung anhand der Daten von Epilepsieanfällen derselben 59 Personen bewertet werden.
Ändern Sie zunächst die Struktur der Daten und geben Sie die Anzahl der Angriffe während des Basiszeitraums in $ y_ {ij} $ ein.
import pandas as pd
yy=pd.DataFrame(np.arange(59)+1,columns=['subject'])
yy['period']=np.log(8)
yy['y']=np.NaN
yy['trt']=0#np.NaN
yy['btrt']=0
data['btrt']=1
data['period']=np.log(2)
ii=0
for i in range(len(data)):
if data.subject.iloc[i]==yy.subject.iloc[ii]:
yy.iloc[ii,2]=data.base.iloc[i]
yy.iloc[ii,3]=data.trt.iloc[i]
yy.iloc[ii,0]=data.subject.iloc[i]
ii+=1
if ii>=59:
break
yy=pd.concat([data,yy],axis=0,join='inner')
for i in range(len(yy)):
if yy.iloc[i,1]=='progabide':
yy.iloc[i:i+1,1]=1
if yy.iloc[i,1]=='placebo':
yy.iloc[i:i+1,1]=0
mod = smf.gee("y ~ btrt + trt ", "subject", yy,
cov_struct=ex, family=fam,offset='period')
res = mod.fit()
print(res.summary())
GEE Regression Results
===================================================================================
Dep. Variable: y No. Observations: 295
Model: GEE No. clusters: 59
Method: Generalized Min. cluster size: 5
Estimating Equations Max. cluster size: 5
Family: Poisson Mean cluster size: 5.0
Dependence structure: Exchangeable Num. iterations: 6
Date: Sat, 02 Nov 2019 Scale: 1.000
Covariance type: robust Time: 23:30:28
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 1.3240 0.168 7.883 0.000 0.995 1.653
btrt 0.0560 0.106 0.530 0.596 -0.151 0.263
trt 0.0705 0.213 0.331 0.740 -0.346 0.487
==============================================================================
Skew: 3.3538 Kurtosis: 15.0311
Centered skew: 2.0882 Centered kurtosis: 6.6169
==============================================================================
Es werden keine Ergebnisse erhalten, dass die Studiengruppe der klinischen Studien signifikant ist. Dort sieht es etwas anders aus. Fügen Sie das Alter des Subjekts hinzu, das eine Kovariate ist.
import pandas as pd
yy=pd.DataFrame(np.arange(59)+1,columns=['subject'])
yy['period']=0
yy['y']=np.NaN
yy['trt']=0#np.NaN
yy['btrt']=0
yy['age']=0
data = sm.datasets.get_rdataset('epil', package='MASS').data
data['btrt']=1
ii=0
for i in range(len(data)):
if data.subject.iloc[i]==yy.subject.iloc[ii]:
yy.iloc[ii,2]=data.base.iloc[i]
#yy.iloc[ii,3]=data.trt.iloc[i]
yy.iloc[ii,0]=data.subject.iloc[i]
yy.iloc[ii,5]=data.age.iloc[i]
ii+=1
if ii>=59:
break
#print(yy.head(200))
yy=pd.concat([data,yy],axis=0,join='inner')
for i in range(len(yy)):
if yy.iloc[i,1]=='progabide':
yy.iloc[i,1]=1
if yy.iloc[i,1]=='placebo':
yy.iloc[i,1]=0
yy[yy.trt==0].count()
mod = smf.gee("y ~ age + btrt + trt ", "subject", yy,
cov_struct=ind, family=fam,offset='period')
res = mod.fit()
print(res.summary())
GEE Regression Results
===================================================================================
Dep. Variable: y No. Observations: 295
Model: GEE No. clusters: 59
Method: Generalized Min. cluster size: 5
Estimating Equations Max. cluster size: 5
Family: Poisson Mean cluster size: 5.0
Dependence structure: Independence Num. iterations: 2
Date: Sat, 02 Nov 2019 Scale: 1.000
Covariance type: robust Time: 23:48:30
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 3.9937 0.614 6.503 0.000 2.790 5.197
age -0.0198 0.021 -0.954 0.340 -0.060 0.021
btrt -4.3316 0.154 -28.050 0.000 -4.634 -4.029
trt -0.1014 0.335 -0.303 0.762 -0.758 0.555
==============================================================================
Skew: 2.6891 Kurtosis: 13.3017
Centered skew: 0.3066 Centered kurtosis: 3.7456
==============================================================================
Der Punkt, dass es nicht signifikant wird, bleibt unverändert.
Verweise:
Analyse von Wiederholungsereignisdaten Wiki notebooks for GEE repeated measures of epileptic seizure counts
Recommended Posts