Membandingkan tingkat faktor setelah GLM di R

25

Berikut adalah sedikit latar belakang tentang situasi saya: data saya merujuk pada jumlah mangsa yang berhasil dimakan oleh pemangsa. Karena jumlah mangsa terbatas (25 tersedia) di setiap percobaan, saya memiliki kolom "Sampel" yang mewakili jumlah mangsa yang tersedia (jadi, 25 di setiap percobaan), dan yang lain disebut "Hitungan" yang merupakan jumlah keberhasilan ( berapa banyak mangsa yang dimakan). Saya mendasarkan analisis saya pada contoh dari buku R pada data proporsi (halaman 578). Variabel penjelas adalah Temperatur (4 level, yang saya perlakukan sebagai faktor), dan Jenis Kelamin pemangsa (jelas, pria atau wanita). Jadi saya berakhir dengan model ini:

model <- glm(y ~ Temperature+Sex+Temperature*Sex data=predator, family=quasibinomial) 

Setelah mendapatkan tabel Analysis of Deviance, ternyata Suhu dan Jenis Kelamin (tetapi bukan interaksi) memiliki pengaruh yang signifikan terhadap konsumsi mangsa. Sekarang, masalah saya: Saya perlu tahu suhu mana yang berbeda, yaitu, saya harus membandingkan 4 suhu satu sama lain. Jika saya memiliki model linier, saya akan menggunakan fungsi TukeyHSD, tetapi karena saya menggunakan GLM saya tidak bisa. Saya telah mencari melalui paket MASS dan mencoba untuk mengatur matriks kontras tetapi untuk beberapa alasan itu tidak berhasil. Ada saran atau referensi?

Inilah ringkasan yang saya dapatkan dari model saya, jika itu membantu membuatnya lebih jelas ...

y <- cbind(data$Count, data$Sample-data$Count)
model <- glm(y ~ Temperature+Sex+Temperature*Sex data=predator, family=quasibinomial) 
> summary(model)

# Call:
# glm(formula = y ~ Temperature + Sex + Temperature * Sex, family=quasibinomial, data=data)

# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -3.7926  -1.4308  -0.3098   0.9438   3.6831  

# Coefficients:
#                                        Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                             -1.6094     0.2672  -6.024 3.86e-08 ***
# Temperature8                             0.3438     0.3594   0.957   0.3414    
# Temperature11                           -1.0296     0.4803  -2.144   0.0348 *  
# Temperature15                           -1.2669     0.5174  -2.449   0.0163 *  
# SexMale                                    0.3822     0.3577   1.069   0.2882    
# Temperature8:SexMale                    -0.2152     0.4884  -0.441   0.6606    
# Temperature11:SexMale                    0.4136     0.6093   0.679   0.4990    
# Temperature15:SexMale                    0.4370     0.6503   0.672   0.5033    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

# (Dispersion parameter for quasibinomial family taken to be 2.97372)    
#     Null deviance: 384.54  on 95  degrees of freedom
# Residual deviance: 289.45  on 88  degrees of freedom
# AIC: NA   
# Number of Fisher Scoring iterations: 5
Anne
sumber
2
Hai @Anne dan selamat datang. Anda dapat mencoba menggunakan glhtfungsi dalam multcomppaket . Untuk melakukan tes TukeyHSD untuk suhu, gunakan seperti itu glht(my.glm, mcp(Temperature="Tukey")). Dan btw: Rumus Model Anda dapat disingkat: model<-glm(y ~ Temperature*Sex data=predator, family=quasibinomial). Dengan tanda bintang ( ) interaksi dan efek utama dipasang.
COOLSerdash
Hai, terima kasih atas balasan cepat Anda! Namun saya harus melakukan sesuatu yang salah karena saya hanya mendapatkan pesan kesalahan ... Saya berasumsi bahwa my.glm adalah glm yang saya lakukan sebelumnya (karena itu, "model" dalam kasus ini). Apa yang dimaksud dengan mcp? Saya mendapatkan pesan kesalahan yang mengatakan bahwa dimensi koefisien dan matriks kovarian tidak cocok ...?
Anne
Akan sangat membantu jika Anda akan mengedit pertanyaan Anda dan memasukkan keluaran model.
COOLSerdash
3
Mengapa Anda menjadi model Temperaturesebagai faktor? Apakah Anda tidak memiliki nilai numerik yang sebenarnya? Saya akan menggunakannya sebagai variabel kontinu & kemudian seluruh masalah ini diperdebatkan.
gung - Reinstate Monica
3
Sangat masuk akal untuk ingin tahu bagaimana melakukan ini secara umum; pertanyaan Anda berdiri. Namun, dengan memperhatikan situasi spesifik Anda, saya akan menggunakan temp sebagai variabel kontinu bahkan jika Anda awalnya menganggapnya sebagai faktor. Mengesampingkan masalah dengan beberapa perbandingan, memodelkan temp sebagai faktor adalah penggunaan info yang Anda miliki tidak efisien.
gung - Reinstate Monica

Jawaban:

15

Anne, saya akan menjelaskan bagaimana melakukan beberapa perbandingan secara umum. Mengapa ini tidak berhasil dalam kasus spesifik Anda, saya tidak tahu; Maafkan saya.

Tetapi biasanya, Anda bisa melakukannya dengan multcomppaket dan fungsinya glht. Berikut ini sebuah contoh:

mydata      <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
mydata$rank <- factor(mydata$rank)
my.mod      <- glm(admit~gre+gpa*rank, data=mydata, family=quasibinomial)

summary(my.mod)
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)  
# (Intercept) -4.985768   2.498395  -1.996   0.0467 *
# gre          0.002287   0.001110   2.060   0.0400 *
# gpa          1.089088   0.731319   1.489   0.1372  
# rank2        0.503294   2.982966   0.169   0.8661  
# rank3        0.450796   3.266665   0.138   0.8903  
# rank4       -1.508472   4.202000  -0.359   0.7198  
# gpa:rank2   -0.342951   0.864575  -0.397   0.6918  
# gpa:rank3   -0.515245   0.935922  -0.551   0.5823  
# gpa:rank4   -0.009246   1.220757  -0.008   0.9940  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Jika Anda ingin menghitung perbandingan berpasangan antara rankmenggunakan Tukey's HSD, Anda bisa melakukannya dengan cara ini:

library(multcomp)
summary(glht(my.mod, mcp(rank="Tukey")))
# 
#    Simultaneous Tests for General Linear Hypotheses
# 
# Multiple Comparisons of Means: Tukey Contrasts
# 
# Fit: glm(formula = admit ~ gre + gpa * rank, family = quasibinomial, data = mydata)   
# 
# Linear Hypotheses:
#            Estimate Std. Error z value Pr(>|z|)
# 2 - 1 == 0   0.5033     2.9830   0.169    0.998
# 3 - 1 == 0   0.4508     3.2667   0.138    0.999
# 4 - 1 == 0  -1.5085     4.2020  -0.359    0.984
# 3 - 2 == 0  -0.0525     2.6880  -0.020    1.000
# 4 - 2 == 0  -2.0118     3.7540  -0.536    0.949
# 4 - 3 == 0  -1.9593     3.9972  -0.490    0.960
# (Adjusted p values reported -- single-step method)
# 
# Warning message:
# In mcp2matrix(model, linfct = linfct) :
#   covariate interactions found -- default contrast might be inappropriate

hal

Catatan: Seperti yang dicatat @gung dalam komentar, Anda harus - bila memungkinkan - memasukkan suhu sebagai variabel kontinu daripada variabel kategorikal. Mengenai interaksi: Anda dapat melakukan tes rasio kemungkinan untuk memeriksa apakah istilah interaksi secara signifikan meningkatkan kesesuaian model. Dalam kasus Anda, kode akan terlihat seperti itu:

# Original model
model <- glm(y ~ Temperature+Sex+Temperature*Sex, data=predator, family=quasibinomial) 

# Model without an interaction
model2 <- glm(y ~ Temperature+Sex data=predator, family=quasibinomial) 

# Likelihood ratio test
anova(model, model2, test="LRT")

Jika tes ini tidak signifikan, Anda dapat menghapus interaksi dari model Anda. Mungkin glhtakan berhasil?

COOLSerdash
sumber
1
Ya Tuhan, terima kasih banyak !! Saya sudah bisa menulis perintah dengan benar kali ini dan berhasil! Terima kasih lagi !
Anne
1
Pertanyaan tambahan: apakah ada cara untuk mendapatkan beberapa perbandingan pada interaksi? Saya punya data yang serupa, di mana interaksi (dari pertanyaan awal, itu akan menjadi Suhu * Jenis Kelamin) adalah signifikan, dan saya bertanya-tanya apakah mungkin untuk membandingkannya bersama-sama ...
Anne
1
Apakah maksud Anda beberapa perbandingan untuk setiap tingkat interaksi? Jika ya, Anda mungkin menemukan situs ini menarik (paragraf terakhir menunjukkan cara menguji semua kemungkinan kombinasi berpasangan).
COOLSerdash
Anda bisa membuat variabel yang sesuai dengan interaksi untuk variabel dan menggunakan variabel ini untuk menjalankan mcp. Anda melakukannya seperti ini. mydata $ gparank <- interaksi (mydata $ gpa, mydata $ rank)
Notquitesure
1
@Nova tautan mana yang Anda maksud? Yang ada di komentar? Ini tautan baru ke situs itu.
COOLSerdash