Est-il tout à fait défendable de stratifier un ensemble de données en fonction de la taille du résidu et de faire une comparaison à deux échantillons?

C'est quelque chose que je vois comme une sorte de méthode ad hoc et cela me semble très louche, mais peut-être que je manque quelque chose. J'ai vu cela en régression multiple, mais restons simples:

Maintenant, prenez les résidus du modèle ajusté

et stratifier l'échantillon en fonction de la taille des résidus. Par exemple, supposons que le premier échantillon soit les 90% inférieurs des résidus et le deuxième échantillon les 10% supérieurs, puis procédez à deux comparaisons d'échantillons - j'ai vu cela fait à la fois sur le prédicteur du modèle, $x$ et sur des variables absentes du modèle. La logique informelle utilisée est que peut-être les points qui ont des valeurs bien supérieures à ce que vous attendez du modèle (c'est-à-dire un grand résidu) sont différents d'une certaine manière, et cette différence est étudiée de cette manière.

Mes réflexions sur le sujet sont:

• Si vous voyez une différence de 2 échantillons sur un prédicteur dans le modèle, alors il y a des effets du prédicteur non pris en compte par le modèle dans son état actuel (c'est-à-dire des effets non linéaires).
• Si vous voyez une différence de 2 échantillons sur une variable qui n'est pas dans le modèle, alors peut-être qu'elle aurait dû être dans le modèle en premier lieu.

Une chose que j'ai trouvée empiriquement (à travers des simulations) est que, si vous comparez la moyenne d'un prédicteur dans le modèle et stratifiez de cette façon pour produire les deux moyennes d'échantillon, ¯ x 1 et ¯ x 2 , elles sont corrélé positivement les uns avec les autres. Cela est logique puisque les deux échantillons dépendent de la ¯ y , ¯ x , σ x , σ y et ρ x $x$$\overline{x}_{1}$$\overline{x}_{2}$$\overline{y}, \overline{x}, \hat{\sigma}_{x}, \hat{\sigma}_{y}$$\hat{\rho}_{xy}$. Cette corrélation augmente à mesure que vous abaissez le seuil (c'est-à-dire le% que vous utilisez pour diviser l'échantillon). Donc, à tout le moins, si vous voulez faire une comparaison à deux échantillons, l'erreur standard dans le dénominateur de la statistique doit être ajustée pour tenir compte de la corrélation (bien que je n'ai pas dérivé de formule explicite pour la covariance).$t$

Quoi qu'il en soit, ma question fondamentale est la suivante: y a-t-il une justification à cela? Si oui, dans quelles situations cela pourrait-il être utile de le faire? De toute évidence, je ne pense pas qu'il y en ait, mais il y a peut-être quelque chose auquel je ne pense pas de la bonne façon.

Comparer les moyennes est trop faible: comparez plutôt les distributions.

Il y a aussi une question de savoir s'il est plus souhaitable de comparer les tailles of the residuals (as stated) or to compare the residuals themselves. Therefore, I evaluate both.

Pour être précis sur ce que l'on veut dire, voici du Rcode à comparer$(x,y)$données (données dans des tableaux parallèles xet y) en régressant$y$ sur $x$, divisant les résidus en trois groupes en les coupant en dessous du quantile$q_0$ et au-dessus du quantile $q_1\gt q_0$, et (au moyen d'un tracé qq) comparer les distributions de $x$ valeurs associées à ces deux groupes.

test <- function(y, x, q0, q1, abs0=abs, ...) {
  y.res <- abs0(residuals(lm(y~x)))
  y.groups <- cut(y.res, quantile(y.res, c(0,q0,q1,1)))
  x.groups <- split(x, y.groups)
  xy <- qqplot(x.groups[[1]], x.groups[[3]], plot.it=FALSE)
  lines(xy, xlab="Low residual", ylab="High residual", ...)
}

Le cinquième argument de cette fonction, abs0 utilise par défaut les tailles (valeurs absolues) des résidus pour former les groupes. Plus tard, nous pouvons remplacer cela par une fonction qui utilise les résidus eux-mêmes.

Les résidus sont utilisés pour détecter de nombreuses choses: valeurs aberrantes, corrélations possibles avec des variables exogènes, qualité de l'ajustement et homoscédasticité. Les valeurs aberrantes, de par leur nature, devraient être peu nombreuses et isolées et ne joueront donc pas un rôle significatif ici. Pour garder cette analyse simple, explorons les deux derniers: qualité de l'ajustement (c'est-à-dire, linéarité du$x$-$y$relation) et l’homoscédasticité (c’est-à-dire la constance de la taille des résidus). Nous pouvons le faire grâce à la simulation:

simulate <- function(n, beta0=0, beta1=1, beta2=0, sd=1, q0=1/3, q1=2/3, abs0=abs,
                     n.trials=99, ...) {
  x <- 1:n - (n+1)/2
  y <- beta0 + beta1 * x + beta2 * x^2 + rnorm(n, sd=sd)
  plot(x,y, ylab="y", cex=0.8, pch=19, ...)
  plot(x, res <- residuals(lm(y ~ x)), cex=0.8, col="Gray", ylab="", main="Residuals")
  res.abs <- abs0(res)
  r0 <- quantile(res.abs, q0); r1 <- quantile(res.abs, q1)
  points(x[res.abs < r0], res[res.abs < r0], col="Blue")
  points(x[res.abs > r1], res[res.abs > r1], col="Red")
  plot(x,x, main="QQ Plot of X",
       xlab="Low residual", ylab="High residual",
       type="n")
  abline(0,1, col="Red", lwd=2)
  temp <- replicate(n.trials, test(beta0 + beta1 * x + beta2 * x^2 + rnorm(n, sd=sd), 
                             x, q0=q0, q1=q1, abs0=abs0, lwd=1.25, lty=3, col="Gray"))
  test(y, x, q0=q0, q1=q1, abs0=abs0, lwd=2, col="Black")
}

Ce code accepte des arguments déterminant le modèle linéaire: ses coefficients $y \sim \beta_0 + \beta_1 x + \beta_2 x^2$, the standard deviations of the error terms sd, the quantiles $q_0$ and $q_1$, the size function abs0, and the number of independent trials in the simulation, n.trials. The first argument n is the amount of data to simulate in each trial. It produces a set of plots--of the $(x,y)$ data, of their residuals, and qq plots of multiple trials--to help us understand how the proposed tests work for a given model (as determined by n, the beta,s and sd). Examples of these plots appear below.

Let us now use these tools to explore some realistic combinations of nonlinearity and heteroscedasticity, using the absolute values of the residuals:

n <- 100
beta0 <- 1
beta1 <- -1/n
sigma <- 1/n

size <- function(x) abs(x)
set.seed(17)
par(mfcol=c(3,4))
simulate(n, beta0, beta1, 0, sigma*sqrt(n), abs0=size, main="Linear Homoscedastic")
simulate(n, beta0, beta1, 0, 0.5*sigma*(n:1), abs0=size, main="Linear Heteroscedastic")
simulate(n, beta0, beta1, 1/n^2, sigma*sqrt(n), abs0=size, main="Quadratic Homoscedastic")
simulate(n, beta0, beta1, 1/n^2, 5*sigma*sqrt(1:n), abs0=size, main="Quadratic Heteroscedastic")

The output is a set of plots. The top row shows one simulated dataset, the second row shows a scatterplot of its residuals against $x$ (color-coded by quantile: red for large values, blue for small values, gray for any intermediate values not further used), and the third row shows the qq plots for all trials, with the qq plot for the one simulated dataset shown in black. An individual qq plot compares the $x$ valeurs associées à des résidus élevés à la $x$valeurs associées à de faibles résidus; après de nombreux essais, une enveloppe grise de parcelles qq probables émerge. Nous nous intéressons à la manière et à la force avec laquelle ces enveloppes varient en fonction des écarts par rapport au modèle linéaire de base: une forte variation implique une bonne discrimination.

Les différences entre les trois dernières et les premières colonnes montrent clairement que cette méthode est capable de détecter l'hétéroscédasticité, mais qu'elle pourrait ne pas être aussi efficace pour identifier une non-linéarité modérée. Cela pourrait facilement confondre la non-linéarité avec l'hétéroscédasticité. En effet, la forme d'hétéroscédasticité simulée ici (qui est courante) est celle où les tailles attendues de la tendance des résidus avec$x$. Cette tendance est facile à détecter. La non-linéarité quadratique, d'autre part, créera de grands résidus aux deux extrémités et au milieu de la plage de$x$valeurs. Ceci est difficile à distinguer simplement en regardant les distributions des$x$ valeurs.

Faisons la même chose, en utilisant exactement les mêmes données , mais en analysant les résidus eux-mêmes. Pour ce faire, le bloc de code précédent a été réexécuté après avoir effectué cette modification:

size <- function(x) x

Cette variation ne détecte pas bien l'hétéroscédasticité: voyez à quel point les graphes qq sont similaires dans les deux premières colonnes. Cependant, il détecte bien la non-linéarité. En effet, les résidus séparent les$x$est dans une partie médiane et une partie extérieure, qui seront assez différentes. Comme indiqué dans la colonne de droite, cependant, l'hétéroscédasticité peut masquer les non-linéarités.

Peut-être que la combinaison de ces deux techniques fonctionnerait. Ces simulations (et leurs variantes, que le lecteur intéressé peut exécuter à loisir) démontrent que ces techniques ne sont pas sans mérite.

En général, cependant, on est beaucoup mieux servi en examinant les résidus de manière standard. Pour le travail automatisé, des tests formels ont été développés pour détecter le genre de choses que nous recherchons dans les parcelles résiduelles. Par exemple, le test de Breusch-Pagan régresse les résidus au carré (plutôt que leurs valeurs absolues) contre$x$. Les tests proposés dans cette question peuvent être compris dans le même esprit. Cependant, en regroupant les données en deux groupes seulement et en négligeant ainsi la plupart des informations bivariées fournies par le$(x, \hat{y}-x)$paires, nous pouvons nous attendre à ce que les tests proposés soient moins puissants que les tests basés sur la régression comme le Breusch-Pagan .

Je suis d'accord avec vos deux points. Si le modèle est inadéquat, les résidus peuvent ne pas être approximativement indépendants et répartis de manière identique. Des variables importantes auraient pu être omises ou la forme fonctionnelle des variables régressives pourrait être incorrecte. Si tel est le cas, j'utiliserais les diagnostics de régression standard pour identifier le problème plutôt que cela. Vous pouvez également avoir les bonnes variables dans le modèle avec la bonne forme fonctionnelle mais avoir toujours une variance non constante. Cela pourrait être apparent simplement en traçant le$e_{i}$ contre $x_i$. Je peux voir un point à dire vouloir trouver des valeurs aberrantes dans le modèle à travers une certaine forme de résidu, mais je recommanderais alors une approche de fonction d'influence pour les détecter. Je ne vois pas ce que cette procédure accomplit.

D'autres ont indiqué qu'il ne s'agissait peut-être que d'un outil exploratoire pour voir si les deux ensembles de données devaient être modélisés séparément. Si tel est le cas, cela et peut-être d'autres approches exploratoires pourraient convenir. Mais la question devient alors que faites-vous ensuite? Si vous allez ensuite faire deux régressions distinctes et tirer des conclusions sur les échantillons, je pense que vous devez en quelque sorte tenir compte de la façon dont vous divisez l'échantillon.

Je suppose qu'il peut y avoir plusieurs motivations à le faire, par exemple en supposant que les résidus sont cohérents, la méthode que vous mentionnez peut aider à identifier les observations éloignées, donc la deuxième étape fournit des estimateurs «corrigés». Mais, il existe des techniques plus rigoureuses qui effectuent la détection des extérieurs ou qui fournissent des estimateurs robustes à la présence de telles observations, comme les régressions quantiles, le LMS (le moins médian des carrés) ou les estimateurs M, etc., où toutes ces méthodes ont bien défini et les propriétés statistiques connues. (Ceci a été adressé par @Michael Chernik)

Une autre motivation pourrait être l'identification de grappe, mais elle est primitive par rapport aux techniques disponibles pour la détection de grappe qui sont également bien définies et largement mises en œuvre.

Dans les deux cas, l'utilisation des résidus semble informelle et primitive, mais pourrait encore être tolérée comme un outil d'exploration. Cela dépend également du domaine des lecteurs. Je trouverais cela acceptable pour certaines sciences sociales où les outils quantitatifs peuvent être moins populaires.