Pourquoi certaines estimations de régression diffèrent-elles par un changement de signe, mais d'autres non, lorsque je change de niveau de référence?

8

Supposons que j'ai un résultat continu yet deux prédicteurs factoriels, chacun avec deux niveaux. L'un de mes prédicteurs catégoriques drug, peut avoir deux niveaux ("A" ou "B"), l'autre l'est smokeYes. Lorsque j'exécute un modèle de régression, je peux choisir la ligne de base ou le niveau de référence d' drugêtre soit "A", comme je l'ai fait dans model1:

set.seed(123)
y<-rnorm(100, 100, 10)
drug.ab<-factor(sample(c("A", "B"), 100, T), levels=c("A", "B"))
drug.ba<-factor(drug.ab, levels=c("B", "A"))
smoke<-factor(sample(c("Yes", "No"), 100, T), levels=c("No", "Yes"))

#model1:
coef(summary(lm(y~drug.ab*smoke)))
                     Estimate Std. Error    t value     Pr(>|t|)
(Intercept)       100.7484158   2.065091 48.7864379 1.465848e-69
drug.abB            0.9030541   2.796146  0.3229639 7.474250e-01
smokeYes           -0.8693598   2.632484 -0.3302431 7.419359e-01
drug.abB:smokeYes   0.8709116   3.746684  0.2324487 8.166844e-01

Ou je peux définir la ligne de base sur "B", comme je l'ai fait dans model2:

#model2:
coef(summary(lm(y~drug.ba*smoke)))
                       Estimate Std. Error       t value     Pr(>|t|)
(Intercept)       101.651469922   1.885161 53.9218978856 1.377147e-73
drug.baA           -0.903054145   2.796146 -0.3229638818 7.474250e-01
smokeYes            0.001551843   2.666021  0.0005820821 9.995368e-01
drug.baA:smokeYes  -0.870911601   3.746684 -0.2324486531 8.166844e-01

Ma question est pourquoi l'estimation de smokeYesdiffère entre model1et model2? Pourquoi ne diffère-t-il pas par un changement de signe drug.baAet le terme d'interaction?

regression anova multiple-regression contrasts David Z
la source

3

Recherchez une bonne explication du contraste du traitement. En bref, si vous calculez la prédiction pour drugB et smokeOui: (mod1) 100,75 + 0,90 - 0,87 + 0,87 = 101,65 | (mod2) 101,65 + 0,00 = 101,65

Roland

Je pensais qu'il était sans doute sur le sujet pour quand j'ai vu là la double question, car il y a une ligne de R très simple de code qui calcule tous les moyens de groupe: tapply( y, interaction( drug.ab, smoke) ,mean). Une explication plus approfondie pourrait impliquer de démontrer la différence entre les contrastes de traitement et les contrastes de somme.

DWin

@Dwin, même avec les deux réponses déjà publiées, je pense qu'il y a certainement de la place pour une autre réponse traitant précisément des problèmes de contraste. J'espère que quelqu'un publiera une réponse en adoptant cette approche.

Silverfish

8

Permettez-moi de vous donner un exemple simple pour expliquer le concept, puis nous pourrons le comparer à vos coefficients.

Notez qu'en incluant à la fois la variable fictive "A / B" et le terme d'interaction, vous donnez effectivement à votre modèle la flexibilité de s'adapter à une interception (en utilisant le modèle) et une pente (en utilisant l'interaction) différentes sur les données "A" et les données "B". Dans ce qui suit, peu importe si l'autre prédicteur $x$ est une variable continue ou, comme dans votre cas, une autre variable fictive. Si je parle en termes d '"interception" et de "pente", cela peut être interprété comme "niveau lorsque le mannequin est nul" et "changement de niveau lorsque le mannequin est changé de $0$ à $1$ " si tu préfères.

Supposons que le modèle ajusté OLS sur les seules données "A" soit $\hat y = 12 + 5x$ et sur les seules données "B" est $\hat y = 11 + 7x$ . Les données peuvent ressembler à ceci:

Supposons maintenant que nous prenons "A" comme niveau de référence et utilisons une variable fictive $b$ pour que $b=1$ pour les observations du groupe B mais $b=0$ dans le groupe A. Le modèle ajusté sur l'ensemble de données est

{\hat{y}}_{je} = {\hat{β}}_{0} + {\hat{β}}_{1} X_{je} + {\hat{β}}_{2} b_{je} + {\hat{β}}_{3} X_{je} b_{je}

$\hat y_i = \hat \beta_0 + \hat \beta_1 x_i + \hat \beta_2 b_i + \hat \beta_3 x_ib_i$

Pour les observations du groupe A, nous avons $\hat y_i = \hat \beta_0 + \hat \beta_1 x_i$ et nous pouvons minimiser leur somme des résidus carrés en définissant $\hat \beta_0 = 12$ et $\hat \beta_1 = 5$ . Pour les données du groupe B, $\hat y_i = (\hat \beta_0 + \hat \beta_2) + (\hat \beta_1 + \hat \beta_3) x_i$ et nous pouvons minimiser leur somme de résidus carrés en prenant $\hat \beta_0 + \hat \beta_2 = 11$ et $\hat \beta_1 + \hat \beta_3 = 7$ . Il est clair que nous pouvons minimiser la somme des résidus au carré dans la régression globale en minimisant les sommes pour les deux groupes, et que cela peut être réalisé en définissant $\hat \beta_0 = 12$ et $\hat \beta_1 = 5$ (du groupe A) et $\hat \beta_2 = -1$ et $\hat \beta_3 = 2$ (puisque les données "B" devraient avoir une interception inférieure et une pente supérieure). Observez comment la présence d'un terme d'interaction était nécessaire pour que nous ayons une flexibilité suffisante pour minimiser la somme des résidus au carré pour les deux groupes à la fois . Mon modèle ajusté sera:

{\hat{y}}_{je} = 12 + 5 X_{je} - 1 b_{je} + 2 X_{je} b_{je}

$\hat y_i = 12 + 5 x_i - 1 b_i +2 x_i b_i$

Basculez tout autour pour que "B" soit le niveau de référence et $a$ est une variable fictive codant pour le groupe A. Pouvez-vous voir que je dois maintenant adapter le modèle

{\hat{y}}_{je} = 11 + 7 X_{je} + 1 {une}_{je} - 2 X_{je} {une}_{je}

$\hat y_i = 11 + 7 x_i + 1 a_i -2 x_i a_i$

Autrement dit, je prends l'interception ( $11$ ) et pente ( $7$ ) de mon groupe de référence «B», et utilisez le terme fictif et d'interaction pour les ajuster pour mon groupe «A». Ces ajustements sont cette fois en sens inverse (j'ai besoin d'une interception plus haute et d'une pente deux plus bas ) donc les signes sont inversés par rapport à quand j'ai pris "A" comme groupe de référence, mais il devrait être clair pourquoi les autres coefficients ont pas simplement changé de signe.

Comparons cela à votre sortie. Dans une notation similaire à celle ci-dessus, votre premier modèle ajusté avec la ligne de base "A" est:

{\hat{y}}_{je} = 100.7484158 + 0.9030541 b_{je} - 0,8693598 X_{je} + 0,8709116 X_{je} b_{je}

$\hat y_i = 100.7484158 + 0.9030541 b_i -0.8693598 x_i + 0.8709116 x_i b_i$

Votre deuxième modèle équipé avec la ligne de base "B" est:

{\hat{y}}_{je} = 101,651469922 - 0.903054145 {une}_{je} + 0,001551843 X_{je} - 0,870911601 X_{je} {une}_{je}

$\hat y_i = 101.651469922 -0.903054145 a_i + 0.001551843 x_i -0.870911601 x_i a_i$

Tout d'abord, vérifions que ces deux modèles vont donner les mêmes résultats. Mettons $b_i = 1 - a_i$ dans la première équation, et on obtient:

{\hat{y}}_{je} = 100.7484158 + 0.9030541 (1 - {une}_{je}) - 0,8693598 X_{je} + 0,8709116 X_{je} (1 - {une}_{je})

$\hat y_i = 100.7484158 + 0.9030541 (1-a_i) -0.8693598 x_i + 0.8709116 x_i (1-a_i)$

Cela simplifie:

{\hat{y}}_{je} = (100.7484158 + 0.9030541) - 0.9030541 {une}_{je} + (- 0,8693598 + 0,8709116) X_{je} - 0,8709116 X_{je} {une}_{je}

$\hat y_i = (100.7484158 + 0.9030541) - 0.9030541 a_i + (-0.8693598 + 0.8709116) x_i - 0.8709116 x_i a_i$

Un peu d'arithmétique confirme que c'est la même chose que le deuxième modèle ajusté; en outre, il devrait maintenant être clair quels coefficients ont changé de signe et quels coefficients ont simplement été ajustés à l'autre référence!

Deuxièmement, voyons quels sont les différents modèles équipés sur les groupes "A" et "B". Votre premier modèle donne immédiatement $\hat y_i = 100.7484158 -0.8693598 x_i$ pour le groupe "A", et votre deuxième modèle donne immédiatement $\hat y_i = 101.651469922 + 0.001551843 x_i$ pour le groupe "B". Vous pouvez vérifier que le premier modèle donne le résultat correct pour le groupe "B" en remplaçant $b_i = 1$ dans son équation; l'algèbre, bien sûr, fonctionne de la même manière que l'exemple plus général ci-dessus. De même, vous pouvez vérifier que le deuxième modèle donne le résultat correct pour le groupe "A" en définissant $a_i = 1$ .

Troisièmement, étant donné que dans votre cas, l'autre régresseur était également une variable fictive, je vous suggère de calculer les moyennes conditionnelles ajustées pour les quatre catégories ("A" avec $x=0$ , "A" avec $x=1$ , "B" avec $x=0$ , "B" avec $x=1$ ) sous les deux modèles et vérifiez que vous comprenez pourquoi ils sont d'accord. À strictement parler, cela n'est pas nécessaire, car nous avons déjà effectué l'algèbre plus générale ci-dessus pour montrer que les résultats seront cohérents même si $x$ est continu , mais je pense que cela reste un exercice précieux. Je ne remplirai pas les détails car l'arithmétique est simple et plus conforme à l'esprit de la très belle réponse de JonB. Un point clé à comprendre est que, quel que soit le groupe de référence que vous utilisez, votre modèle a suffisamment de flexibilité pour s'adapter séparément à chaque moyenne conditionnelle. (C'est là que cela fait une différence que votre $x$ est un mannequin pour un facteur binaire plutôt que pour une variable continue - avec des prédicteurs continus, nous ne nous attendons généralement pas à la moyenne conditionnelle estimée $\hat y$ pour faire correspondre la moyenne de l'échantillon pour chaque combinaison de prédicteurs observée.) Calculez la moyenne de l'échantillon pour chacune de ces quatre combinaisons de catégories, et vous devriez trouver qu'elles correspondent à vos moyennes conditionnelles ajustées.

Code R pour tracer le tracé et explorer les modèles ajustés, prévu $\hat y$ et groupe signifie

#Make data set with desired conditional means
data.df <- data.frame(
  x = c(0,0,0,        1,1,1,        0,0,0,        1,1,1),
  b = c(0,0,0,        0,0,0,        1,1,1,        1,1,1),
  y = c(11.8,12,12.2, 16.8,17,17.2, 10.8,11,11.2, 17.8,18,18.2)
)
data.df$a <- 1 - data.df$b

baselineA.lm <- lm(y ~ x * b, data.df)
summary(baselineA.lm) #check this matches y = 12 + 5x - 1b + 2xb

baselineB.lm <- lm(y ~ x * a, data.df)
summary(baselineB.lm) #check this matches y = 11 + 7x + 1a - 2xa

fitted(baselineA.lm)
fitted(baselineB.lm) #check the two models give the same fitted values for y...
with(data.df, tapply(y, interaction(x, b), mean)) #...which are the group sample means

colorSet <- c("red", "blue")
symbolSet <- c(19,17)
with(data.df, plot(x, y, yaxt="n", col=colorSet[b+1], pch=symbolSet[b+1],
                   main="Response y against other predictor x",
                   panel.first = {
                     axis(2, at=10:20)
                     abline(h = 10:20, col="gray70")
                     abline(v = 0:1,  col="gray70")
                   }))
abline(lm(y ~ x, data.df[data.df$b==0,]), col=colorSet[1])
abline(lm(y ~ x, data.df[data.df$b==1,]), col=colorSet[2])
legend(0.1, 17, c("Group A", "Group B"), col = colorSet,
       pch = symbolSet, bg = "gray95")

Silverfish
la source

Oui, bonne explication donc je vais voter pour ça!

JonB

@DavidZ Merci! Je suggère de ne pas décocher la case "accepter" trop tôt, car il pourrait y avoir d'autres réponses à venir. Mon explication est assez soignée, mais j'ai abordé un aspect assez général qui aurait également fonctionné pour une

x

$x$ . Il est possible d'aborder votre question d'une manière qui accorde plus d'attention à la nature catégorique de vos prédicteurs, je vous suggère donc de décocher afin d'encourager plus de participation des autres. (@JobB Merci! J'ai aussi aimé votre réponse, +1)

Silverfish

3

Cela a à voir avec la façon dont l'interception est définie. Dans le premier exemple, l'interception est définie comme ceux qui ne fument pas et qui ont la drogue A. Les fumeurs, qui ont aussi la drogue A, auront une valeur de 100,75 - 0,87 = 99,9 tandis que les fumeurs qui ont la drogue B auront une valeur de 100,75 + 0,90 - 0,87 + 0,87 = 101,65.

Dans le deuxième exemple, l'interception est définie comme ceux qui ne fument pas et qui ont le médicament B. Les fumeurs avec le médicament B auront alors une valeur de 101,65 + 0,001 = 101,65, et les fumeurs avec le médicament A auront une valeur de 100,65 - 0,90 + 0,001-0,87 = 99,9.

Donc, tout cela ajoute upp, c'est juste une question de définition de l'interception, c'est-à-dire du niveau lorsque tous les facteurs sont définis sur la catégorie de référence.

JonB
la source

Pourquoi certaines estimations de régression diffèrent-elles par un changement de signe, mais d'autres non, lorsque je change de niveau de référence?

Réponses: