ANCOVA dalam R menunjukkan penyadapan yang berbeda, tetapi 95% CI tumpang tindih ... bagaimana ini mungkin?

Kami memiliki kumpulan data dengan dua kovariat dan variabel pengelompokan kategoris dan ingin tahu apakah ada perbedaan yang signifikan antara kemiringan atau intersep di antara kovariat yang terkait dengan variabel pengelompokan yang berbeda. Kami telah menggunakan anova () dan lm () untuk membandingkan kecocokan dari tiga model yang berbeda: 1) dengan kemiringan tunggal dan intersep, 2) dengan intersep yang berbeda untuk setiap grup, dan 3) dengan kemiringan dan intersep untuk setiap kelompok . Menurut anova () tes linier umum, model kedua adalah yang paling tepat dari ketiganya, ada peningkatan yang signifikan untuk model dengan memasukkan intersep terpisah untuk masing-masing kelompok. Namun, ketika kita melihat interval kepercayaan 95% untuk intersep ini - semuanya tumpang tindih, menunjukkan tidak ada perbedaan yang signifikan antara intersep. Bagaimana kedua hasil ini dapat direkonsiliasi? Kami pikir cara lain untuk menginterpretasikan hasil dari metode pemilihan model adalah bahwa harus ada setidaknya satu perbedaan signifikan di antara penyadapan ... tapi mungkin ini tidak benar?

Di bawah ini adalah kode R untuk mereplikasi analisis ini. Kami telah menggunakan fungsi dput () sehingga Anda dapat bekerja dengan data yang persis sama dengan yang kami garap.

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

Singkatnya - masalahnya adalah bahwa 95% CI untuk intersep semuanya tumpang tindih, tetapi metode pemilihan model menunjukkan bahwa model terbaik adalah yang cocok dengan intersep yang berbeda. Jadi saya cenderung berpikir bahwa metode pemilihan model kami cacat atau 95% CI untuk perkiraan intersep dihitung secara tidak benar. Pikiran apa pun akan sangat dihargai!

Ingatlah bahwa perbedaan antara signifikan dan non-signifikan tidak (selalu) signifikan secara statistik

Sekarang, lebih ke titik pertanyaan Anda, model 1 disebut regresi gabungan, dan model 2 regresi tidak dikumpulkan. Seperti yang Anda catat, dalam regresi gabungan, Anda menganggap bahwa grup tidak relevan, yang berarti bahwa varians antara grup diatur ke nol.

Dalam regresi yang tidak disatukan, dengan intersep per grup, Anda mengatur varians menjadi tak terbatas.

Secara umum, saya lebih menyukai solusi perantara, yang merupakan model hierarkis atau regresi gabungan parsial (atau estimator penyusutan). Anda dapat memuat model ini dalam R dengan paket lmer4.

Akhirnya, lihatlah makalah ini oleh Gelman , di mana ia berargumen mengapa model hierarkis membantu dengan berbagai masalah perbandingan (dalam kasus Anda, apakah koefisien per kelompok berbeda? Bagaimana kita mengoreksi nilai-p untuk beberapa perbandingan).

Misalnya, dalam kasus Anda,

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

Jika Anda ingin menyesuaikan kemiringan bervariasi, kemiringan bervariasi (model ketiga), jalankan saja

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

Kemudian Anda dapat melihat varians grup dan melihat apakah itu berbeda dari nol (regresi yang dikumpulkan bukan model yang lebih baik) dan jauh dari infinity (regresi yang tidak dikumpulkan).

pembaruan: Setelah komentar (lihat di bawah), saya memutuskan untuk memperluas jawaban saya.

Tujuan dari model hierarkis, khususnya dalam kasus seperti ini, adalah untuk memodelkan variasi berdasarkan kelompok (dalam hal ini, Situs). Jadi, daripada menjalankan ANOVA untuk menguji apakah suatu model berbeda dari yang lain, saya akan melihat prediksi model saya dan melihat apakah prediksi berdasarkan kelompok lebih baik dalam model hierarkis vs regresi yang dikumpulkan (regresi klasik) .

Sekarang, saya menjalankan sugestions saya di atas dan foudn itu

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

Akan mengembalikan nol sebagai taksiran kemiringan yang berbeda-beda (efek yang bervariasi dari kepala ke situs). lalu aku berlari

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

Dan saya mendapat perkiraan nihil untuk berbagai efek kepala. Saya belum tahu mengapa ini terjadi, karena ini adalah pertama kalinya saya menemukan ini. Saya benar-benar minta maaf untuk masalah ini, tetapi, dalam pembelaan saya, saya hanya mengikuti spesifikasi yang diuraikan dalam bantuan fungsi lmer. (Lihat contoh dengan sleepstudy data). Saya akan mencoba memahami apa yang terjadi dan saya akan melaporkan di sini ketika (jika) saya mengerti apa yang terjadi.

Sebelum ada intervensi moderator, Anda bisa melihatnya

library(car)

crPlots(lm.2,terms=~Site)

Ini adalah Plot Komponen + Sisa (Sisa Sebagian)

Saya pikir, antara lain, bahwa Anda salah menghitung interval kepercayaan. Berikut adalah dua cara untuk melihatnya:

(1) perbedaan dari setiap situs dari situs baseline (ANZ) [Anda juga bisa menghitung perbedaan dari rata-rata keseluruhan dengan mengubah ke jumlah-ke-nol-kontras

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

atau (2) semua perbandingan berpasangan (Saya tidak menyukai pendekatan ini, tapi ini biasa):

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

Perhatikan bahwa semua Headnilai Anda berada dalam kisaran 1,7 - 2,4, sementara intersep mencoba memperkirakan Legnilai pada Head=0. Ini adalah ekstrapolasi utama, jadi ada banyak ketidakpastian. Jika Anda memusatkan Headnilai - nilai, dan mengulangi analisis ini, interval kepercayaan akan semakin ketat.

Selain itu, interval kepercayaan 95% yang tumpang tindih tidak menyiratkan kurangnya perbedaan yang signifikan secara statistik. Bahkan, untuk dua kelompok, interval kepercayaan 84% yang tidak tumpang tindih mendekati yang memiliki perbedaan signifikan pada tingkat 5%. Tentu saja, karena banyak pengujian, ini tidak cukup untuk beberapa grup.

Selain jawaban lain, berikut adalah beberapa tautan dari Unit Konsultasi Statistik Cornell yang relevan dengan interval kepercayaan yang tumpang tindih dan berfungsi sebagai pengingat yang baik dan singkat tentang apa yang mereka lakukan dan tidak maksudkan.

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

Inilah poin utamanya:

Jika dua statistik memiliki interval kepercayaan yang tidak tumpang tindih, mereka tentu berbeda secara signifikan tetapi jika mereka memiliki interval kepercayaan yang tumpang tindih, itu tidak selalu benar bahwa mereka tidak berbeda secara signifikan.

Inilah teks yang relevan dari tautan pertama:

Kita dapat menggambarkan ini dengan contoh sederhana. Misalkan kita tertarik untuk membandingkan cara dari dua sampel independen. Mean dari sampel pertama adalah 9 dan rata-rata dari sampel kedua adalah 17. Mari kita asumsikan bahwa dua kelompok berarti memiliki kesalahan standar yang sama, sama dengan 2,5. Interval kepercayaan 95 persen untuk rata-rata kelompok pertama dapat dihitung sebagai: ± × 5,296.19 di mana 1,96 adalah nilai t kritis. Interval kepercayaan untuk rata-rata kelompok pertama demikian (4.1, 13.9). Demikian pula untuk kelompok kedua, interval kepercayaan untuk rata-rata adalah (12.1, 21.9). Perhatikan bahwa kedua interval tumpang tindih. Namun, t-statistik untuk membandingkan dua cara adalah:

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

yang mencerminkan bahwa hipotesis nol, bahwa rata-rata dari kedua kelompok adalah sama, harus ditolak pada tingkat α = 0,05. Untuk memverifikasi kesimpulan di atas, pertimbangkan interval kepercayaan 95 persen untuk perbedaan antara rata-rata dua kelompok: (17-9) ± 1,96 x √ (2,5² + 2,5²) yang menghasilkan (1,09, 14,91). Interval tidak mengandung nol, maka kami menolak hipotesis nol bahwa mean kelompok adalah sama.

Secara umum, ketika membandingkan dua estimasi parameter, selalu benar bahwa jika interval kepercayaan tidak tumpang tindih, maka statistik akan berbeda secara statistik. Namun, kebalikannya tidak benar. Artinya, adalah salah untuk menentukan signifikansi statistik dari perbedaan antara dua statistik berdasarkan interval kepercayaan yang tumpang tindih. Untuk penjelasan mengapa ini berlaku untuk kasus perbandingan dua sampel cara, lihat tautan berikut: http://www.cscu.cornell.edu/news/statnews/Stnews73insert.pdf

Ini info dari tautan lain: