Mit welchem Test kann ich Steigungen aus zwei oder mehr Regressionsmodellen vergleichen?

Ich möchte den Unterschied in der Reaktion zweier Variablen auf einen Prädiktor testen. Hier ist ein minimal reproduzierbares Beispiel.

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")

Ich kann sehen, dass die Steigungskoeffizienten unterschiedlich sind:

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

Ich habe drei Fragen:

Wie kann ich den Unterschied zwischen Pisten testen?
Wie kann ich den Unterschied zwischen Restabweichungen testen?
Wie lassen sich diese Vergleiche auf einfache und effektive Weise darstellen?

Eine verwandte Frage, Methode zum Vergleichen des Variablenkoeffizienten in zwei Regressionsmodellen , schlägt vor, das Modell mit einer Dummy-Variablen erneut auszuführen, um die Steigungen zu unterscheiden. Gibt es Optionen, die die Verwendung unabhängiger Datensätze ermöglichen würden?

r data-visualization multivariate-analysis hypothesis-testing Abe
quelle

Bezüglich der ersten Frage siehe stats.stackexchange.com/questions/55501/… .

Russellpierce

Antworten:

Wie kann ich den Unterschied zwischen Pisten testen?

Fügen Sie einen Dummy für Arten hinzu, lassen Sie ihn mit interagieren und , ob dieser Dummy von Bedeutung ist. Sei die Kelchblattlänge und die Pedalbreite und die Dummy-Variablen für die drei Arten. Die vergleichen das Modell $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$

mit dem Modell, mit dem die Wirkung von für jede Spezies unterschiedlich sein kann: $P_i$

E (L_{ich}) = α_{0} + α_{1} S_{2} + α_{2} S_{3} + α_{4} P_{ich} + α_{5} P_{ich} S_{2} + α_{6} P_{ich} 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$

logLik $4$ $4$

Was ist ein einfacher und effektiver Weg, um den Vergleich zu präsentieren?

$P_i$

Bearbeiten: Mir ist aufgefallen, dass dem Text eine weitere Frage hinzugefügt wurde. Also füge ich eine Antwort hinzu:

Wie kann ich den Unterschied zwischen Restabweichungen testen?

Dazu müssen Sie den Datensatz schichten und separate Modelle anpassen, da das von mir vorgeschlagene interaktionsbasierte Modell die Restvarianz so einschränkt, dass sie in jeder Gruppe gleich ist. Wenn Sie separate Modelle anpassen, wird diese Einschränkung aufgehoben. In diesem Fall können Sie den Likelihood-Ratio-Test weiterhin verwenden (die Wahrscheinlichkeit für das größere Modell wird jetzt berechnet, indem die Wahrscheinlichkeiten aus den drei separaten Modellen summiert werden). Das "Null" -Modell hängt davon ab, womit Sie es vergleichen möchten

Makro
quelle

S_{1}

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

S_{1}

$S_1$

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

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

@Macro Nette Antwort (+1)! Ich frage mich, ob Sie das glsModell passen könnten, aber mit der Option weights=varIdent(form=~1|Species)(in Bezug auf die zweite Frage) unterschiedliche Restabweichungen für jede Art berücksichtigen könnten .

COOLSerdash

Um diese Fragen mit R-Code zu beantworten, verwenden Sie Folgendes:
1. Wie kann ich den Unterschied zwischen Steigungen testen?
Antwort: Untersuchen Sie den ANOVA-p-Wert aus der Interaktion von Petal.Width nach Species und vergleichen Sie dann die Steigungen mit lsmeans :: lstrends wie folgt.

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. Wie kann ich die Differenz zwischen Restabweichungen testen?
Wenn ich die Frage verstehe, können Sie Pearson-Korrelationen wie folgt mit einer Fisher-Transformation vergleichen, die auch als "Fisher's r-to-z" bezeichnet wird.

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. Wie lassen sich diese Vergleiche auf einfache und effektive Weise darstellen?
Wir verwendeten eine lineare Regression, um das Verhältnis von Kelchblattlänge zu Blütenblattbreite für jede Art zu vergleichen. Wir fanden keine signifikante Wechselwirkung in den Beziehungen von Kelchblattlänge zu Blütenblattbreite für I. Setosa (B = 0,9), I. Versicolor (B = 1,4), noch I. Virginica (B = 0,6), F (2, 144) = 1,6, p = 0,19 Ein Fisher-r-zu-z-Vergleich zeigte, dass die Pearson-Korrelation für I. Setosa (r = 0,28) war signifikant niedriger (p = 0,02) als I. Versicolor (r = 0,55) Ebenso war die Korrelation für I. Virginica (r = 0,28) signifikant schwächer (p = 0,02) als die für I. Versicolor beobachtete. "

Schließlich immer visualisieren Sie Ihre Ergebnisse!

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
quelle

Ich stimme dem vorherigen Vorschlag zu. Sie sollten ein Modell mit mehreren Regressionen mit einer Dummy-Variablen für jeden Datensatz anpassen. Auf diese Weise können Sie testen, ob sich die Abschnitte unterscheiden. Wenn Sie auch wissen möchten, ob sich die Steigungen unterscheiden, müssen Sie auch Interaktionen zwischen den Dummies und der betreffenden Variablen einbeziehen. Es ist kein Problem damit, dass die Daten unabhängig sind. Beachten Sie, dass , wenn sie sowohl unabhängig und (zum Beispiel) verschiedene Arten, dann würden Sie nicht in der Lage sein zu sagen , ob die Differenz Sie finden aufgrund der unterschiedlichen Spezies ist oder die unterschiedlichen Datensätzen, wie sie perfekt verwechselt werden. Es gibt jedoch keinen Test / Get-out-of-Jail-Free-Karte, mit der Sie dieses Problem umgehen können, ohne eine neue Probe zu sammeln und Ihre Studie erneut durchzuführen.

gung - Wiedereinsetzung von Monica
quelle

Es sieht so aus, als hätten wir ziemlich ähnliche Antworten fast zur gleichen Zeit gepostet. +1

Makro

@Macro, ja, aber deins ist meistens besser (+1 früher); Sie haben alle drei Fragen angesprochen, die mir bei meiner ersten (nicht gründlichen) Lektüre der Frage entgangen sind. Mein Beitrag hier ist der Teil über das Verwirren.

gung - Wiedereinsetzung von Monica

ja das ist ein guter punkt. Ich nehme an, wenn Sie diese Untersuchung überhaupt durchführen würden, müssten Sie unter der Annahme arbeiten, dass die Datensätze dasselbe messen, etc ... mit dem einzigen Unterschied, dass die Arten unterschiedlich sind.

Makro

Aus meiner Sicht solltet ihr beide Upvotes bekommen, was ich tue.

Michael R. Chernick

Der Dummy-Variablenvorschlag ist gut, vorausgesetzt, die Fehlervarianz unterscheidet sich nicht nennenswert zwischen den Modellen. Andernfalls können Sie einen Satterthwaite-Welch-T-Test anwenden (der den einzigartigen Vorteil hat, verfügbar zu sein, wenn nur zusammenfassende Statistiken bekannt sind, wie dies häufig beim Lesen von veröffentlichten Artikeln der Fall ist) oder gewichtete kleinste Quadrate verwenden, um das kombinierte Modell anzupassen.

whuber

Mit welchem ​​Test kann ich Steigungen aus zwei oder mehr Regressionsmodellen vergleichen?

Antworten:

Mit welchem Test kann ich Steigungen aus zwei oder mehr Regressionsmodellen vergleichen?