Introduction à la modélisation statistique pour le modèle linéaire généralisé d'analyse de données (GLM)
Continuant de la fois précédente (Chapitre 2), Chapitre 3. Etudier un modèle statistique qui intègre des ** variables explicatives **.
3.1 Exemple: lorsque le nombre moyen de graines varie d'un individu à l'autre
Comme la différence avec la fois précédente, la taille corporelle $ x_i $, qui est l'un des attributs de l'individu, et la présence ou l'absence de contrôle (ne rien faire: ** C **, ajouter de l'engrais: ** T **) sont utilisées comme variables explicatives.
Lire les données
[Ici] Téléchargez et utilisez data3a.csv dans le chapitre 3 de ().
Les données utilisées cette fois sont appelées type DataFrame.
# d <- read.csv("data3a.csv")
>>> import pandas
>>> d = pandas.read_csv("data3a.csv")
>>> d
y x f
0 6 8.31 C
1 6 9.44 C
2 6 9.50 C
3 12 9.07 C
...(Omis)...
98 7 10.86 T
99 9 9.97 T
[100 rows x 3 columns]
# d$x
>>> d.x
0 8.31
1 9.44
2 9.50
...(Omis)...
98 10.86
99 9.97
Name: x, Length: 100, dtype: float64
# d$y
>>> d.y
0 6
1 6
2 6
...(Omis)...
98 7
99 9
Name: y, Length: 100, dtype: int64
# d$f
>>> d.f = d.f.astype('category')
>>> d.f
0 C
1 C
2 C
3 C
...(Omis)...
98 T
99 T
Name: f, Length: 100, dtype: category
Categories (2, object): [C < T]
Dans R, la colonne f est appelée ** facteur ** (facteur), En python (pandas), c'est une catégorie.
Pour vérifier le type de données
# class(d)Impossible de confirmer pendant un instant.
# class(d$y)Classe Yara(d$x)Classe Yara(d$f)
>>> d.y.dtype
dtype('int64')
>>> d.x.dtype
dtype('float64')
>>> d.f.dtype
dtype('O')
Pour un aperçu de la trame de données, voir Série (identique à la dernière fois)
# summary(d)
>>> d.describe()
y x
count 100.000000 100.000000
mean 7.830000 10.089100
std 2.624881 1.008049
min 2.000000 7.190000
25% 6.000000 9.427500
50% 8.000000 10.155000
75% 10.000000 10.685000
max 15.000000 12.400000
>>> d.f.describe()
count 100
unique 2
top T
freq 50
Name: f, dtype: object
Visualisation
Je n'ai pas trouvé de moyen de coder en couleur le diagramme de dispersion en une seule ligne comme R. .. ..
# plot(d$x, d$y, pch=c(21, 19)[d$f])
>>> plt.plot(d.x[d.f == 'C'], d.y[d.f == 'C'], 'bo')
[<matplotlib.lines.Line2D at 0x1184d0490>]
>>> plt.plot(d.x[d.f == 'T'], d.y[d.f == 'T'], 'ro')
[<matplotlib.lines.Line2D at 0x111b58dd0>]
>>> plt.show()
>>> d.boxplot(column='y', by='f')
<matplotlib.axes._subplots.AxesSubplot at 0x11891a9d0>
Ce que vous pouvez voir sur ces chiffres
--Il semble que le nombre de graines $ y $ augmente à mesure que la taille du corps $ x_i $ augmente.
- L'effet de l'engrais ne semble pas être $ f $
Modèle statistique de régression de Poisson
Créez un modèle statistique qui semble être en mesure d'exprimer les données de nombre de graines, qui sont les données de comptage.
Modèle statistique qui dépend de la taille corporelle $ x_i $ de l'individu $ i $
** Variable explicative **: $ x_i $ ** Variable de réponse **: $ y_i $
La probabilité $ p (y_i | \ lambda_i) $ que le nombre de graines soit $ y_i $ chez un individu $ i $ suit la distribution de Poisson.
p(y_i|\lambda_i) = \frac{\lambda _i ^{y_i} \exp (-\lambda_i)}{y_i !}
Supposer. Ici, le nombre moyen de graines $ \ lambda_i $ d'un individu $ i $ est
\begin{eqnarray}
\lambda_i &=& \exp(\beta_1 + \beta_2 x_i )
\log \lambda_i &=& \beta_1 + \beta_2 x_i
\end{eqnarray}
Supposer que. $ \ Beta_1 $ s'appelle ** section **, $ \ beta_2 $ s'appelle ** tilt **, Le côté droit $ \ beta_1 + \ beta_2 x_i $ est appelé ** prédicteur linéaire **.
Si ($ \ lambda_i $ function) = (prédicteur linéaire), La fonction sur le côté gauche est appelée ** fonction de lien **.
Retour de Poisson
La régression de Poisson trouve $ \ beta_1, \ beta_2 $ qui maximise l'équation suivante.
\log L(\beta_1, \beta_2) = \sum_i \log \frac{\lambda_i^{y_i}\exp(-\lambda_i)}{y_i !}
Dans R, il semble que ce soit un coup si vous utilisez quelque chose comme glm.
Avec python, je dois travailler sans cesse en utilisant la bibliothèque statsmodels ... (Reference)
# fit <- glm(y ~ x, data=d, familiy=poisson)
>>> import statsmodels.api as sm
>>> import statsmodel.formula.api as smf
>>> model = smf.glm('y ~ x', data=d family=sm.families.Poisson())
>>> results = model.fit()
>>> print(results.summary())
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 98
Model Family: Poisson Df Model: 1
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -235.39
Date:Mois, 01 6 2015 Deviance: 84.993
Time: 23:06:45 Pearson chi2: 83.8
No. Iterations: 7
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
const 1.2917 0.364 3.552 0.000 0.579 2.005
x 0.0757 0.036 2.125 0.034 0.006 0.145
==============================================================================
>>> results.llf
-235.38625076986077
Il semble qu'il s'appelle Intercept en R, En Python, le const du bas semble représenter $ \ beta_1 $ (section) et x représente $ \ beta_2 $ (inclinaison).
En d'autres termes, les estimations les plus probables sont $ \ beta_1 = 1,2917 $ et $ \ beta_2 = 0,0757 $.
L'erreur std indique l '** erreur standard **, qui indique la variation de $ \ beta_1 et \ beta_2 $.
Si le même nombre de données différentes est rééchantillonné par la même méthode d'enquête, le montant d'estimation le plus probable changera également et le degré de variation changera.
z est une statistique appelée ** valeur z **, qui est l'estimation la plus probable divisée par l'erreur standard.
Il semble que ce z puisse être utilisé pour estimer l '** intervalle de confiance de Wald **. Apparemment, cet intervalle de confiance n'inclut pas 0 dans la valeur que l'estimation la plus probable peut prendre [voir](http://hosho.ees.hokudai.ac.jp/~kubo/log/2009/0101. html # 08)
Prédiction par modèle de régression de Poisson
Le résultat de l'estimation de la régression de Poisson est utilisé pour prédire le nombre moyen de graines $ \ lambda $ à la taille corporelle $ x $.
\lambda = \exp(1.29 + 0.0757x)
Est illustré en utilisant.
# plot(d$x, d$y, pch=c(21, 19)[d$f])
# xx <- seq(min(d$x), max(d$x), length=100)
# lines(xx, exp(1.29 + 0.0757 * xx), lwd = 2)
>>> plt.plot(d.x[d.f == 'C'], d.y[d.f == 'C'], 'bo')
>>> plt.plot(d.x[d.f == 'T'], d.y[d.f == 'T'], 'ro')
>>> xx = [d.x.min() + i * ((d.x.max() - d.x.min()) / 100.) for i in range(100)]
>>> yy = [numpy.exp(1.29 + 0.0757 * i )for i in xx]
>>> plt.plot(xx, yy)
>>> plt.show()
3.5 Modèle statistique avec variables explicatives de type facteur
Considérons un modèle statistique qui incorpore l'application d'engrais $ f_i $ comme variable explicative.
Dans GLM, ** variable factice ** est utilisée comme variable explicative du type de facteur. Autrement dit, la valeur moyenne du modèle,
\begin{eqnarray}
\lambda_i &=& \exp (\beta_1 + \beta_3 d_i) \\
\therefore d_i &=& \left\{ \begin{array}{ll}
0 & (f_i =Pour C) \\
1 & (f_i =Pour T)
\end{array} \right.
\end{eqnarray}
Sera écrit.
# fit.f <- glm(y ~ f, data=d, family=poisson)
>>> model = smf.glm('y ~ f', data=d, familiy=sm.families.Poisson())
>>> model.fit().summary()
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 98
Model Family: Poisson Df Model: 1
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -237.63
Date:bois, 04 6 2015 Deviance: 89.475
Time: 17:20:11 Pearson chi2: 87.1
No. Iterations: 7
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept 2.0516 0.051 40.463 0.000 1.952 2.151
f[T.T] 0.0128 0.071 0.179 0.858 -0.127 0.153
==============================================================================
>>> model.llf
-237.62725696068682
À partir de ce résultat, si $ f_i $ de l'individu $ i $ est C,
\lambda_i = \exp(2.05 + 0.0128 * 0) = \exp(2.05) = 7.77
Si T
\lambda_i = \exp(2.05 + 0.0128 * 1) = \exp(2.0628) = 7.87
On prévoit que le nombre moyen de graines n'augmentera que légèrement lorsque l'engrais est administré.
On peut dire que la vraisemblance logarithmique (llf) n'est pas applicable car elle est devenue plus petite.
3.6 Modèle statistique avec variables explicatives de type quantitatif + factoriel
Considérons un modèle statistique qui intègre à la fois la taille corporelle individuelle $ x_i $ et le traitement de fertilisation $ f_i $.
\log \lambda_i = \beta_1 + \beta_2 x_i + \beta_3 d_i
>>> model = smf.glm('y ~ x + f', data=d, family=sm.families.Poisson())
>>> result = model.fit()
>>> result.summary()
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 97
Model Family: Poisson Df Model: 2
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -235.29
Date:bois, 04 6 2015 Deviance: 84.808
Time: 17:31:34 Pearson chi2: 83.8
No. Iterations: 7
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept 1.2631 0.370 3.417 0.001 0.539 1.988
f[T.T] -0.0320 0.074 -0.430 0.667 -0.178 0.114
x 0.0801 0.037 2.162 0.031 0.007 0.153
==============================================================================
>>> result.llf
-235.29371924249369
Puisque la vraisemblance logarithmique maximale est la plus grande, on peut dire qu'elle s'adapte mieux que les deux modèles statistiques ci-dessus.
Facilité de compréhension de la fonction de lien logarithmique: effet à multiplier
En transformant le modèle quantitatif ci-dessus + modèle statistique de type facteur,
\begin{eqnarray}
\log \lambda_i &=& \beta_1 + \beta_2 x_i + \beta_3 d_i \\
\lambda_i &=& \exp(\beta_1 + \beta_2 x_i + \beta_3 d_i) \\
\lambda_i &=& \exp(\beta_1) * \exp(\beta_2 x_i) * \exp(\beta_3 d_i) \\
\lambda_i &=& \exp(constant) * \exp(Effet de la taille du corps) * \exp(Effet du traitement de fertilisation)
\end{eqnarray}
On voit que les effets se multiplient.
S'il n'y a pas de fonction de lien, elle est appelée ** fonction de lien constant **.
Dans GLM
Le cas de la distribution normale et de la fonction de lien oral est appelé ** modèle linéaire ** (** modèle linéaire, LM ) ou ** modèle linéaire général ** ( modèle linéaire général **).
Recommended Posts