12 Régression de Poisson

Exercice 1 (Questions de cours) C, A, B, A, B, B, C, A

Exercice 2  

• Le modèle linéaire gaussien s’écrit $$\begin{array}{r}{Y}_i=X_{i.}^{\prime }\beta +{\epsilon }_i,{\epsilon }_i\sim \mathcal{N}(0,{\sigma }^2)\end{array}$$ c’est-à-dire $$\begin{array}{r}{Y}_i\sim \mathcal{N}(X_{i.}^{\prime }\beta ,{\sigma }^2)\end{array}$$ La densité en $y$ est $$\begin{array}{rl}f(y)& =\frac{1}{\sqrt{2\pi }}\exp (\frac{-(y-x^{\prime }\beta {)}^2}{2{\sigma }^2})=\exp (\frac{(y(x^{\prime }\beta )-(x^{\prime }\beta {)}^2/2)}{{\sigma }^2}-\frac{1}{2}(\frac{y^2}{{\sigma }^2}+\log (2\pi {\sigma }^2)))\end{array}$$ On a donc $$\begin{array}{rl}\alpha & =x^{\prime }\beta \\ b(\alpha )& ={\alpha }^2/2\\ \varphi & ={\sigma }^2\\ a(\varphi )& ={\sigma }^2\\ c(y,\varphi )& =-\frac{1}{2}(\frac{y^2}{\varphi }+\log (2\pi \varphi )).\end{array}$$ La fonction de lien canonique est par définition $$\begin{array}{r}g(u)=(b^{\prime }{)}^{-1}(u)b^{\prime }(\alpha )=({\alpha }^2/2{)}^{\prime }=\alpha \end{array}$$ L’inverse de $b^{\prime }(.):\alpha \mapsto u=b^{\prime }(\alpha )=\alpha$ est donc $b^{\prime }(.{)}^{-1}:u\mapsto \alpha =u$, c’est le lien identité.

• Le modèle logistique avec répétition en ${\tilde{x}}_t$ (voir paragraphe 11.4.1, p. 270) s’écrit $$\begin{array}{r}{\tilde{Y}}_t={\tilde{y}}_t|X={\tilde{x}}_t\sim \mathcal{B}(n_t,p_{\beta }(x_t))\end{array}$$ où $p_{\beta }(x_t)=\text{logit}(x_t\beta ).$ La probabilité que ${\tilde{Y}}_t$ vaille ${\tilde{y}}_t$ ou densité par rapport à la mesure de comptage au point ${\tilde{y}}_t$ en s’écrit $$\begin{array}{r}f({\tilde{y}}_t)=P({\tilde{Y}}_t={\tilde{y}}_t)=(\genfrac{}{}{0ex}{}{n_t}{{\tilde{y}}_t})p_{\beta }(x_t{)}^{{\tilde{y}}_t}(1-p_{\beta }(x_t){)}^{{\tilde{n}}_t-{\tilde{y}}_t}\end{array}$$ Si il n’y a pas de répétitions en ${\tilde{x}}_t$ on a $n_t=1$ et ${\tilde{y}}_t=0$ ou $1$.

  En modèle GLM on s’intéresse non pas au comptage ${\tilde{Y}}_t$ mais à la moyenne ${\overline{Y}}_t={\tilde{Y}}_t/n_t$, nous avons donc $$\begin{array}{r}f({\overline{y}}_t)=P({\overline{Y}}_t={\overline{y}}_t)=(\genfrac{}{}{0ex}{}{n_t}{n_t{\overline{y}}_t})p_{\beta }(x_t{)}^{n_t{\overline{y}}_t}(1-p_{\beta }(x_t){)}^{n_t-n_t{\overline{y}}_t}\end{array}$$ En introduisant l’exponentielle: $$\begin{array}{rl}f({\overline{y}}_t)& =\exp (\log (\genfrac{}{}{0ex}{}{n_t}{n_t{\overline{y}}_t}))\exp (\log (p_{\beta }(x_t{)}^{n_t{\overline{y}}_t}))\exp (\log ((1-p_{\beta }(x_t){)}^{n_t-n_t{\overline{y}}_t}))\\ & =\exp (n_t{\overline{y}}_t\log p_{\beta }(x_t)+(n_t-n_t{\overline{y}}_t)\log (1-p_{\beta }(x_t))+\log (\genfrac{}{}{0ex}{}{n_t}{n_t{\overline{y}}_t}))\\ & =\exp (n_t{\overline{y}}_t(\log p_{\beta }(x_t)-\log (1-p_{\beta }(x_t)))-n_t\log (1-p_{\beta }(x_t))+\log (\genfrac{}{}{0ex}{}{n_t}{n_t{\overline{y}}_t}))\\ & =\exp (n_t({\overline{y}}_t\log \frac{p_{\beta }(x_t)}{1-p_{\beta }(x_t)}-\log (1-p_{\beta }(x_t)))+\log (\genfrac{}{}{0ex}{}{n_t}{n_t{\overline{y}}_t}))\end{array}$$ Ici nous avons $\alpha =$ ou $$\begin{array}{rl}\alpha & =\log \frac{p_{\beta }(x_t)}{1-p_{\beta }(x_t)}\\ b(\alpha )& =\log (1+\exp \alpha )=\log (1-p_{\beta }(x_t))\\ \varphi & =1/n_t\\ a(\varphi )& =\varphi \\ c({\overline{y}}_t,\varphi )& =\log (\genfrac{}{}{0ex}{}{1/\varphi }{{\overline{y}}_t/\varphi }).\end{array}$$ Pour calculer la fonction de lien canonique nous devons dériver $b(.)$ et inverser la fonction $$\begin{array}{r}{b}^{\prime }(\alpha )=\log (1+\exp \alpha {)}^{\prime }=\frac{\exp (\alpha )}{1+\exp \alpha }=u\\ (b^{\prime }{)}^{-1}(u)=\log \frac{u}{1-u}=\alpha \end{array}$$ c’est le lien logistique.

• Le modèle de poisson s’écrit $$\begin{array}{r}{Y}_i=y_i|X=x_i\sim \mathcal{P}({\lambda }_{\beta }(x_i))\end{array}$$ La probabilité que $Y_i$ vaille $y_i ou densité par rapport à la mesure de comptage au point $y_i$ en s’écrit $$\begin{array}{rl}f(y_i)& =P({\tilde{Y}}_i={\tilde{y}}_i)=\exp (-{\lambda }_{\beta }(x_i))\frac{{\lambda }_{\beta }(x_i{)}^{y_i}}{y_i!}\\ & =\exp (-{\lambda }_{\beta }(x_i)+\log ({\lambda }_{\beta }(x_i{)}^{y_i})-\log (y_i!))\\ & =\exp (y_i\log ({\lambda }_{\beta }(x_i))-{\lambda }_{\beta }(x_i)-\log (y_i!))\end{array}$$ Ici nous avons $$\begin{array}{rl}\alpha & =\log ({\lambda }_{\beta }(x_i))\\ b(\alpha )& =\exp (\alpha )={\lambda }_{\beta }(x_i)\\ \varphi & =1\\ a(\varphi )& =1\\ c({\overline{y}}_t,\varphi )& =-\log (y_i!).\end{array}$$ Pour calculer la fonction de lien canonique nous devons dériver $b(.)$ et inverser la fonction $$\begin{array}{r}{b}^{\prime }(\alpha )=\exp (\alpha )=u\\ (b^{\prime }{)}^{-1}(u)=\log u=\alpha \end{array}$$ c’est le lien log.

Exercice 3  

1. Les moments factoriels d’ordre $r$ d’une variable aléatoire suivant une loi de Poisson valent: $$\begin{array}{rl}\mathbf{E}[X(X-1)\dots (X-r+1)]& =\sum _{k=0}^{\mathrm{\infty }}k(k-1)\dots (k-r+1)\frac{{\lambda }^k}{k!}{e}^{-\lambda }\\ & =e^{-\lambda }\sum _{k=r}^{\mathrm{\infty }}\frac{k(k-1)\dots (k-r+1)}{k!}{\lambda }^k\\ & =e^{-\lambda }{\lambda }^r\sum _{k=r}^{\mathrm{\infty }}\frac{1}{(k-r)!}{\lambda }^{k-r}\\ & ={\lambda }^r.\end{array}$$
   2. Pour $r=1$ nous avons donc $$\begin{array}{r}\mathbf{E}(X)=\lambda \end{array}$$ Pour $r=2$ nous avons donc $$\begin{array}{rl}\mathbf{E}(X(X-1))& =\mathbf{E}(X^2)-\mathbf{E}(X)={\lambda }^2\\ \mathbf{E}(X^2)& ={\lambda }^2+\lambda \end{array}$$ Et nous avons la variance $$\begin{array}{rl}\mathrm{V}(X)& =\mathbf{E}(X^2)-\mathbf{E}(X{)}^2={\lambda }^2+\lambda -{\lambda }^2=\lambda .\end{array}$$

Exercice 4 (Stabilisation de la variance) Tirons notre échantillon de taille $n={10}^6$ selon des loi de poisson de paramètre $\lambda \in \{1,2,\dots ,20\}$

n <- 1e7
lambdas <- 1:20
varX <- rep(0, length(lambdas))
varZ <- rep(0, length(lambdas))
for (l in 1:length(lambdas)){
    ech <- rpois(n, lambdas[l])
    varX[l] <- var(ech)
    ech <- sqrt(ech)
    varZ[l] <- var(ech)
}
print(varX)

 [1]  0.9997369  1.9997818  2.9984626  4.0031435  4.9959479  5.9992567
 [7]  7.0043691  7.9954055  8.9991719  9.9990725 10.9885709 12.0038720
[13] 12.9821786 13.9998158 14.9955777 16.0014936 16.9962143 17.9969313
[19] 19.0054412 19.9929461

print(varZ)

 [1] 0.4022116 0.3897837 0.3396863 0.3059119 0.2856383 0.2753928 0.2693276
 [8] 0.2652610 0.2629826 0.2613118 0.2598299 0.2590222 0.2578215 0.2575848
[15] 0.2568114 0.2565663 0.2560164 0.2556442 0.2554479 0.2549726

Exercice 5 (Stabilisation de la variance) Comme $X_i\stackrel{\mathrm{iid}}{\sim }\mathcal{P}(\lambda )$ nous avons par le TCL que la moyenne d’un $n$ échantillon ${\overline{X}}_n$ converge vers la loi normale et plus exactement $\sqrt{n}({\overline{X}}_n-\lambda )\stackrel{\mathcal{D}}{\to }\mathcal{N}(0,\lambda )$. Posons la fonction $g(.)$ de classe $C^1$ et utilisons un développement de Taylor autour de $\lambda$: $$\begin{array}{rl}g(x)& =g(\lambda )+(x-\lambda )g^{\prime }(\lambda )+o(x-\lambda )\end{array}$$ Ici nous nous intéressons à stabiliser ${\overline{X}}_n$ et nous posons la transformation $g({\overline{X}}_n)$ ce qui donne $$\begin{array}{rl}g({\overline{X}}_n)& =g(\lambda )+({\overline{X}}_n-\lambda )g^{\prime }(\lambda )+o({\overline{X}}_n-\lambda )\\ \sqrt{n}(g({\overline{X}}_n)-g(\lambda ))& =\sqrt{n}({\overline{X}}_n-\lambda )g^{\prime }(\lambda )+\sqrt{n}o({\overline{X}}_n-\lambda )\end{array}$$ Comme $\sqrt{n}({\overline{X}}_n-\lambda )$ converge en distribution donc $({\overline{X}}_n-\lambda )=O_P(1/\sqrt{n})$ donc $\sqrt{n}o({\overline{X}}_n-\lambda )=o_P(1)$ et donc $$\begin{array}{rl}\sqrt{n}(g({\overline{X}}_n)-g(\lambda ))& =\sqrt{n}({\overline{X}}_n-\lambda )g^{\prime }(\lambda )+o_P(1)\end{array}$$ Si $g^{\prime }(\lambda )\ne 0$ on a $$\begin{array}{r}\sqrt{n}(g({\overline{X}}_n)-g(\lambda ))\stackrel{\mathcal{D}}{\to }\mathcal{N}(0,\lambda g^{\prime }(\lambda {)}^2)\end{array}$$ Pour avoir la variance unité, on choisit $\lambda g^{\prime }(\lambda {)}^2)=1$ c’est à dire $$\begin{array}{rl}{g}^{\prime }(\lambda )& =\frac{1}{\sqrt{\lambda }}\\ g(\lambda )& =2\sqrt{\lambda }\end{array}$$ On peut donc prendre la transformée $g(X)=2\sqrt{X}$ (pour avoir la variance unité) ou $g(X)=\sqrt{X}$ et on aura une variance de $1/4$.

Exercice 6  

1. Graphique

   Malaria <- read.csv("../donnees/poissonData3.csv")
   tab <- lapply( split( Malaria$N.malaria, Malaria$Sexe ), table )
   Tab <- matrix( 0,2, max( Malaria$N.malaria )+1 )
   colnames(Tab) <- 0:max( Malaria$N.malaria )
   Tab[1, names(tab[[1]]) ] <- tab[[1]]
   Tab[2, names(tab[[2]]) ] <- tab[[2]]
   barplot( Tab, offset = -Tab[1,])

   

   Les barres sont superposées : les filles puis les garçons. Comme on soustrait à la hauteur totale les garçons (argument offset), on a les effectifs des filles en dessous et ceux des garçons au dessus.

2. Les moyennes par groupe

   aggregate(Malaria$N.malaria, list(Malaria$Sexe), mean)

     Group.1        x
   1       F 4.579012
   2       M 4.794370

   et la différence des logarithme népérien

   diff(log(aggregate(Malaria$N.malaria, list(Malaria$Sexe), mean)[,"x"]))

   [1] 0.04595891

3. Un autre code

   mm <- sapply( split( Malaria$N.malaria, list(Malaria$Sexe) ),mean,na.rm=T)
   round( c( log( mm[1] ), diff( log( mm ) ) ), 5)

         F       M 
   1.52148 0.04596 

4. Régression de Poisson

   mod <- glm(N.malaria ~ 1 + Sexe, data=Malaria, family = poisson)
   mod

   
   Call:  glm(formula = N.malaria ~ 1 + Sexe, family = poisson, data = Malaria)
   
   Coefficients:
   (Intercept)        SexeM  
       1.52148      0.04596  
   
   Degrees of Freedom: 1626 Total (i.e. Null);  1625 Residual
   Null Deviance:      5710 
   Residual Deviance: 5706     AIC: 10510

   Nous retrouvons que le coefficient constant (Intercept) est le logarithme népérien de la moyenne du nombre de visites chez les filles. La modalité fille est la première modalité de la variable Sexe par ordre alphabétique et constitue la modalité de référence. Le coefficient constant est ici le logarithme (qui est la fonction de lien) de la moyenne du nombre de visites chez les filles. L’effet Age est ici la différence des logarithmes ce que nous retrouvons dans le second coefficient.

Exercice 7 (Table de contingence et loi de Poisson)  

1. Importons les données

   Malaria <- read.table("../donnees/poissonData3.csv", sep=",", header=TRUE, stringsAsFactors=TRUE)

   et le tableau de contingence est

   table(Malaria$Sexe, Malaria$Prevention)

   
       Autre Moustiquaire Rien Serpentin/Spray
     F     2          557  223              28
     M     6          543  233              35

2. Lien entre Sexe et Prévention

   chisq.test(Malaria$Sexe,Malaria$Prevention)

   
       Pearson's Chi-squared test
   
   data:  Malaria$Sexe and Malaria$Prevention
   X-squared = 3.1452, df = 3, p-value = 0.3698

   La probabilité critique est de 0.3698 donc nous conservons ${\mathrm{H}}_0$ il y a indépendance entre Sexe et Prévention.

3. Constitution du data-frame

   Y <- as.vector(table(Malaria$Sexe, Malaria$Prevention))
   eff <- table(Malaria$Sexe, Malaria$Prevention)
   Sexe <- factor(rep(c("F", "M"), 4))
   Prevention <- factor(rep(levels(Malaria$Prevention), each=2))
   don <- data.frame(Y, Sexe, Prevention)

4. Modèle de Poisson avec interaction

   mod1 <- glm( Y ~ -1 + Sexe:Prevention, data=don, family=poisson)

   • Il n’y a pas de contraintes ici. On a $2\times 4$ niveaux d’interactions et 8 coefficients estimés. Il s’agit d’une situation analogue à une analyse de variance à 2 facteurs (à 2 et 4 niveaux) avec interaction et sans répétitions. Dans chaque niveau d’interaction $i\times j$ on observe un comptage $Y_{ij}$ (l’effectif) qui est modélisé par un modèle de Poisson. Les comptages suivent une loi de Poisson de paramètre ${\lambda }_{ij}$ dont le log est fonction du niveau de l’interaction $$\begin{array}{r}\log {\lambda }_{ij}={\gamma }_{ij}\end{array}$$ En ANOVA c’est un modèle normal et le lien est l’identité.

   • Le modèle est saturé, il y a une $n=8$ observations et 8 points distincts dans le design: 8 niveaux d’interactions. Par ailleurs on a bien $p=8$ paramètres.

   • Comme le modèle est saturé nous avons que les items ci-dessous représentent la même chose:

     • les ajustements du modèle mod1
     • les observations en chaque point du design (chaque niveau d’interaction)
     • la moyenne en chaque point du design car il n’y a qu’un point par niveau d’interaction

     all(abs(fitted(mod1) - don$Y)<1e-10)

     [1] TRUE

     moyparcellule <- aggregate(don$Y, list(don$Sexe,don$Prevention), mean)[,"x"]
     all(abs(fitted(mod1) - moyparcellule)<1e-10)

     [1] TRUE

     En passant au log, les ajustements (ou les comptages, ou la moyenne par niveau) valent les coefficients ${\beta }_{ij}$:

     all(abs(log(don$Y) - coef(mod1))<1e-10)

     [1] TRUE

5. Modèle de Poisson sans interaction

   mod2 <- glm( Y ~ 1 + Sexe + Prevention, data=don, family=poisson)

   Le comptage $Y_{ij}$ (l’effectif) est modélisé par un modèle de Poisson : les comptages suivent une loi de Poisson de paramètre ${\lambda }_{ij}$ dont le log est modélisé par $$\begin{array}{r}\log {\lambda }_{ij}=\mu +{\alpha }_i+{\beta }_j\end{array}$$

   • Il y a des contraintes identifiantes qui sont celles choisies par défaut: ${\beta }_1={\alpha }_1=0$

         mod2

   
   Call:  glm(formula = Y ~ 1 + Sexe + Prevention, family = poisson, data = don)
   
   Coefficients:
                 (Intercept)                      SexeM  
                    1.381983                   0.008605  
      PreventionMoustiquaire             PreventionRien  
                    4.923624                   4.043051  
   PreventionSerpentin/Spray  
                    2.063693  
   
   Degrees of Freedom: 7 Total (i.e. Null);  3 Residual
   Null Deviance:      1998 
   Residual Deviance: 3.24     AIC: 60.92

   • Le modèle n’est pas saturé, il y a toujours 8 points distincts dans le design et $p=5$ paramètres.
   • Les estimations du modèle tendent à modéliser les effectifs par un effet Sexe et Prevention sans interaction.

6. Comparaison via AIC

   AIC(mod1)

   [1] 63.67658

   AIC(mod2)

   [1] 60.91634

   Le modèle avec le plus petit AIC est retenu: le modèle sans interaction. Ce choix est en accord avec la décision prise en question 2.

Exercice 8 (Table de contingence et probabilité)  

1. La probabilité d’être une femme $$\begin{array}{rl}{\pi }_{1.}& =Pr(F)=Pr(F\cup \mathrm{Autre})+Pr(F\cup \mathrm{Moust}+Pr(F\cup \mathrm{Rien}+Pr(F\cup \mathrm{Ser}/\mathrm{Sp}\\ & =\sum _{j=1}^4{\pi }_{1j}\end{array}$$
2. La probabilité d’utiliser au moyen Autre est $$\begin{array}{rl}{\pi }_{.1}& =Pr(\mathrm{Autre})=Pr(F\cup \mathrm{Autre})+Pr(M\cup \mathrm{Autre})\\ & =\sum _{i=1}^2{\pi }_{i1}\end{array}$$ Toutes les autres lignes ou colonnes sont traitées de manière similaire.
3. Contraintes Il s’agit de probabilité donc leur somme vaut 1: $$\begin{array}{rl}{\pi }_{.1}+{\pi }_{.2}+{\pi }_{.3}+{\pi }_{.4}& =1\\ {\pi }_{1.}+{\pi }_{2.}& =1\end{array}$$ On a aussi d’autres probabilités : $$\begin{array}{rl}\sum i,j{\pi }_{ij}& =1\\ \sum i{\pi }_{i|j}& =1,\mathrm{\forall }j\in \{1,2,3,4\}\\ \sum j{\pi }_{j|i}& =1,\mathrm{\forall }i\in \{1,2\}\end{array}$$
4. Indépendance $$\begin{array}{rl}{\pi }_{ij}& =Pr(\mathrm{Sexe}=i\cup \mathrm{Pr}\mathrm{é}\mathrm{vention}=j)=Pr(\mathrm{Sexe}=i)Pr(\mathrm{Pr}\mathrm{é}\mathrm{vention}=j){\pi }_{i.}{\pi }_{.j}\end{array}$$
5. On part de l’écart proposé et on multiplie par $N$: $$\begin{array}{rl}{\pi }_{ij}& ={\pi }_{i.}{\pi }_{.j}\frac{{\pi }_{ij}}{{\pi }_{i.}{\pi }_{.j}}\\ N{\pi }_{ij}& =N{\pi }_{i.}N{\pi }_{.j}\frac{N{\pi }_{ij}}{N{\pi }_{i.}N{\pi }_{.j}}\log (N{\pi }_{ij})& =\log N{\pi }_{i.}+\log N{\pi }_{.j}+\log \frac{N{\pi }_{ij}}{N{\pi }_{i.}N{\pi }_{.j}}\\ \log ({\lambda }_{ij})& ={\alpha }_i+{\beta }_j+{\gamma }_{ij}.\end{array}$$ On a comme contrainte $$\begin{array}{rlr}{\pi }_{1.}+{\pi }_{2.}& =1=\frac{1}{N}(\exp ({\alpha }_1)+\exp ({\alpha }_2)){\pi }_{.1}+{\pi }_{.2}+{\pi }_{.3}+{\pi }_{.4}& =1=\frac{1}{N}(\sum _{j=1}^4\exp ({\alpha }_1{\pi }_{.j})\end{array}$$ On peut aussi exprimer les contraintes sur les $\{{\gamma }_{ij}\}$ avec les contraintes sur les probabilités conditionnelles.

Exercice 9 (Loi Multinomiale) Soit $N$ qui suit une multinomiale : $$\begin{array}{rl}Pr(N_1=n_1,\dots ,N_K=n_K)& =\frac{n!}{n_1!\dots n_K!}\prod _{k=1}^K{\pi }_k^{n_k}.\end{array}$$ en introduisant l’exponentielle $$\begin{array}{rl}f(n1,\dots ,n_K)& =\exp (\log (\frac{n!}{n_1!\dots n_K!}))\exp (\log (\prod _{k=1}^K{\pi }_k^{n_k}))\\ f(n1,\dots ,n_K)& =\exp (\sum _{k=1}^K{n}_k\log {\pi }_k+\log (\frac{n!}{n_1!\dots n_K!}))\end{array}$$

Exercice 10  

1. Nous souhaitons modéliser des effectifs répartis dans des cases $i\times j$ ce qui représente le champs d’action de la loi multinomiale.

   • Un premier modèle possible est de modéliser toute la table d’un coup. On désigne par $N_{ij},i\in \{1,2\},j\in \{1,2,3\}$ l’effectif (aléatoire) du croisement $(i,j)$ et par ${\pi }_{ij}$ la probabilité (inconnue) d’appartenir à ce croisement (voir le tableau de l’énoncé). Ainsi $N_{1,3}$ correspondra par exemple au croisement (,). Nous avons donc l’effectif total de la table $n=631$ et nous utilisons toute la table d’un coup et nous modélisons les effectifs $N_{ij}$ par une loi multinomiale $\mathcal{M}(n,{\pi }_{11},{\pi }_{12},\dots ,{\pi }_{23})$. Pour une case $i,j$ nous avons pour les probabilités: $$\begin{array}{rl}Pr(\mathtt{Gr}=i,\mathtt{In}=j)& ={\pi }_{ij}.\end{array}$$ Si nous soupçonnons qu’il n’y ait pas de lien entre Groupe et Intention nous pouvons simplifier le modèle via l’indépendance et avoir pour une case $i,j$ $$\begin{array}{rl}Pr(\mathtt{Gr}=i,\mathtt{In}=j)& =Pr(\mathtt{Gr}=i)Pr(\mathtt{In}=j)={\pi }_{i.}{\pi }_{.j}.\end{array}$$ Cette approche est celle du test du ${\chi }^2$ et nous ne la développerons pas dans cet exercice.
   • Une seconde approche (utilisée ici) est de modéliser chaque ligne du tableau.
     • Si nous soupçonnons qu’il n’y ait pas de lien entre Groupe et Intention nous pouvons dire que chaque ligne est modélisée par une loi multinomiale $\mathcal{M}(n_i,{\pi }_{.1},{\pi }_{.2},{\pi }_{.3})$. D’une ligne à l’autre la loi est presque la même (l’effectif $n_i$ change): la répartition dans les cases suit les mêmes probabilité puisque l’on soupçonne qu’il n’y a pas de lien entre Groupe et Intention.
     • Si nous soupçonnons qu’il y ait un lien entre Groupe et Intention nous pouvons dire que chaque ligne est modélisée par une loi multinomiale $\mathcal{M}(n_i,{\pi }_{i1},{\pi }_{i2},{\pi }_{i3})$. D’une ligne à l’autre la loi est différente: l’effectif $n_i$ change et la répartition dans les cases (${\pi }_{i1},{\pi }_{i2},{\pi }_{i3}$) change aussi ; c’est l’objet de la question 2.

   Pour la ligne $i$ on a comme observation les effectifs (d’une multinomiale) $n_{i1},n_{i2},n_{i3}$. Leur somme vaut $n_{i.}$. On doit répartir dans 3 groupes mais les probabilités de répartitions sont les mêmes d’une ligne à l’autre elles sont notées ${\pi }_{.j}$. On a donc la probabilité d’observer $n_{i1},n_{i2},n_{i3}$: $$\begin{array}{rl}Pr(N_{i1}=n_{i1},N_{i2}=n_{i2},N_{i3}=n_{i3})& =\frac{n_{i.}!}{n_{i1}!n_{i1}!n_{i3}!}\prod _j({\pi }_{.j}{)}^{n_{ij}}.\end{array}$$ Pour toutes les données nous avons alors $$\begin{array}{rl}Pr(N_{11}=n_{11},N_{12}=n_{12},,N_{13}=n_{13})Pr(N_{21}=n_{11},N_{22}=n_{22},N_{23}=n_{23})& ==\prod _{i=1}^2\frac{n_{i.}!}{n_{i1}!n_{i1}!n_{i3}!}\prod _j({\pi }_{.j}{)}^{n_{ij}}.\end{array}$$

2. Le modèle avec lien entre Groupe et Intention permet d’écrire la probabilité de l’échantillon comme: $$\begin{array}{rl}Pr(N_{11}=n_{11},N_{12}=n_{12},,N_{13}=n_{13})Pr(N_{21}=n_{11},N_{22}=n_{22},N_{23}=n_{23})& ==\prod _{i=1}^2\frac{n_{i.}!}{n_{i1}!n_{i1}!n_{i3}!}\prod _j({\pi }_{ij}{)}^{n_{ij}}.\end{array}$$

3. Calculons les log-vraisemblances pour les deux modèles.

   • Modèle sans effet Groupe: $$\begin{array}{rl}\mathcal{L}(Y,{\pi }_{ij})=& \sum _{i,j}{n}_{ij}\log {\pi }_{.j}+\mathrm{cte}.\end{array}$$ sous les contraintes $\sum _j{\pi }_{.j}=1$ et les probabilités sont positives. La positivité peut être imposée en écrivant ${\pi }_{.j}=\exp ({\mu }_j)$ et la somme à 1 impose alors $${\pi }_{.j}=\frac{\exp ({\mu }_j)}{\sum _j\exp ({\mu }_j)}.$$ En utilisant cette paramétrisation nous avons $$\begin{array}{rl}\mathcal{L}(Y,{\mu }_j)=& \sum _{i,j}{n}_{ij}{\mu }_j-\sum _{i,j}{n}_{ij}\log (\sum _j\exp ({\mu }_j))+\mathrm{cte}\\ \mathcal{L}(Y,{\mu }_j)=& \sum _j{n}_{i.}{\mu }_j-\sum _i{n}_{i.}\log (\sum _j\exp ({\mu }_j))+\mathrm{cte}\end{array}$$

   • Modèle avec effet Groupe: $$\begin{array}{rl}\mathcal{L}(Y,\pi )=& \sum _{i,j}{n}_{ij}\log {\pi }_{ij}+\sum _{i,j}\log (n!)-\sum _{i,j}\log (n_{ij}!)=\sum _{i,j}{n}_{ij}\log {\pi }_{ij}+\mathrm{cte}.,\end{array}$$ sous les contraintes $\sum _j{\pi }_{ij}=1,\mathrm{\forall }i$ et les probabilités sont positives. La positivité peut être imposée en écrivant ${\pi }_{ij}=\exp ({\eta }_{ij})$ et la somme à 1 impose alors $${\pi }_{ij}=\frac{\exp ({\eta }_{ij})}{\sum _j\exp ({\eta }_{ij})}.$$ En utilisant cette paramétrisation nous avons $$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathcal{L}(Y,{\eta }_{ij})=& \sum _{i,j}{n}_{ij}{\eta }_{ij}-\sum _{i,j}{n}_{ij}\log (\sum _j\exp ({\eta }_{ij}))+\mathrm{cte}\\ \mathcal{L}(Y,{\eta }_{ij})=& \sum _j{n}_{ij}{\eta }_{ij}-\sum _i{n}_{i.}\log (\sum _j\exp ({\eta }_{ij}))+\mathrm{cte}\end{array}\end{array}$$

4. Le modèle de poisson complet est $$\begin{array}{rl}Y& \sim \mathcal{P}({\lambda }_{ij})\\ \log ({\lambda }_{ij})& =\mu +{\alpha }_i+{\beta }_j+{\gamma }_{ij}\end{array}$$ Ce modèle est sur-paramétré et nécessite des contraintes pour être estimable. On peut choisir $$\begin{array}{rl}{\alpha }_1& =0,\\ {\beta }_1& =0,\\ {\gamma }_{11}={\gamma }_{12}=\dots =\gamma 31=0\end{array}$$ Il donne les mêmes prévisions que le modèle plus simple suivant: $$\begin{array}{rl}Y& \sim \mathcal{P}({\lambda }_{ij})\\ \log ({\lambda }_{ij})& =\mu +{\alpha }_i+{\gamma }_{ij}\end{array}$$ Ce modèle est sur-paramétré et nécessite des contraintes pour être estimable. On peut choisir $$% $$\begin{array}{rl}{\alpha }_1& =0,\\ {\gamma }_{11}={\gamma }_{12}=\dots =\gamma 31=0{\gamma }_{11}={\gamma }_{21}=0\end{array}$$

   Ils donnent les mêmes prévisions que le modèle sans contraintes suivant: $$\begin{array}{rl}Y& \sim \mathcal{P}({\lambda }_{ij})\\ \log ({\lambda }_{ij})& =m_{ij}\end{array}$$ Pour ce dernier modèle, nous avons l’écriture suivante (qui correspond bien à l’écriture de l’équation ci-dessus) et où l’on voit apparaître sous forme matricielle le modèle saturé (la matrice de design est $6\times 6$ et de plein rang): $$\begin{array}{rl}\left[\begin{array}{c}\log ({\lambda }_{11})\\ \log ({\lambda }_{21})\\ \log ({\lambda }_{12})\\ \log ({\lambda }_{22})\\ \log ({\lambda }_{13})\\ \log ({\lambda }_{23})\end{array}\right]& =\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{m}_{11}\\ m_{21}\\ m_{12}\\ m_{22}\\ m_{13}\\ m_{23}\end{array}\right]\end{array}$$ Pour le second modèle nous avons $$\begin{array}{rl}\left[\begin{array}{c}\log ({\lambda }_{11})\\ \log ({\lambda }_{21})\\ \log ({\lambda }_{12})\\ \log ({\lambda }_{22})\\ \log ({\lambda }_{13})\\ \log ({\lambda }_{23})\end{array}\right]& =\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 1& 1& 0& 0& 0& 0\\ 1& 0& 1& 0& 0& 0\\ 1& 1& 0& 1& 0& 0\\ 1& 0& 0& 0& 1& 0\\ 1& 1& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}\mu \\ {\alpha }_2\\ {\gamma }_{12}\\ {\gamma }_{22}\\ {\gamma }_{13}\\ {\gamma }_{23}\end{array}\right]=\left[\begin{array}{c}\mu +{\alpha }_1+{\gamma }_{11}\\ \mu +{\alpha }_1+{\gamma }_{21}\\ \mu +{\alpha }_2+{\gamma }_{12}\\ \mu +{\alpha }_2+{\gamma }_{22}\\ \mu +{\alpha }_2+{\gamma }_{13}\\ \mu +{\alpha }_3+{\gamma }_{23}\end{array}\right]\end{array}$$ Pour le premier modèle nous avons $$\begin{array}{rl}\left[\begin{array}{c}\log ({\lambda }_{11})\\ \log ({\lambda }_{21})\\ \log ({\lambda }_{12})\\ \log ({\lambda }_{22})\\ \log ({\lambda }_{13})\\ \log ({\lambda }_{23})\end{array}\right]& =\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 1& 1& 0& 0& 0& 0\\ 1& 0& 1& 0& 0& 0\\ 1& 1& 1& 0& 1& 0\\ 1& 0& 0& 1& 0& 0\\ 1& 1& 0& 1& 0& 1\end{array}\right]\left[\begin{array}{c}\mu \\ {\alpha }_2\\ {\beta }_2\\ {\beta }_3\\ {\gamma }_{13}\\ {\gamma }_{23}\end{array}\right]=\left[\begin{array}{c}\mu +{\alpha }_1+{\beta }_1+{\gamma }_{11}\\ \mu +{\alpha }_1+{\beta }_2+{\gamma }_{21}\\ \mu +{\alpha }_2+{\beta }_3+{\gamma }_{12}\\ \mu +{\alpha }_2+{\beta }_1+{\gamma }_{22}\\ \mu +{\alpha }_2+{\beta }_2+{\gamma }_{13}\\ \mu +{\alpha }_3+{\beta }_3+{\gamma }_{23}\end{array}\right]\end{array}$$

5. Calculons la vraisemblance du modèle avec contrainte pour les observations $\{n_{ij}\}$: $$\begin{array}{r}\prod _{i,j}{e}^{-{\lambda }_{ij}}\frac{{\lambda }_{ij}^{n_{ij}}}{n_{ij}}\end{array}$$ dont le logarithme est $$\begin{array}{rl}\mathcal{L}(n,\lambda )& =-\sum _{i,j}{\lambda }_{ij}+\sum _{i,j}{n}_{ij}\log {\lambda }_{ij}+cte.\end{array}$$ Comme $\log ({\lambda }_{ij})=\mu +{\alpha }_i+{\gamma }_{ij}$, on en déduit $$\begin{array}{rl}\mathcal{L}(n,\mu ,{\alpha }_i,{\gamma }_{ij})& =-\sum _i\sum _j\exp (\mu +{\alpha }_i+{\gamma }_{ij})+\sum _{i,j}{n}_{ij}(\mu +{\alpha }_i+{\gamma }_{ij})+cte.\end{array}$$ Posons ${\tau }_i=\sum _j\exp ((\mu +{\alpha }_i)+{\gamma }_{ij})$ et on a donc $$\begin{array}{rl}\mathcal{L}(n,\mu ,{\tau }_i,{\gamma }_{ij})& =-\sum _i{\tau }_i+\sum _{i,j}{n}_{ij}(\mu +{\alpha }_i)+\sum _{i,j}{n}_{ij}{\gamma }_{ij}+cte.\end{array}$$ Exprimons $(\mu +{\alpha }_i)$ en fonction de ${\tau }_i$: $$\begin{array}{rl}{\tau }_i& =\sum _j\exp ((\mu +{\alpha }_i)+{\gamma }_{ij})=\exp (\mu +{\alpha }_i)\sum _j\exp ({\gamma }_{ij})\\ \exp (\mu +{\alpha }_i)& =\frac{{\tau }_i}{\sum _j\exp ({\gamma }_{ij})}\\ \mu +{\alpha }_i& =\log \frac{{\tau }_i}{\sum _j\exp ({\gamma }_{ij})}=\log {\tau }_i-\log (\sum _j\exp ({\gamma }_{ij}))\end{array}$$ et nous introduisons cette dernière équation dans la log-vraisemblance: $$\begin{array}{rl}\mathcal{L}(n,\mu ,{\tau }_i,{\gamma }_{ij})& =-\sum _i{\tau }_i+\sum _{i,j}{n}_{ij}(\log {\tau }_i-\log (\sum _j\exp ({\gamma }_{ij})))+\sum _{i,j}{n}_{ij}{\gamma }_{ij}+cte\\ \mathcal{L}(n,\mu ,{\tau }_i,{\gamma }_{ij})& =\sum _i(n_{i.}\log {\tau }_i-{\tau }_i)+\sum _{i,j}{n}_{ij}{\gamma }_{ij}-\sum _i{n}_{i.}\log (\sum _j\exp ({\gamma }_{ij}))+cte.\end{array}$$ La vraisemblance est séparée en 2 parties: le premier terme $\sum _i(n_{i.}\log {\tau }_i-{\tau }_i)$ dépend de ${\tau }_i$ alors que les 2 derniers dépendent de ${\gamma }_{ij}$. On reconnaît dans les 2 derniers termes la vraisemblance de la modélisation multinomiale avec effet Groupe (équation 1). La maximisation donnerait les mêmes ${\gamma }_{ij}$ si il n’y avait pas les contraintes… Regardons cela:

   • Du côté de la régression Poisson on a pour chaque cellule $(i,j)$ la modélisation $\mu +{\alpha }_i+{\gamma }_{ij}$ avec comme contrainte par exemple $$\begin{array}{rl}{\alpha }_1& =0,\\ {\beta }_1& =0,\\ {\gamma }_{11}={\gamma }_{12}={\gamma }_{21}={\gamma }_{22}=0\end{array}$$ Pour passer sur l’échelle des probabilités on utilise toujours l’exponentielle et on divise par la somme sur chaque ligne $i$ (les probabilités doivent sommer à 1 pour chaque ligne) $$\begin{array}{rl}{\pi }_{ij}& =\frac{\exp (\mu +{\alpha }_i+{\gamma }_{ij})}{\sum _j\exp (\mu +{\alpha }_i+{\gamma }_{ij})}\end{array}$$
   • Du côté de la régression Multinomiale on a pour chaque cellule $(i,j)$ directement ${\gamma }_{ij}$ et pour passer sur l’échelle des probabilités on utilise toujours l’exponentielle et on divise par la somme sur chaque ligne $i$ (les probabilités doivent sommer à 1 pour chaque ligne) $$\begin{array}{rl}{\pi }_{ij}& =\frac{\exp ({\gamma }_{ij})}{\sum _j\exp ({\gamma }_{ij})}\end{array}$$

   Le même calcul que ci-dessus peut être conduit avec le modèle sans contrainte avec les $m_{ij}$ et on obtient (non demandé) $$\begin{array}{rl}\mathcal{L}(n,m_{ij})& =-\sum _{i,j}\exp (m_{ij})+\sum _{i,j}{n}_{ij}{m}_{ij}+cte\end{array}$$ Pour passer sur l’échelle des probabilités on utilise toujours l’exponentielle et on divise par la somme sur chaque ligne $i$ (les probabilités doivent sommer à 1 pour chaque ligne) $$\begin{array}{rl}{\pi }_{ij}& =\frac{\exp (m_{ij})}{\sum _j\exp (m_{ij})}\end{array}$$