[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
- La variable dépendante $ y_i $ suit la famille de distribution exponentielle.
- Chaque variable dépendante est indépendante.
Dans l'hypothèse
- En supposant la distribution de probabilité de la variable de probabilité $ y_i $. -Définir un modèle de régression linéaire avec la valeur attendue de $ y_i $ $ \ mu $. --Créez une fonction de vraisemblance pour la distribution de probabilité supposée de $ y_i $ et trouvez $ \ beta (\ hat {\ beta}) $ pour la maximiser.
- L'estimation la plus probable $ \ hat {\ beta} $ a une cohérence et une quasi-normalité.
Estimez le modèle comme.
Dans GEE
--Relaxez l'hypothèse que la variable dépendante $ y_i $ suit la famille de distribution exponentielle.
- Redéfinissez l'hypothèse que la variable dépendante est indépendante. --Correspond à une distribution inconnue en utilisant la pseudo-probabilité. --Introduire la matrice de corrélation de travail $ R_i (\ alpha) $ pour représenter la structure de corrélation.
- L'équation d'estimation généralisée est divisée en deux étapes et résolue par calcul numérique.
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)
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()
data.lage.hist()
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)
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)
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)
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() 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)
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)
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}
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
