Contrastes polynomiaux pour la régression

Juste pour récapituler (et au cas où les hyperliens OP échoueraient à l'avenir), nous examinons un ensemble hsb2de données en tant que tel:

   id     female race ses schtyp prog read write math science socst
1  70        0    4   1      1    1   57    52   41      47    57
2 121        1    4   2      1    3   68    59   53      63    61
...
199 118      1    4   2      1    1   55    62   58      58    61
200 137      1    4   3      1    2   63    65   65      53    61

qui peut être importé ici .

Nous transformons la variable readen variable ordonnée / ordinale:

hsb2$readcat<-cut(hsb2$read, 4, ordered = TRUE)
(means = tapply(hsb2$write, hsb2$readcat, mean))
 (28,40]  (40,52]  (52,64]  (64,76] 
42.77273 49.97849 56.56364 61.83333

Maintenant , nous sommes tous mis à courir juste un ANOVA régulier - oui, il est R, et nous avons essentiellement une variable dépendante continue, writeet une variable explicative avec plusieurs niveaux, readcat. Dans R, nous pouvons utiliserlm(write ~ readcat, hsb2)

1. Génération de la matrice de contraste:

Il existe quatre niveaux différents pour la variable ordonnée readcat, nous aurons donc contrastes. $n-1=3$

table(hsb2$readcat)

(28,40] (40,52] (52,64] (64,76] 
     22      93      55      30

Tout d'abord, allons-y pour l'argent, et jetons un coup d'œil à la fonction R intégrée:

contr.poly(4)
             .L   .Q         .C
[1,] -0.6708204  0.5 -0.2236068
[2,] -0.2236068 -0.5  0.6708204
[3,]  0.2236068 -0.5 -0.6708204
[4,]  0.6708204  0.5  0.2236068

Décortiquons maintenant ce qui s'est passé sous le capot:

scores = 1:4  # 1 2 3 4 These are the four levels of the explanatory variable.
y = scores - mean(scores) # scores - 2.5

$y = \small [-1.5, -0.5, 0.5, 1.5]$

$\small \text{seq_len(n) - 1} = [0, 1, 2, 3]$

n = 4; X <- outer(y, seq_len(n) - 1, "^") # n = 4 in this case

$\small\begin{bmatrix} 1&-1.5&2.25&-3.375\\1&-0.5&0.25&-0.125\\1&0.5&0.25&0.125\\1&1.5&2.25&3.375 \end{bmatrix}$

Que s'est-il passé là-bas? le outer(a, b, "^")élève les éléments de aaux éléments de b, de sorte que la première colonne résulte des opérations, , , et ; la deuxième colonne de , , et ; le troisième de $\small(-1.5)^0$ $\small(-0.5)^0$ $\small 0.5^0$ $\small 1.5^0$ $\small(-1.5)^1$ $\small(-0.5)^1$ $\small0.5^1$ $\small1.5^1$ $\small(-1.5)^2=2.25$ , , et ; et le quatrième, , , et . $\small(-0.5)^2 = 0.25$ $\small0.5^2 = 0.25$ $\small1.5^2 = 2.25$ $\small(-1.5)^3=-3.375$ $\small(-0.5)^3=-0.125$ $\small0.5^3=0.125$ $\small1.5^3=3.375$

Ensuite, nous faisons une décomposition orthonormée de cette matrice et prenons la représentation compacte de Q ( ). Certains des fonctionnements internes des fonctions utilisées dans la factorisation QR dans R utilisées dans ce post sont expliqués plus en détail ici . $QR$ c_Q = qr(X)$qr

$\small\begin{bmatrix} -2&0&-2.5&0\\0.5&-2.236&0&-4.584\\0.5&0.447&2&0\\0.5&0.894&-0.9296&-1.342 \end{bmatrix}$

... dont nous enregistrons uniquement la diagonale ( z = c_Q * (row(c_Q) == col(c_Q))). Ce qui se trouve dans la diagonale: juste les entrées "inférieures" de la partie de la décomposition Juste? eh bien, non ... Il s'avère que la diagonale d'une matrice triangulaire supérieure contient les valeurs propres de la matrice! $\bf R$ $QR$

Ensuite, nous appelons la fonction suivante:, raw = qr.qy(qr(X), z)dont le résultat peut être répliqué "manuellement" par deux opérations: 1. Transformer la forme compacte de , c'est -à- dire en , une transformation qui peut être réalisée avec , et 2. Effectuer la multiplication matricielle , comme dans . $Q$ qr(X)$qr $Q$ Q = qr.Q(qr(X)) $Qz$ Q %*% z

Crucialement, la multiplication de par les valeurs propres de ne change pas l'orthogonalité des vecteurs de colonne constitutifs, mais étant donné que la valeur absolue des valeurs propres apparaît dans l'ordre décroissant du haut à gauche vers le bas à droite, la multiplication de aura tendance à diminuer la valeurs dans les colonnes polynomiales d'ordre supérieur: $\bf Q$ $\bf R$ $Qz$

Matrix of Eigenvalues of R
     [,1]      [,2] [,3]      [,4]
[1,]   -2  0.000000    0  0.000000
[2,]    0 -2.236068    0  0.000000
[3,]    0  0.000000    2  0.000000
[4,]    0  0.000000    0 -1.341641

Comparez les valeurs des derniers vecteurs de colonne (quadratique et cubique) avant et après les opérations de factorisation , et aux deux premières colonnes non affectées. $QR$

Before QR factorization operations (orthogonal col. vec.)
     [,1] [,2] [,3]   [,4]
[1,]    1 -1.5 2.25 -3.375
[2,]    1 -0.5 0.25 -0.125
[3,]    1  0.5 0.25  0.125
[4,]    1  1.5 2.25  3.375


After QR operations (equally orthogonal col. vec.)
     [,1] [,2] [,3]   [,4]
[1,]    1 -1.5    1 -0.295
[2,]    1 -0.5   -1  0.885
[3,]    1  0.5   -1 -0.885
[4,]    1  1.5    1  0.295

Enfin, nous appelons la (Z <- sweep(raw, 2L, apply(raw, 2L, function(x) sqrt(sum(x^2))), "/", check.margin = FALSE))transformation de la matrice rawen vecteurs orthonormés :

Orthonormal vectors (orthonormal basis of R^4)
     [,1]       [,2] [,3]       [,4]
[1,]  0.5 -0.6708204  0.5 -0.2236068
[2,]  0.5 -0.2236068 -0.5  0.6708204
[3,]  0.5  0.2236068 -0.5 -0.6708204
[4,]  0.5  0.6708204  0.5  0.2236068

Cette fonction "normalise" simplement la matrice en divisant ( "/") dans les colonnes chaque élément par le . Il peut donc être décomposé en deux étapes:, cequi donne les dénominateurs pour chaque colonne deoù chaque élément d'une colonne est divisé par la valeur correspondante de. $\small\sqrt{\sum_\text{col.} x_i^2}$ $(\text{i})$ apply(raw, 2, function(x)sqrt(sum(x^2)))2 2.236 2 1.341 $(\text{ii})$ $(\text{i})$

À ce stade, les vecteurs de colonne forment une base orthonormée de , jusqu'à ce que nous nous débarrassions de la première colonne, qui sera l'ordonnée à l'origine, et nous avons reproduit le résultat de : $\mathbb{R}^4$ contr.poly(4)

$\small\begin{bmatrix} -0.6708204&0.5&-0.2236068\\-0.2236068&-0.5&0.6708204\\0.2236068&-0.5&-0.6708204\\0.6708204&0.5&0.2236068 \end{bmatrix}$

Les colonnes de cette matrice sont orthonormées , comme le montrent (sum(Z[,3]^2))^(1/4) = 1et z[,3]%*%z[,4] = 0par exemple (incidemment, il en va de même pour les lignes). Et, chaque colonne est le résultat de l' élévation des premiers pour le -st, -ND et alimentation -rd, respectivement - c. -à- linéaire, quadratique et cubique . $\text{scores - mean}$ $1$ $2$ $3$

2. Quels contrastes (colonnes) contribuent de manière significative à expliquer les différences entre les niveaux de la variable explicative?

Nous pouvons simplement exécuter l'ANOVA et regarder le résumé ...

summary(lm(write ~ readcat, hsb2))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  52.7870     0.6339  83.268   <2e-16 ***
readcat.L    14.2587     1.4841   9.607   <2e-16 ***
readcat.Q    -0.9680     1.2679  -0.764    0.446    
readcat.C    -0.1554     1.0062  -0.154    0.877

... pour voir qu'il y a un effet linéaire de readcaton write, de sorte que les valeurs d'origine (dans le troisième morceau de code au début de l'article) peuvent être reproduites comme:

coeff = coefficients(lm(write ~ readcat, hsb2))
C = contr.poly(4)
(recovered = c(coeff %*% c(1, C[1,]),
               coeff %*% c(1, C[2,]),
               coeff %*% c(1, C[3,]),
               coeff %*% c(1, C[4,])))
[1] 42.77273 49.97849 56.56364 61.83333

... ou...

... ou bien mieux ...

Étant des contrastes orthogonaux, la somme de leurs composantes ajoute à zéro pour constantes, et le produit scalaire de deux d'entre elles est nul. Si nous pouvions les visualiser, ils ressembleraient à ceci: $\displaystyle \sum_{i=1}^t a_i = 0$ $a_1,\cdots,a_t$

$X^0, X^1, \cdots. X^n$

Graphiquement, c'est beaucoup plus facile à comprendre. Comparez les moyennes réelles par groupes dans de grands blocs noirs carrés aux valeurs prédites, et voyez pourquoi une approximation en ligne droite avec une contribution minimale des polynômes quadratiques et cubiques (avec des courbes uniquement approximées avec du loess) est optimale:

Si, juste pour l'effet, les coefficients de l'ANOVA avaient été aussi importants pour le contraste linéaire pour les autres approximations (quadratique et cubique), le tracé absurde qui suit représenterait plus clairement les diagrammes polynomiaux de chaque "contribution":

Le code est ici .

Antoni Parellada
la source

+1 Wow. Cette réponse (je ne l'ai pas lue jusqu'à la fin jusqu'à présent) peut-elle également être considérée comme une réponse à ma vieille question oubliée stats.stackexchange.com/q/63639/3277 ?

ttnphns

(+1) @ttnphns: On peut dire que ça irait encore mieux là-bas.

Scortchi - Réintégrer Monica

Juste un conseil: vous voudrez peut-être me commenter là-bas avec un lien vers ici; ou émettre une réponse - que je suis susceptible d'accepter.

ttnphns

@ttnphns et @Scortchi Merci! J'ai passé pas mal de temps à essayer de donner un sens à ces concepts et je ne m'attendais pas à beaucoup de réactions. C'est donc une surprise très positive. Je pense qu'il y a quelques rides à corriger en ce qui concerne l'explication de la qr.qy()fonction, mais je vais certainement essayer de voir si je peux dire quelque chose de peu cohérent à propos de votre question dès que j'aurai du temps.

Antoni Parellada

@Elvis J'ai essayé de choisir une bonne phrase récapitulative et de la placer quelque part dans le message. Je pense que c'est un bon point, et appelle une belle explication mathématique, mais il peut être trop à ce stade de développer davantage.

Antoni Parellada

Je vais utiliser votre exemple pour expliquer comment cela fonctionne. L'utilisation de contrastes polynomiaux avec quatre groupes donne les résultats suivants.

\begin{aligned} E w r je t e_{1} & = μ - 0,67 L + 0,5 Q - 0,22 C \\ E w r je t e_{2} & = μ - 0,22 L - 0,5 Q + 0,67 C \\ E w r je t e_{3} & = μ + 0,22 L - 0,5 Q - 0,67 C \\ E w r je t e_{4} & = μ + 0,67 L + 0,5 Q + 0,22 C \end{aligned}

$\begin{align} E\,write_1 &= \mu -0.67L + 0.5Q -0.22C\\ E\,write_2 &= \mu -0.22L -0.5Q + 0.67C\\ E\,write_3 &= \mu + 0.22L -0.5Q -0.67C\\ E\,write_4 &= \mu + 0.67L + 0.5Q + 0.22C \end{align}$

Où la première équation fonctionne pour le groupe des scores de lecture les plus bas et la quatrième pour le groupe des meilleurs scores de lecture. nous pouvons comparer ces équations à celle donnée en utilisant une régression linéaire normale (en supposant $read_i$ est continu)

E w r je t e_{je} = μ + r e une {ré}_{je} L + r e une {ré}_{je}^{2} Q + r e une {ré}_{je}^{3} C

$E\,write_i=\mu+read_iL + read_i^2Q+read_i^3C$

Habituellement au lieu de $L,Q,C$ tu aurais $\beta_1, \beta_2, \beta_3$ et écrit en première position. Mais cette écriture ressemble à celle aux contrastes polynomiaux. Donc, les chiffres devant $L, Q, C$ sont en fait au lieu de $read_i, read_i^2, read_i^3$ . You can see that coefficients before $L$ have linear trend, before $Q$ quadratic and before $C$ cubic.

Then R estimates parameters $\mu, L,Q,C$ and gives you

\hat{μ} = 52.79, \hat{L} = 14.26, \hat{Q} = - 0.97, \hat{C} = - 0.16

$\widehat{\mu}=52.79, \widehat{L}=14.26, \widehat{Q}=−0.97, \widehat{C}=−0.16$ Where

\hat{μ} = \frac{1}{4} \sum_{i = 1}^{4} E w r i t e_{i}

$\widehat{\mu}=\frac{1}{4}\sum_{i=1}^4E\,write_i$ and estimated coefficients

\hat{μ}, \hat{L}, \hat{Q}, \hat{C}

$\widehat{\mu}, \widehat{L}, \widehat{Q}, \widehat{C}$ are something like estimates at normal linear regression. So from the output you can see if estimated coefficients are significantly different from zero, so you could anticipate some kind of linear, quadratic or cubic trend.

In that example is significantly non-zero only $\widehat{L}$ . So your conclusion could be: We see that the better scoring in writing depends linearly on reading score, but there is no significant quadratic or cubic effect.

Fimba
la source

Contrastes polynomiaux pour la régression

Réponses: