Tes apa yang dapat saya gunakan untuk membandingkan kemiringan dari dua atau lebih model regresi?

Saya ingin menguji perbedaan respons dua variabel terhadap satu prediktor. Berikut adalah contoh minimal yang dapat direproduksi.

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

Saya dapat melihat bahwa koefisien kemiringan berbeda:

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 

Saya punya tiga pertanyaan:

1. Bagaimana saya bisa menguji perbedaan antara lereng?
2. Bagaimana saya bisa menguji perbedaan antara varian residual?
3. Apa cara sederhana dan efektif untuk menyajikan perbandingan ini?

Sebuah pertanyaan terkait, Metode untuk membandingkan koefisien variabel dalam dua model regresi , menyarankan menjalankan kembali model dengan variabel dummy untuk membedakan lereng, apakah ada opsi yang akan memungkinkan penggunaan set data independen?

Bagaimana saya bisa menguji perbedaan antara lereng?

Sertakan boneka untuk spesies, biarkan berinteraksi dengan , dan lihat apakah boneka ini penting. Biarkan menjadi panjang sepal dan menjadi lebar pedal dan menjadi variabel dummy untuk ketiga spesies. Bandingkan modelnya$P_i$$L_i$$P_i$$S_1, S_2, S_3$

dengan model yang memungkinkan efek berbeda untuk setiap spesies:$P_i$

Estimator GLS adalah MLE dan model pertama adalah submodel pada model kedua, sehingga Anda dapat menggunakan uji rasio kemungkinan di sini. Kemungkinan dapat diekstraksi menggunakan logLikfungsi dan derajat kebebasan untuk tes akan menjadi karena Anda telah menghapus parameter untuk sampai pada submodel.$4$$4$

Apa cara sederhana, efektif untuk menyajikan perbandingan?

Saya pikir cara yang paling menarik adalah memplot garis regresi untuk masing-masing spesies pada sumbu yang sama, mungkin dengan bar kesalahan berdasarkan kesalahan standar. Ini akan membuat perbedaan (atau non-perbedaan) antara spesies dan hubungannya dengan sangat jelas.$P_i$

Sunting: Saya perhatikan pertanyaan lain telah ditambahkan ke badan. Jadi, saya menambahkan jawaban untuk itu:

Bagaimana saya bisa menguji perbedaan antara varian residual?

Untuk ini, Anda perlu membuat stratifikasi set data dan menyesuaikan model yang terpisah karena, model berbasis interaksi yang saya sarankan akan membatasi varian residual menjadi sama di setiap kelompok. Jika Anda cocok dengan model terpisah, kendala ini hilang. Dalam hal ini, Anda masih dapat menggunakan uji rasio kemungkinan (kemungkinan untuk model yang lebih besar sekarang dihitung dengan menjumlahkan kemungkinan dari tiga model terpisah). Model "null" tergantung pada apa yang ingin Anda bandingkan

• jika Anda hanya ingin menguji varians, sambil meninggalkan efek utama, maka model "null" harus menjadi model dengan interaksi yang saya tulis di atas. Tingkat kebebasan untuk tes tersebut adalah .$2$

• Jika Anda ingin menguji varians bersama dengan koefisien, maka model nol harus menjadi model pertama yang saya tulis di atas. Tingkat kebebasan untuk tes tersebut adalah .$6$

Untuk menjawab pertanyaan ini dengan kode R, gunakan yang berikut:
1. Bagaimana saya bisa menguji perbedaan antara lereng?
Jawab: Periksa nilai p ANOVA dari interaksi Petal.Width by Species, lalu bandingkan lereng menggunakan lsmeans :: lstrends, sebagai berikut.

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. Bagaimana saya bisa menguji perbedaan antara varian residual?
Jika saya memahami pertanyaan ini, Anda dapat membandingkan korelasi Pearson dengan transformasi Fisher, juga disebut "Fisher's r-to-z", sebagai berikut.

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. Apa cara sederhana, efektif untuk menyajikan perbandingan ini?
"Kami menggunakan regresi linier untuk membandingkan hubungan Panjang Sepal dengan Lebar Petal untuk masing-masing Spesies. Kami tidak menemukan interaksi yang signifikan dalam hubungan Panjang Sepal dengan Lebar Petal untuk I. Setosa (B = 0,9), I. Versicolor (B = 1,4), atau I. Virginica (B = 0,6); F (2, 144) = 1,6, p = 0,19. Perbandingan Fisher r-to-z menunjukkan bahwa korelasi Pearson untuk I. Setosa (r = 0,28) adalah secara signifikan lebih rendah (p = 0,02) dari I. Versicolor (r = 0,55). Demikian pula, korelasi untuk I. Virginica (r = 0,28) secara signifikan lebih lemah (p = 0,02) daripada yang diamati untuk I. Versicolor. "

Akhirnya, selalu visualisasikan hasil Anda!

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

Saya setuju dengan saran sebelumnya. Anda harus memasukkan model regresi berganda dengan variabel dummy untuk setiap set data. Ini akan memungkinkan Anda untuk menguji apakah intersepsi berbeda. Jika Anda juga ingin tahu apakah kemiringannya berbeda, maka Anda juga harus menyertakan interaksi antara boneka dan variabel yang dipertanyakan. Tidak ada masalah dengan fakta bahwa data independen. Perhatikan bahwa jika mereka baik independen dan (misalnya) spesies yang berbeda, maka Anda tidak akan dapat mengatakan apakah perbedaan yang Anda temukan adalah karena spesies yang berbeda atau set data yang berbeda, seperti yang sempurna bingung. Namun, tidak ada tes / bebaskan kartu bebas penjara yang akan membantu Anda mengatasi masalah itu tanpa mengumpulkan sampel baru dan menjalankan studi Anda lagi.