Jakiego testu mogę użyć do porównania nachyleń z dwóch lub więcej modeli regresji?

Chciałbym przetestować różnicę w odpowiedzi dwóch zmiennych na jeden predyktor. Oto minimalny odtwarzalny przykład.

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

Widzę, że współczynniki nachylenia są różne:

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

Mam trzy pytania:

Jak mogę sprawdzić różnicę między stokami?
Jak mogę sprawdzić różnicę między różnicami rezydualnymi?
Jaki jest prosty, skuteczny sposób przedstawienia tych porównań?

Powiązane pytanie, Metoda porównywania współczynnika zmiennej w dwóch modelach regresji , sugeruje ponowne uruchomienie modelu ze zmienną fikcyjną w celu zróżnicowania nachyleń, czy istnieją opcje, które pozwoliłyby na użycie niezależnych zestawów danych?

r data-visualization multivariate-analysis hypothesis-testing Abe
źródło

Jeśli chodzi o pierwsze pytanie, patrz stats.stackexchange.com/questions/55501/… .

russellpierce

Odpowiedzi:

Jak mogę sprawdzić różnicę między stokami?

Dołącz manekina do gatunku, pozwól mu wchodzić w interakcje z i sprawdź, czy ten manekin jest znaczący. Niech mieć długość sepal i mieć szerokość pedał i być zmienne obojętne dla trzech gatunków. Porównaj model $P_i$ $L_i$ $P_i$ $S_1, S_2, S_3$

mi ({L.}_{ja}) = β_{0} + β_{1} {P.}_{ja}

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

z modelem, który pozwala na wpływ dla każdego gatunku: $P_i$

mi ({L.}_{ja}) = α_{0} + α_{1} {S.}_{2)} + α_{2)} {S.}_{3)} + α_{4} {P.}_{ja} + α_{5} {P.}_{ja} {S.}_{2)} + α_{6} {P.}_{ja} {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$

Estymatorami GLS są MLE, a pierwszy model jest podmodelem na drugim, więc możesz tutaj użyć testu współczynnika wiarygodności. Prawdopodobieństwa można wyodrębnić za pomocą logLikfunkcji, a stopnie swobody dla testu wynoszą ponieważ usunięto parametry, aby dotrzeć do podmodelu. $4$ $4$

Jaki jest prosty, skuteczny sposób przedstawienia porównania?

Myślę, że najbardziej atrakcyjnym sposobem byłoby wykreślenie linii regresji dla każdego gatunku na tych samych osiach, być może z paskami błędów opartymi na standardowych błędach. To sprawiłoby, że różnica (lub brak różnicy) między gatunkiem a ich związkiem z bardzo widoczna. $P_i$

Edycja: Zauważyłem, że do ciała dodano kolejne pytanie. Dodaję odpowiedź na to:

Jak mogę sprawdzić różnicę między różnicami rezydualnymi?

W tym celu należy stratyfikować zestaw danych i dopasować osobne modele, ponieważ sugerowany przeze mnie model oparty na interakcji ograniczy resztkową wariancję, aby była taka sama w każdej grupie. Jeśli pasujesz do osobnych modeli, to ograniczenie znika. W takim przypadku nadal można użyć testu współczynnika wiarygodności (prawdopodobieństwo dla większego modelu jest teraz obliczane poprzez zsumowanie prawdopodobieństwa z trzech oddzielnych modeli). Model „zerowy” zależy od tego, z czym chcesz go porównać

jeśli chcesz tylko przetestować wariancję, pozostawiając główne efekty, model „zerowy” powinien być modelem z interakcjami, które napisałem powyżej. Stopnie swobody dla testu wynoszą wówczas . $2$
Jeśli chcesz przetestować wariancję łącznie ze współczynnikami, to model zerowy powinien być pierwszym modelem, który napisałem powyżej. Stopnie swobody dla testu wynoszą wówczas . $6$

Makro
źródło

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 Ładna odpowiedź (+1)! Zastanawiam się, czy mógłbyś dopasować glsmodel, ale pozwalając na różne odchylenia resztkowe dla każdego gatunku z opcją weights=varIdent(form=~1|Species)(w odniesieniu do drugiego pytania)?

COOLSerdash

Aby odpowiedzieć na te pytania za pomocą kodu R, wykonaj następujące czynności:
1. Jak mogę sprawdzić różnicę między zboczami?
Odpowiedź: Zbadaj wartość p ANOVA z interakcji Petal.Width według gatunków, a następnie porównaj stoki za pomocą lsmeans :: lstrends w następujący sposób.

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. Jak mogę sprawdzić różnicę między różnicami rezydualnymi?
Jeśli rozumiem pytanie, możesz porównać korelacje Pearsona z transformacją Fishera, zwaną także „r-to-z Fishera”, jak następuje.

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. Jaki jest prosty, skuteczny sposób prezentacji tych porównań?
„Zastosowaliśmy regresję liniową do porównania zależności długości oddzielenia od szerokości płatka dla każdego gatunku. Nie znaleźliśmy istotnej interakcji w zależnościach długości oddzielenia od szerokości płatka dla I. Setosa (B = 0,9), I. Versicolor (B = 1,4), ani I. Virginica (B = 0,6); F (2, 144) = 1,6, p = 0,19. Porównanie r-to-z Fishera wykazało, że korelacja Pearsona dla I. Setosa (r = 0,28) wynosiła istotnie niższa (p = 0,02) niż I. Versicolor (r = 0,55) Podobnie korelacja dla I. Virginica (r = 0,28) była znacznie słabsza (p = 0,02) niż obserwowana dla I. Versicolor. ”

Wreszcie, zawsze wizualizuj swoje wyniki!

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
źródło

Zgadzam się z poprzednią sugestią. Należy dopasować model regresji wielokrotnej ze zmienną fikcyjną dla każdego zestawu danych. Umożliwi to sprawdzenie, czy przechwyty różnią się. Jeśli chcesz również wiedzieć, czy stoki różnią się, musisz również uwzględnić interakcje między manekinami i zmienną, o której mowa. Nie ma problemu z tym, że dane są niezależne. Pamiętaj, że jeśli są one zarówno niezależne, jak i (na przykład) różne gatunki, nie byłbyś w stanie stwierdzić, czy zauważona różnica wynika z różnych gatunków, czy z różnych zestawów danych, ponieważ są one całkowicie pomieszane. Jednak nie ma testu / karta „wyjdź z więzienia”, która obejdzie ten problem bez zebrania nowej próbki i ponownego uruchomienia badania.

gung - Przywróć Monikę
źródło

Wygląda na to, że opublikowaliśmy dość podobne odpowiedzi w tym samym czasie. +1

Makro

@Macro, tak, ale twój jest w większości lepszy (+1 wcześniej); odpowiedziałeś na wszystkie 3 pytania, które przeoczyłem podczas pierwszego (nie dogłębnego) przeczytania pytania. Mój wkład tutaj dotyczy części o zamieszaniu.

gung - Przywróć Monikę

tak, to dobra uwaga. Przypuszczam, że gdybyś w ogóle przeprowadzał to zapytanie, musiałbyś działać przy założeniu, że zbiory danych mierzyły to samo itp., Z tą różnicą, że gatunki były różne.

Makro

Z mojego sposobu myślenia oboje powinniście uzyskać pozytywne głosy, co właśnie robię.

Michael R. Chernick

Sugestia zmiennej manekina jest dobra, pod warunkiem, że wariancja błędu nie różni się znacznie między modelami. W przeciwnym razie możesz zastosować test t Satterthwaite-Welch (który ma tę szczególną zaletę, że jest dostępny, gdy znane są tylko statystyki podsumowujące, jak to często ma miejsce podczas czytania opublikowanych artykułów) lub użyć ważonych najmniejszych kwadratów, aby dopasować model łączony.

whuber