Quel test puis-je utiliser pour comparer les pentes de deux modèles de régression ou plus?

Je voudrais tester la différence de réponse de deux variables à un prédicteur. Voici un exemple reproductible minimal.

library(nlme) 
## gls is used in the application; lm would suffice for this example
m.set <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "setosa")
m.vir <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "virginica")
m.ver <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "versicolor")

Je peux voir que les coefficients de pente sont différents:

m.set$coefficients
(Intercept) Petal.Width 
  4.7771775   0.9301727
m.vir$coefficients
(Intercept) Petal.Width 
  5.2694172   0.6508306 
m.ver$coefficients
(Intercept) Petal.Width 
   4.044640    1.426365

J'ai trois questions:

Comment tester la différence entre les pentes?
Comment puis-je tester la différence entre les variances résiduelles?
Quelle est une manière simple et efficace de présenter ces comparaisons?

Une question connexe, Méthode pour comparer le coefficient variable dans deux modèles de régression , suggère de relancer le modèle avec une variable fictive pour différencier les pentes. Existe-t-il des options qui permettraient d'utiliser des ensembles de données indépendants?

r data-visualization multivariate-analysis hypothesis-testing Abe
la source

En ce qui concerne la première question, voir stats.stackexchange.com/questions/55501/… .

russellpierce

Réponses:

Comment tester la différence entre les pentes?

Incluez un mannequin pour les espèces, laissez-le interagir avec et voyez si ce mannequin est significatif. Soit la longueur du sépale et la largeur de la pédale et les variables muettes pour les trois espèces. Le comparer le modèle $P_i$ $L_i$ $P_i$ $S_1, S_2, S_3$

E (L_{i}) = β_{0} + β_{1} P_{i}

$E(L_i) = \beta_0 + \beta_1 P_i$

avec le modèle qui permet à l'effet de d'être différent pour chaque espèce: $P_i$

E (L_{i}) = α_{0} + α_{1} S_{2} + α_{2} S_{3} + α_{4} P_{i} + α_{5} P_{i} S_{2} + α_{6} P_{i} S_{3}

$E(L_i) = \alpha_0 + \alpha_1 S_2 + \alpha_2 S_3 + \alpha_4 P_i + \alpha_5 P_iS_2 + \alpha_6 P_i S_3$

Les estimateurs GLS sont des MLE et le premier modèle est un sous-modèle du second, vous pouvez donc utiliser le test du rapport de vraisemblance ici. Les probabilités peuvent être extraites à l'aide de la logLikfonction et les degrés de liberté pour le test seront de puisque vous avez supprimé paramètres pour arriver au sous-modèle. $4$ $4$

Quelle est une manière simple et efficace de présenter la comparaison?

Je pense que la façon la plus intéressante serait de tracer les lignes de régression pour chaque espèce sur les mêmes axes, peut-être avec des barres d'erreur basées sur les erreurs standard. Cela rendrait la différence (ou la non-différence) entre les espèces et leur relation avec très apparente. $P_i$

Edit: j'ai remarqué qu'une autre question a été ajoutée au corps. J'ajoute donc une réponse à cela:

Comment puis-je tester la différence entre les variances résiduelles?

Pour cela, vous devrez stratifier l'ensemble de données et ajuster des modèles distincts, car le modèle basé sur l'interaction que j'ai suggéré contraindra la variance résiduelle à être la même dans chaque groupe. Si vous ajustez des modèles distincts, cette contrainte disparaît. Dans ce cas, vous pouvez toujours utiliser le test du rapport de vraisemblance (la probabilité pour le modèle plus grand est maintenant calculée en additionnant les probabilités des trois modèles distincts). Le modèle "nul" dépend de ce avec quoi vous voulez le comparer

si vous voulez seulement tester la variance, tout en laissant les effets principaux dedans, alors le modèle "nul" devrait être le modèle avec les interactions que j'ai écrites ci-dessus. Les degrés de liberté pour le test sont alors de . $2$
Si vous voulez tester la variance conjointement avec les coefficients, alors le modèle nul devrait être le premier modèle que j'ai écrit ci-dessus. Le degré de liberté du test est alors de . $6$

Macro
la source

Pourquoi n'y a-t-il pas de dans le deuxième modèle? La mise en œuvre correcte du modèle dans R?

S_{1}

$S_1$ gls(Sepal.Length ~ species:Petal.Width, data = iris)

Abe

Salut @Abe. est l'espèce "de référence" - la droite de régression pour cette espèce est donnée par . Si est une variable catégorielle, je pense que ce serait la syntaxe.

S_{1}

$S_1$

α_{0} + α_{4} P_{i}

$\alpha_0 + \alpha_4 P_i$ speciesgls(Sepal.Length ~ species*Petal.Width, data=iris)

Macro

@Macro Belle réponse (+1)! Je me demande si vous pourriez adapter le glsmodèle mais en permettant différentes variances résiduelles pour chaque espèce avec l'option weights=varIdent(form=~1|Species)(concernant la deuxième question)?

COOLSerdash

Pour répondre à ces questions avec le code R, utilisez ce qui suit:
1. Comment puis-je tester la différence entre les pentes?
Réponse: Examinez la valeur p de l'ANOVA à partir de l'interaction de Petal.Width par espèce, puis comparez les pentes à l'aide de lsmeans :: lstrends, comme suit.

library(lsmeans)
m.interaction <- lm(Sepal.Length ~ Petal.Width*Species, data = iris)
anova(m.interaction)
 Analysis of Variance Table

 Response: Sepal.Length
                      Df Sum Sq Mean Sq  F value Pr(>F)    
 Petal.Width           1 68.353  68.353 298.0784 <2e-16 ***
 Species               2  0.035   0.017   0.0754 0.9274    
 Petal.Width:Species   2  0.759   0.380   1.6552 0.1947    
 Residuals           144 33.021   0.229                    
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Obtain slopes
m.interaction$coefficients
m.lst <- lstrends(m.interaction, "Species", var="Petal.Width")
 Species    Petal.Width.trend        SE  df   lower.CL upper.CL
 setosa             0.9301727 0.6491360 144 -0.3528933 2.213239
 versicolor         1.4263647 0.3459350 144  0.7425981 2.110131
 virginica          0.6508306 0.2490791 144  0.1585071 1.143154

# Compare slopes
pairs(m.lst)
 contrast                 estimate        SE  df t.ratio p.value
 setosa - versicolor    -0.4961919 0.7355601 144  -0.675  0.7786
 setosa - virginica      0.2793421 0.6952826 144   0.402  0.9149
 versicolor - virginica  0.7755341 0.4262762 144   1.819  0.1669

2. Comment puis-je tester la différence entre les variances résiduelles?
Si je comprends bien la question, vous pouvez comparer les corrélations de Pearson avec une transformée de Fisher, également appelée "r-to-z de Fisher", comme suit.

library(psych)
library(data.table)
iris <- as.data.table(iris)
# Calculate Pearson's R
m.correlations <- iris[, cor(Sepal.Length, Petal.Width), by = Species]
m.correlations
# Compare R values with Fisher's R to Z
paired.r(m.correlations[Species=="setosa", V1], m.correlations[Species=="versicolor", V1], 
         n = iris[Species %in% c("setosa", "versicolor"), .N])
paired.r(m.correlations[Species=="setosa", V1], m.correlations[Species=="virginica", V1], 
         n = iris[Species %in% c("setosa", "virginica"), .N])
paired.r(m.correlations[Species=="virginica", V1], m.correlations[Species=="versicolor", V1], 
         n = iris[Species %in% c("virginica", "versicolor"), .N])

3. Quelle est une manière simple et efficace de présenter ces comparaisons?
"Nous avons utilisé une régression linéaire pour comparer la relation entre la longueur du sépale et la largeur du pétale pour chaque espèce. Nous n'avons trouvé aucune interaction significative dans les relations entre la longueur du sépale et la largeur du pétale pour I. Setosa (B = 0,9), I. Versicolor (B = 1,4), ni I. Virginica (B = 0,6); F (2, 144) = 1,6, p = 0,19. Une comparaison r-à-z de Fisher a indiqué que la corrélation de Pearson pour I. Setosa (r = 0,28) était significativement plus faible (p = 0,02) que I. Versicolor (r = 0,55). De même, la corrélation pour I. Virginica (r = 0,28) était significativement plus faible (p = 0,02) que celle observée pour I. Versicolor. "

Enfin, visualisez toujours vos résultats!

plotly_interaction <- function(data, x, y, category, colors = col2rgb(viridis(nlevels(as.factor(data[[category]])))), ...) {
  # Create Plotly scatter plot of x vs y, with separate lines for each level of the categorical variable. 
  # In other words, create an interaction scatter plot.
  # The "colors" must be supplied in a RGB triplet, as produced by col2rgb().

  require(plotly)
  require(viridis)
  require(broom)

  groups <- unique(data[[category]])

  p <- plot_ly(...)

  for (i in 1:length(groups)) {
    groupData = data[which(data[[category]]==groups[[i]]), ]
    p <- add_lines(p, data = groupData,
                   y = fitted(lm(data = groupData, groupData[[y]] ~ groupData[[x]])),
                   x = groupData[[x]],
                   line = list(color = paste('rgb', '(', paste(colors[, i], collapse = ", "), ')')),
                   name = groups[[i]],
                   showlegend = FALSE)
    p <- add_ribbons(p, data = augment(lm(data = groupData, groupData[[y]] ~ groupData[[x]])),
                     y = groupData[[y]],
                     x = groupData[[x]],
                     ymin = ~.fitted - 1.96 * .se.fit,
                     ymax = ~.fitted + 1.96 * .se.fit,
                     line = list(color = paste('rgba','(', paste(colors[, i], collapse = ", "), ', 0.05)')), 
                     fillcolor = paste('rgba', '(', paste(colors[, i], collapse = ", "), ', 0.1)'),
                     showlegend = FALSE)
    p <- add_markers(p, data = groupData, 
                     x = groupData[[x]], 
                     y = groupData[[y]],
                     symbol = groupData[[category]],
                     marker = list(color=paste('rgb','(', paste(colors[, i], collapse = ", "))))
  }
  p <- layout(p, xaxis = list(title = x), yaxis = list(title = y))
  return(p)
}

plotly_interaction(iris, "Sepal.Length", "Petal.Width", "Species")

Kayle Sawyer
la source

Je suis d'accord avec la suggestion précédente. Vous devez adapter un modèle de régression multiple avec une variable fictive pour chaque ensemble de données. Cela vous permettra de tester si les interceptions diffèrent. Si vous souhaitez également savoir si les pentes diffèrent, vous devez également inclure les interactions entre les variables muettes et la variable en question. Il n'y a aucun problème avec le fait que les données sont indépendantes. Notez que s'ils sont à la fois des espèces indépendantes et (par exemple) différentes, vous ne pourrez pas dire si la différence que vous trouvez est due aux espèces ou aux ensembles de données différents, car ils sont parfaitement confondus. Cependant, il n'y a pas de test / carte de sortie de prison qui vous permettra de contourner ce problème sans rassembler un nouvel échantillon et relancer votre étude.

gung - Réintégrer Monica
la source

Il semble que nous ayons publié des réponses assez similaires presque en même temps. +1

Macro

@Macro, oui, mais le vôtre est généralement meilleur (+1 plus tôt); vous avez répondu aux 3 questions qui m'ont échappé lors de ma première lecture (non approfondie) de la question. Ma contribution ici est la partie sur la confusion.

gung - Réintégrer Monica

oui c'est un bon point. Je suppose que si vous faisiez cette enquête, vous devriez fonctionner en supposant que les ensembles de données mesuraient la même chose, etc ... à la seule différence que les espèces étaient différentes.

Macro

D'après ma façon de penser, vous devriez tous les deux obtenir des votes positifs, ce que je fais.

Michael R. Chernick

La suggestion de variable fictive est bonne à condition que la variance d'erreur ne diffère pas sensiblement entre les modèles. Sinon, vous pourriez appliquer un test t de Satterthwaite-Welch (qui a l'avantage singulier d'être disponible lorsque seules les statistiques sommaires sont connues, comme c'est souvent le cas lors de la lecture des articles publiés) ou d'utiliser les moindres carrés pondérés pour s'adapter au modèle combiné.

whuber