[PYTHON] Introduction aux équations d'estimation généralisées par des modèles de statistiques

Statsmodels est un logiciel d'analyse statistique qui s'exécute sur un langage de programmation appelé Python. Python doit être installé sur votre PC pour exécuter l'exemple statsmodels. Si vous ne l'avez pas encore installé, reportez-vous à Installation du notebook Jupyter. Le notebook Jupyter est très utile pour exécuter des modèles de statistiques.

La plupart sont basés sur la référence statsmodels.

Équations d'estimation généralisées

L'équation d'estimation généralisée est un modèle qui vise à estimer l'effet moyen des valeurs observées des panels, des grappes et des données de mesure itératives dans une population de données corrélées au sein du cluster et non corrélées en dehors du cluster. , Modèle linéaire généralisé. Bien que l'on dise qu'il existe une corrélation au sein du cluster, on peut dire que la cause de la corrélation n'a pas été identifiée ou ne peut pas être identifiée. La structure moyenne et la structure de corrélation peuvent être traitées séparément. GLM

Dans l'hypothèse

Estimez le modèle comme.

Dans GEE

--Relaxez l'hypothèse que la variable dépendante $ y_i $ suit la famille de distribution exponentielle.

Estimez le modèle dans cet ordre. Par conséquent, il est possible de gérer la diversité des données, mais cela nécessite également la quantité d'informations.

Exemple: nombre de crises chez les patients épileptiques

Thall et Vail (1990) ont analysé la période de référence de 8 semaines et 4 crises consécutives aux deux semaines chez 59 patients épileptiques. Aucun traitement n'est administré au patient pendant la période de référence. Après la période de référence, les patients reçoivent au hasard un placebo ou un progavide. Les observations post-prescription sont effectuées pendant 4x2 = 8 semaines.

Détails des données: Utilisation: épilepsie Format: 236 lignes et 9 colonnes

y: Nombre d'attaques en 2 semaines trt: traitement, placebo ou progavide base: nombre d'attaques au cours de la période de référence de 8 semaines âge: âge, âge du sujet V4: variable indicatrice 0/1 pour 4 périodes sujet: numéro du sujet, 1 à 59. période: période, 1 à 4. lbase: normaliser le logarithmique et la moyenne du nombre d'attaques pendant la période de référence à zéro lage: Journal des âges, normalisé à une moyenne nulle

Les crises épileptiques correspondent à des événements récurrents dans l'analyse du temps de survie. Les données d'événement de récurrence sont obtenues en mesurant à plusieurs reprises la survenue de crises chez le même sujet, et peuvent être considérées comme un type de données longitudinales et de données de mesures répétées. On considère qu'il n'y a pas de corrélation dans le nombre de crises d'épilepsie chez différents sujets, mais il est naturel de penser qu'il y a une corrélation dans les résultats si les mesures sont faites pour le même sujet.

Un total de 4 observations ont été faites à des intervalles de 8 semaines avant l'administration du médicament et 2 semaines après l'administration du médicament, et le nombre d'attaques a été examiné pour chaque section. Seul l'âge du sujet est une covariable.

Tout d'abord, initialisez pour obtenir les données.

import statsmodels.formula.api as smf
data = sm.datasets.get_rdataset('epil', package='MASS').data

Statistiques descriptives

Examinez l'état des données.

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

Comparons la distribution du nombre de crises d'épilepsie après s'être vu prescrire un médicament ne contenant pas l'ingrédient actif appelé placebo et un antiépileptique de type GABA appelé progavide. Le nombre de données est de 112 pour les premiers et de 124 pour les seconds.

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)

image.png

Faisons un graphique de fréquence du nombre d'attaques et du nombre d'âges au cours de la période de référence standardisée.

data.lbase.hist()

image.png

data.lage.hist()

image.png

Ensuite, nous obtenons des statistiques descriptives. Écart moyen et standard du nombre d'attaques au cours de la période de référence

data.base.mean(),data.base.std()
(31.220338983050848, 26.705051109822254)

Écart moyen et standard du nombre d'attaques après progavide

data[data.trt!="placebo"].y.mean(),data[data.trt!="placebo"].y.std()
(7.959677419354839, 13.92978863002629)

Écart moyen et standard du nombre de crises pendant la période de référence des patients ayant reçu Progavid

data[data.trt!="placebo"].base.mean(),data[data.trt!="placebo"].base.std()
(31.612903225806452, 27.638405492528324)

Écart moyen et standard du nombre d'attaques après application du placebo

data[data.trt=="placebo"].y.mean(),data[data.trt=="placebo"].y.std()
(8.580357142857142, 10.369426989352696)

Écart moyen et standard du nombre d'attaques au cours de la période de référence des patients ayant reçu le placebo

data[data.trt=="placebo"].base.mean(),data[data.trt=="placebo"].base.std()
(30.785714285714285, 25.749111266541444)

Le journal de l'âge du sujet et du nombre d'attaques au cours de la période de référence est tracé.

plt.scatter(data.lage,data.lbase)
rg=sm.OLS(data.lbase,data.lage).fit()
plt.plot(data.lage,rg.fittedvalues)

image.png

Devinez les statistiques

Estimez le modèle. Il essaie de suivre la distribution de Poisson pour le nombre d'attaques.

Différence dans la matrice de corrélation de travail

independence()

Pour la matrice de corrélation de travail, on suppose que les données appariées des mesures répétées de chaque patient sont indépendantes. «sujet» indique des groupes.

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
==============================================================================

Ensuite, prenez la valeur moyenne du nombre d'attaques pour chaque sujet et le logarithme de la variance et effacez-la.

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)

image.png

exchangeable() Dans ce cas, on suppose que les données appariées des mesures répétées de chaque patient pour la matrice de corrélation de travail ont la même corrélation.

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)

image.png

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)

image.png

plt.plot(res.fittedvalues, res.resid, 'o', alpha=0.5)
plt.xlabel("Fitted values", size=17)
plt.ylabel("Residuals", size=17)

image.png

Autoregressive() Dans ce cas, on suppose que les données appariées des mesures répétées de chaque patient pour la matrice de corrélation de travail sont autocorrélées.

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)

image.png

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)

image.png

La question de savoir si l'administration du médicament à l'étude réduisait ou non le nombre d'attaques par rapport au placebo a été testée en modifiant la matrice de corrélation de travail, mais le résultat selon lequel elle est significative lors de l'utilisation de l'équation d'estimation généralisée de cette étude clinique est Je ne peux pas comprendre.

Exemple: l'épilepsie (Leppik et.al.1987)

Regardons le cas où les résultats des essais cliniques sont évalués à l'aide de l'équation d'estimation générale en utilisant les données des crises d'épilepsie des mêmes 59 personnes.

Tout d'abord, changez la structure des données et mettez le nombre d'attaques pendant la période de référence dans $ y_ {ij} . $log E(Y_{ij})=log(T_{ij})+\beta_0+\beta_1x_{ij1}+\beta_2x_{ij2}+\beta_3x_{ij1}x_{ij2}$$ Ici, $ x_ {ij1} $ exprime la période avant (0) et après (1) de la randomisation sous forme de valeurs binaires. $ x_ {ij2} $ indique le groupe de traitement du placebo (0) et du médicament à l'étude (1) en binaire. $ T_ {ij} $ est la durée de la période de référence de 8 semaines et l'intervalle de 2 semaines pour chaque période de mesure post-traitement.

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
==============================================================================

Aucun résultat n'est obtenu indiquant que le groupe d'étude des essais cliniques est significatif. C'est un peu différent là-bas. Ajoutez l'âge du sujet, qui est une covariable.

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
==============================================================================

Le fait qu'il ne devienne pas significatif reste inchangé.

Les références:

Analyse des données d'événements de récurrence Wiki notebooks for GEE repeated measures of epileptic seizure counts

Recommended Posts

Introduction aux équations d'estimation généralisées par des modèles de statistiques
Introduction à MQTT (Introduction)
Introduction à Scrapy (1)
Introduction à Scrapy (3)
Premiers pas avec Supervisor
Introduction à Tkinter 1: Introduction
Introduction à PyQt
Introduction à Scrapy (2)
[Linux] Introduction à Linux
Introduction à Scrapy (4)
Introduction à discord.py (2)
Introduction au modèle linéaire généralisé (GLM) par Python
Introduction à Lightning Pytorch
Premiers pas avec le Web Scraping
Introduction aux baies non paramétriques
Introduction à EV3 / MicroPython
Introduction au langage Python
Introduction à la reconnaissance d'image TensorFlow
Introduction à OpenCV (python) - (2)
Introduction à PyQt4 Partie 1
Introduction à l'injection de dépendances
Introduction à Private Chainer
Introduction à l'apprentissage automatique
[Introduction à la simulation] J'ai essayé de jouer en simulant une infection corona ♬
[Introduction à Pandas] J'ai essayé d'augmenter les données d'échange par interpolation de données ♬
AOJ Introduction à la programmation Sujet 1, Sujet 2, Sujet 3, Sujet 4
Introduction au module de papier électronique
Introduction à l'algorithme de recherche de dictionnaire
Introduction à la méthode Monte Carlo
[Mémorandum d'apprentissage] Introduction à vim
Introduction à PyTorch (1) Différenciation automatique
opencv-python Introduction au traitement d'image
Introduction à l'écriture de Cython [Notes]
Introduction à Private TensorFlow
Une introduction à l'apprentissage automatique
[Introduction à cx_Oracle] Présentation de cx_Oracle
Une super introduction à Linux
AOJ Introduction à la programmation Sujet n ° 7, Sujet n ° 8
Introduction à la détection des anomalies 1 principes de base
Introduction à RDB avec sqlalchemy Ⅰ
[Introduction au système] Retracement de Fibonacci ♬
Introduction à l'optimisation non linéaire (I)
Introduction à la communication série [Python]
AOJ Introduction à la programmation Sujet n ° 5, Sujet n ° 6
Introduction au Deep Learning ~ Règles d'apprentissage ~
[Introduction à Python] <liste> [modifier le 22/02/2020]
Introduction à Python (version Python APG4b)
Une introduction à la programmation Python
[Introduction à cx_Oracle] (8e) version de cx_Oracle 8.0
Introduction à discord.py (3) Utilisation de la voix
Introduction à l'optimisation bayésienne
Apprentissage par renforcement profond 1 Introduction au renforcement de l'apprentissage
Super introduction à l'apprentissage automatique
Introduction à Ansible Part «Inventaire»
Série: Introduction à cx_Oracle Contents
[Introduction] Comment utiliser open3d
Introduction à Python pour, pendant