L'ANCOVA dans R suggère différentes interceptions, mais les IC à 95% se chevauchent… comment est-ce possible?

Nous avons un ensemble de données avec deux covariables et une variable de regroupement catégorique et nous voulons savoir s'il existe des différences significatives entre la pente ou l'ordonnée à l'origine parmi les covariables associées aux différentes variables de regroupement. Nous avons utilisé anova () et lm () pour comparer les ajustements de trois modèles différents: 1) avec une seule pente et interception, 2) avec différentes interceptions pour chaque groupe, et 3) avec une pente et une interception pour chaque groupe . Selon le test linéaire général anova (), le deuxième modèle est le plus approprié des trois, il y a une amélioration significative du modèle en incluant une interception distincte pour chaque groupe. Cependant, lorsque nous examinons les intervalles de confiance à 95% pour ces interceptions - ils se chevauchent tous, ce qui suggère qu'il n'y a pas de différences significatives entre les interceptions. Comment concilier ces deux résultats? Nous avons pensé qu'une autre façon d'interpréter les résultats de la méthode de sélection de modèle était qu'il devait y avoir au moins une différence significative entre les interceptions ... mais peut-être que ce n'est pas correct?

Vous trouverez ci-dessous le code R pour reproduire cette analyse. Nous avons utilisé la fonction dput () pour que vous puissiez travailler avec exactement les mêmes données que celles avec lesquelles nous nous débattons.

# Begin R Script
# > dput(data)
structure(list(Head = c(1.92, 1.93, 1.79, 1.94, 1.91, 1.88, 1.91, 
1.9, 1.97, 1.97, 1.95, 1.93, 1.95, 2, 1.87, 1.88, 1.97, 1.88, 
1.89, 1.86, 1.86, 1.97, 2.02, 2.04, 1.9, 1.83, 1.95, 1.87, 1.93, 
1.94, 1.91, 1.96, 1.89, 1.87, 1.95, 1.86, 2.03, 1.88, 1.98, 1.97, 
1.86, 2.04, 1.86, 1.92, 1.98, 1.86, 1.83, 1.93, 1.9, 1.97, 1.92, 
2.04, 1.92, 1.9, 1.93, 1.96, 1.91, 2.01, 1.97, 1.96, 1.76, 1.84, 
1.92, 1.96, 1.87, 2.1, 2.17, 2.1, 2.11, 2.17, 2.12, 2.06, 2.06, 
2.1, 2.05, 2.07, 2.2, 2.14, 2.02, 2.08, 2.16, 2.11, 2.29, 2.08, 
2.04, 2.12, 2.02, 2.22, 2.22, 2.2, 2.26, 2.15, 2, 2.24, 2.18, 
2.07, 2.06, 2.18, 2.14, 2.13, 2.2, 2.1, 2.13, 2.15, 2.25, 2.14, 
2.07, 1.98, 2.16, 2.11, 2.21, 2.18, 2.13, 2.06, 2.21, 2.08, 1.88, 
1.81, 1.87, 1.88, 1.87, 1.79, 1.99, 1.87, 1.95, 1.91, 1.99, 1.85, 
2.03, 1.88, 1.88, 1.87, 1.85, 1.94, 1.98, 2.01, 1.82, 1.85, 1.75, 
1.95, 1.92, 1.91, 1.98, 1.92, 1.96, 1.9, 1.86, 1.97, 2.06, 1.86, 
1.91, 2.01, 1.73, 1.97, 1.94, 1.81, 1.86, 1.99, 1.96, 1.94, 1.85, 
1.91, 1.96, 1.9, 1.98, 1.89, 1.88, 1.95, 1.9, 1.94, NA, 1.84, 
1.83, 1.84, 1.96, 1.74, 1.91, 1.84, 1.88, 1.83, 1.93, 1.78, 1.88, 
1.93, 2.15, 2.16, 2.23, 2.09, 2.36, 2.31, 2.25, 2.29, 2.3, 2.04, 
2.22, 2.19, 2.25, 2.31, 2.3, 2.28, 2.25, 2.15, 2.29, 2.24, 2.34, 
2.2, 2.24, 2.17, 2.26, 2.18, 2.17, 2.34, 2.23, 2.36, 2.31, 2.13, 
2.2, 2.27, 2.27, 2.2, 2.34, 2.12, 2.26, 2.18, 2.31, 2.24, 2.26, 
2.15, 2.29, 2.14, 2.25, 2.31, 2.13, 2.09, 2.24, 2.26, 2.26, 2.21, 
2.25, 2.29, 2.15, 2.2, 2.18, 2.16, 2.14, 2.26, 2.22, 2.12, 2.12, 
2.16, 2.27, 2.17, 2.27, 2.17, 2.3, 2.25, 2.17, 2.27, 2.06, 2.13, 
2.11, 2.11, 1.97, 2.09, 2.06, 2.11, 2.09, 2.08, 2.17, 2.12, 2.13, 
1.99, 2.08, 2.01, 1.97, 1.97, 2.09, 1.94, 2.06, 2.09, 2.04, 2, 
2.14, 2.07, 1.98, 2, 2.19, 2.12, 2.06, 2, 2.02, 2.16, 2.1, 1.97, 
1.97, 2.1, 2.02, 1.99, 2.13, 2.05, 2.05, 2.16, 2.02, 2.02, 2.08, 
1.98, 2.04, 2.02, 2.07, 2.02, 2.02, 2.02), Site = structure(c(2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("ANZ", "BC", "DV", "MC", 
"RB", "WW"), class = "factor"), Leg = c(2.38, 2.45, 2.22, 2.23, 
2.26, 2.32, 2.28, 2.17, 2.39, 2.27, 2.42, 2.33, 2.31, 2.32, 2.25, 
2.27, 2.38, 2.28, 2.33, 2.24, 2.21, 2.22, 2.42, 2.23, 2.36, 2.2, 
2.28, 2.23, 2.33, 2.35, 2.36, 2.26, 2.26, 2.3, 2.23, 2.31, 2.27, 
2.23, 2.37, 2.27, 2.26, 2.3, 2.33, 2.34, 2.27, 2.4, 2.22, 2.25, 
2.28, 2.33, 2.26, 2.32, 2.29, 2.31, 2.37, 2.24, 2.26, 2.36, 2.32, 
2.32, 2.15, 2.2, 2.29, 2.37, 2.26, 2.24, 2.23, 2.24, 2.26, 2.18, 
2.11, 2.23, 2.31, 2.25, 2.15, 2.3, 2.33, 2.35, 2.21, 2.36, 2.27, 
2.24, 2.35, 2.24, 2.33, 2.32, 2.24, 2.35, 2.36, 2.39, 2.28, 2.36, 
2.19, 2.27, 2.39, 2.23, 2.29, 2.32, 2.3, 2.32, NA, 2.25, 2.24, 
2.21, 2.37, 2.21, 2.21, 2.27, 2.27, 2.26, 2.19, 2.2, 2.25, 2.25, 
2.25, NA, 2.24, 2.17, 2.2, 2.2, 2.18, 2.14, 2.17, 2.27, 2.28, 
2.27, 2.29, 2.23, 2.25, 2.33, 2.22, 2.29, 2.19, 2.15, 2.24, 2.24, 
2.26, 2.25, 2.09, 2.27, 2.18, 2.2, 2.25, 2.24, 2.18, 2.3, 2.26, 
2.18, 2.27, 2.12, 2.18, 2.33, 2.13, 2.28, 2.23, 2.16, 2.2, 2.3, 
2.31, 2.18, 2.33, 2.29, 2.26, 2.21, 2.22, 2.27, 2.32, 2.24, 2.25, 
2.17, 2.2, 2.26, 2.27, 2.24, 2.25, 2.09, 2.25, 2.21, 2.24, 2.21, 
2.22, 2.13, 2.24, 2.21, 2.3, 2.34, 2.35, 2.32, 2.46, 2.43, 2.42, 
2.41, 2.32, 2.25, 2.33, 2.19, 2.45, 2.32, 2.4, 2.38, 2.35, 2.39, 
2.29, 2.35, 2.43, 2.29, 2.33, 2.31, 2.28, 2.38, 2.32, 2.43, 2.27, 
2.4, 2.37, 2.27, 2.41, 2.32, 2.38, 2.23, 2.33, 2.21, 2.34, 2.19, 
2.34, 2.35, 2.35, 2.31, 2.33, 2.41, 2.53, 2.39, 2.17, 2.16, 2.38, 
2.34, 2.33, 2.33, 2.29, 2.43, 2.28, 2.34, 2.38, 2.3, 2.29, 2.43, 
2.36, 2.24, 2.35, 2.38, 2.4, 2.36, 2.42, 2.28, 2.45, 2.33, 2.32, 
2.33, 2.31, 2.44, 2.37, 2.4, 2.35, 2.33, 2.31, 2.36, 2.43, 2.38, 
2.4, 2.38, 2.46, 2.33, 2.38, 2.23, 2.24, 2.39, 2.36, 2.19, 2.32, 
2.37, 2.39, 2.34, 2.39, 2.23, 2.25, 2.29, 2.39, 2.35, NA, 2.28, 
2.35, 2.38, 2.34, 2.17, 2.29, NA, 2.26, NA, NA, NA, 2.24, 2.33, 
2.23, 2.28, 2.29, 2.23, 2.2, 2.27, 2.31, 2.31, 2.26, 2.28)), .Names = c("Head", 
"Site", "Leg"), class = "data.frame", row.names = c(NA, -312L
)) 

# plot graph
library(ggplot2)

qplot(Head, Leg, 
    color=Site, 
    data=data) + 
        stat_smooth(method="lm", alpha=0.2) + 
        theme_bw()

# create linear models
lm.1 <- lm(Leg ~ Head, data)
lm.2 <- lm(Leg ~ Head + Site, data)
lm.3 <- lm(Leg ~ Head*Site, data)

# evaluate linear models
anova(lm.1, lm.2, lm.3)
anova(lm.1, lm.2)

# > anova(lm.1, lm.2)
# Analysis of Variance Table
# Model 1: Leg.3.1 ~ Head.W1
# Model 2: Leg.3.1 ~ Head.W1 + Site
  # Res.Df     RSS Df Sum of Sq     F    Pr(>F)    
# 1    302 1.25589                                 
# 2    297 0.91332  5   0.34257 22.28 < 2.2e-16 ***


# examining the multiple-intercepts model (lm.2)
summary(lm.2)
coef(lm.2)
confint(lm.2)

# extracting the intercepts
intercepts <- coef(lm.2)[c(1, 3:7)]
intercepts.1 <- intercepts[1]
intercepts <- intercepts.1 + intercepts
intercepts[1] <- intercepts.1
intercepts

# extracting the confidence intervals
ci <- confint(lm.2)[c(1, 3:7),]
ci[2:6,] <- ci[2:6,] + confint(lm.2)[1,]
ci[,1]

# putting everything together in a dataframe
labels <- c("ANZ", "BC", "DV", "MC", "RB", "WW")
ci.dataframe <- data.frame(Site=labels, Intercept=intercepts, CI.low = ci[,1], CI.high = ci[,2])
ci.dataframe

# plotting intercepts and 95% CI
qplot(Site, Intercept, geom=c("point", "errorbar"), ymin=CI.low, ymax=CI.high, data=ci.dataframe, ylab="Intercept & 95% CI")

Pour résumer - le problème est que les IC à 95% pour les interceptions se chevauchent tous, mais la méthode de sélection du modèle suggère que le meilleur modèle est celui qui correspond à différentes interceptions. Je suis donc enclin à penser que notre méthode de sélection de modèle est défectueuse ou que les IC à 95% pour les estimations d'interception ont été calculés incorrectement. Toutes les pensées seraient grandement appréciées!

N'oubliez pas que la différence entre significatif et non significatif n'est pas (toujours) statistiquement significative

Maintenant, plus au point de votre question, le modèle 1 est appelé régression groupée et le modèle 2 régression non groupée. Comme vous l'avez noté, dans la régression groupée, vous supposez que les groupes ne sont pas pertinents, ce qui signifie que la variance entre les groupes est définie sur zéro.

Dans la régression non groupée, avec une interception par groupe, vous définissez la variance à l'infini.

En général, je préférerais une solution intermédiaire, qui est un modèle hiérarchique ou une régression groupée partielle (ou un estimateur de retrait). Vous pouvez adapter ce modèle en R avec le package lmer4.

Enfin, jetez un œil à cet article de Gelman , dans lequel il explique pourquoi les modèles hiérarchiques aident à résoudre les problèmes de comparaisons multiples (dans votre cas, les coefficients par groupe sont-ils différents? Comment corriger une valeur de p pour les comparaisons multiples).

Par exemple, dans votre cas,

library(lme4)
summary(lmer( leg ~ head + (1 | site)) # varying intercept model

Si vous souhaitez adapter une pente variable à l'interception variable (le troisième modèle), exécutez simplement

summary(lmer( leg ~ head + (1 | site) + (0+head|site) )) # varying intercept, varying-slope model

Ensuite, vous pouvez regarder la variance de groupe et voir si elle est différente de zéro (la régression groupée n'est pas le meilleur modèle) et loin de l'infini (régression non groupée).

mise à jour: Après les commentaires (voir ci-dessous), j'ai décidé d'élargir ma réponse.

Le but d'un modèle hiérarchique, spécialement dans des cas comme celui-ci, est de modéliser la variation par groupes (dans ce cas, Sites). Donc, au lieu d'exécuter une ANOVA pour tester si un modèle est différent d'un autre, je regarderais les prédictions de mon modèle et verrais si les prédictions par groupe sont meilleures dans les modèles hiérarchiques vs la régression groupée (régression classique) .

Maintenant, j'ai couru mes suggestions ci-dessus et foudn que

ranef(lmer( leg ~ head + (1 | site) + (0+head|site) )

Renverrait zéro comme estimation de la pente variable (effet variable de la hauteur par site). alors j'ai couru

ranef(lmer( leg ~ head + (head| site))

Et j'ai obtenu une estimation non nulle de l'effet variable de la tête. Je ne sais pas encore pourquoi cela s'est produit, car c'est la première fois que je trouve cela. Je suis vraiment désolé pour ce problème, mais, pour ma défense, je viens de suivre les spécifications décrites dans l'aide de la fonction lmer. (Voir l'exemple avec l'étude de sommeil des données). J'essaierai de comprendre ce qui se passe et je ferai rapport ici quand (si) je comprendrai ce qui se passe.

Avant toute intervention du modérateur, vous pouvez regarder

library(car)

crPlots(lm.2,terms=~Site)

Ce sont des parcelles Composant + Résiduel (Résiduel Partiel)

Je pense, entre autres, que vous calculez mal les intervalles de confiance. Voici deux façons de voir les choses:

(1) différences de chaque site par rapport au site de référence (ANZ) [vous pouvez également calculer les différences à partir de la moyenne globale en passant à des contrastes somme-à-zéro

library(coefplot2)  ## on r-forge
coefplot2(lm.2)

ou (2) toutes les comparaisons par paires (je n'aime pas cette approche, mais c'est courant):

library(multcomp)
ci <- confint(glht(lm.2, linfct = mcp(Site = "Tukey")))
ggplot(fortify(ci),aes(lhs,estimate,ymin=lwr,ymax=upr))+
    geom_pointrange()+theme_bw()+geom_hline(yintercept=0,col="red")

Notez que toutes vos Headvaleurs sont comprises entre 1,7 et 2,4, tandis que les interceptions tentent d'estimer la Legvaleur à Head=0. Il s'agit d'une extrapolation majeure, il y a donc beaucoup d'incertitude. Si vous deviez centrer les Headvaleurs et répéter cette analyse, les intervalles de confiance deviendraient beaucoup plus serrés.

De plus, le chevauchement des intervalles de confiance à 95% n'implique pas l'absence de différence statistiquement significative. En fait, pour deux groupes, les intervalles de confiance à 84% non chevauchants se rapprochent et présentent des différences significatives au niveau de 5%. Bien sûr, en raison de plusieurs tests, cela ne fonctionne pas tout à fait avec plusieurs groupes.

En plus des autres réponses, voici quelques liens de la Cornell Statistical Consulting Unit qui sont pertinents pour les intervalles de confiance qui se chevauchent et qui servent de bon et bref rappel de ce qu'ils font et ne signifient pas

http://www.cscu.cornell.edu/news/statnews/stnews73.pdf http://www.cscu.cornell.edu/news/statnews/Stnews73insert.pdf

Voici le point principal:

Si deux statistiques ont des intervalles de confiance qui ne se chevauchent pas, ils sont nécessairement significativement différents mais s'ils ont des intervalles de confiance qui se chevauchent, il n'est pas nécessairement vrai qu'ils ne sont pas significativement différents.

Voici le texte pertinent du premier lien:

Nous pouvons illustrer cela avec un exemple simple. Supposons que nous souhaitons comparer les moyennes de deux échantillons indépendants. La moyenne du premier échantillon est 9 et la moyenne du deuxième échantillon est 17. Supposons que les moyennes des deux groupes ont les mêmes erreurs standard, égales à 2,5. L'intervalle de confiance à 95% pour la moyenne du premier groupe peut être calculé comme suit: ± × 5.296.19 où 1.96 est la valeur t critique. L'intervalle de confiance pour la moyenne du premier groupe est donc (4.1, 13.9). De même pour le deuxième groupe, l'intervalle de confiance pour la moyenne est (12,1, 21,9). Notez que les deux intervalles se chevauchent. Cependant, la statistique t pour comparer deux moyennes est:

t = (17-9)/√(2.5² + 2.5²) = 2.26

ce qui reflète que l'hypothèse nulle, selon laquelle les moyennes des deux groupes sont les mêmes, doit être rejetée au niveau α = 0,05. Pour vérifier la conclusion ci-dessus, considérons l'intervalle de confiance à 95% pour la différence entre les moyennes des deux groupes: (17-9) ± 1,96 x √ (2,5² + 2,5²) qui donne (1,09, 14,91). L'intervalle ne contient pas zéro, donc nous rejetons l'hypothèse nulle que les moyennes de groupe sont les mêmes.

Généralement, lorsque l'on compare deux estimations de paramètres, il est toujours vrai que si les intervalles de confiance ne se chevauchent pas, les statistiques seront statistiquement significativement différentes. Cependant, l'inverse n'est pas vrai. En d'autres termes, il est erroné de déterminer la signification statistique de la différence entre deux statistiques sur la base d'intervalles de confiance qui se chevauchent. Pour une explication de la raison pour laquelle cela est vrai dans le cas de la comparaison des moyennes sur deux échantillons, voir le lien suivant: http://www.cscu.cornell.edu/news/statnews/Stnews73insert.pdf

Voici les informations de l'autre lien: