Menyesuaikan model linear umum heteroskedastik untuk respons binomial

9

Saya memiliki data dari desain eksperimental berikut: pengamatan saya adalah jumlah dari jumlah keberhasilan ( K) dari jumlah percobaan yang sesuai ( N), diukur untuk dua kelompok yang masing-masing terdiri dari Iindividu, dari Tperawatan, di mana dalam setiap kombinasi faktor tersebut terdapat Rulangan . Oleh karena itu, secara keseluruhan saya memiliki 2 * I * T * RK dan N yang sesuai .

Data berasal dari biologi. Setiap individu adalah gen yang saya ukur level ekspresi dari dua bentuk alternatif (karena sebuah fenomena yang disebut splicing alternatif). Karenanya, K adalah level ekspresi dari salah satu bentuk dan N adalah jumlah level ekspresi dari kedua bentuk. Pilihan antara dua bentuk dalam satu salinan yang diekspresikan diasumsikan sebagai eksperimen Bernoulli, karenanya K keluar dari Nsalinan mengikuti binomial. Setiap kelompok terdiri dari ~ 20 gen yang berbeda dan gen dalam setiap kelompok memiliki beberapa fungsi yang sama, yang berbeda antara kedua kelompok. Untuk setiap gen dalam setiap kelompok saya memiliki ~ 30 pengukuran seperti itu dari masing-masing tiga jaringan yang berbeda (perawatan). Saya ingin memperkirakan efek yang dimiliki kelompok dan pengobatan terhadap varian K / N.

Ekspresi gen diketahui overdispersi sehingga penggunaan binomial negatif dalam kode di bawah ini.

Misalnya, Rkode data yang disimulasikan:

library(MASS)
set.seed(1)
I = 20 # individuals in each group
G = 2  # groups
T = 3  # treatments
R = 30 # replicates of each individual, in each group, in each treatment

groups     = letters[1:G]
ids        = c(sapply(groups, function(g){ paste(rep(g, I), 1:I, sep=".") }))
treatments = paste(rep("t", T), 1:T, sep=".")
 # create random mean number of trials for each individual and 
 #  dispersion values to simulate trials from a negative binomial:
mean.trials = rlnorm(length(ids), meanlog=10, sdlog=1)
thetas      = 10^6/mean.trials
 # create the underlying success probability for each individual:
p.vec = runif(length(ids), min=0, max=1)
 # create a dispersion factor for each success probability, where the 
 #  individuals of group 2 have higher dispersion thus creating a group effect:
dispersion.vec = c(runif(length(ids)/2, min=0, max=0.1),
                   runif(length(ids)/2, min=0, max=0.2))
 # create empty an data.frame:
data.df = data.frame(id=rep(sapply(ids, function(i){ rep(i, R) }), T),
                     group=rep(sapply(groups, function(g){ rep(g, I*R) }), T),
                     treatment=c(sapply(treatments, 
                                        function(t){ rep(t, length(ids)*R) })),
                     N=rep(NA, length(ids)*T*R), 
                     K=rep(NA, length(ids)*T*R) )
 # fill N's and K's - trials and successes
for(i in 1:length(ids)){
  N     = rnegbin(T*R, mu=mean.trials[i], theta=thetas[i])
  probs = runif(T*R, min=max((1-dispersion.vec[i])*p.vec[i],0),
                max=min((1+dispersion.vec)*p.vec[i],1))
  K     = rbinom(T*R, N, probs)
  data.df$N[which(as.character(data.df$id) == ids[i])] = N
  data.df$K[which(as.character(data.df$id) == ids[i])] = K
}

Saya tertarik untuk memperkirakan efek yang dimiliki kelompok dan pengobatan terhadap dispersi (atau varians) dari probabilitas keberhasilan (yaitu, K/N). Oleh karena itu saya mencari glm yang sesuai di mana responsnya adalah K / N tetapi selain memodelkan nilai respons yang diharapkan, varian respons juga dimodelkan.

Jelas, varians dari probabilitas keberhasilan binomial dipengaruhi oleh jumlah percobaan dan probabilitas keberhasilan yang mendasarinya (semakin tinggi jumlah percobaan dan semakin ekstrim probabilitas keberhasilan yang mendasarinya adalah (yaitu, dekat 0 atau 1), semakin rendah varians dari probabilitas keberhasilan), jadi saya terutama tertarik pada kontribusi kelompok dan perawatan di luar jumlah percobaan dan probabilitas keberhasilan yang mendasarinya. Saya kira menerapkan transformasi akar arcsin kuadrat untuk respon akan menghilangkan yang terakhir tetapi tidak dari jumlah percobaan.

Meskipun dalam contoh data simulasi di atas desainnya seimbang (jumlah individu yang sama dalam setiap dua kelompok dan jumlah ulangan yang sama pada setiap individu dari setiap kelompok dalam setiap perlakuan), dalam data saya yang sebenarnya tidak - kedua kelompok melakukan tidak memiliki jumlah individu yang sama dan jumlah ulangan bervariasi. Juga, saya membayangkan individu harus ditetapkan sebagai efek acak.

Merencanakan varians sampel vs. rata-rata sampel dari estimasi probabilitas keberhasilan (dilambangkan sebagai p = K / N) dari masing-masing individu menggambarkan bahwa probabilitas keberhasilan ekstrem memiliki varians yang lebih rendah: masukkan deskripsi gambar di sini masukkan deskripsi gambar di sini

Ini dihilangkan ketika probabilitas keberhasilan yang diperkirakan ditransformasikan menggunakan transformasi penstabilan varians akar arcsin kuadrat (dilambangkan sebagai arcsin (sqrt (p hat)): masukkan deskripsi gambar di sini

Merencanakan varians sampel dari probabilitas keberhasilan estimasi yang ditransformasikan vs rata-rata N menunjukkan hubungan negatif yang diharapkan: masukkan deskripsi gambar di sini

Merencanakan varians sampel dari probabilitas keberhasilan estimasi yang ditransformasikan untuk kedua kelompok menunjukkan bahwa kelompok b memiliki varians yang sedikit lebih tinggi, yang merupakan cara saya mensimulasikan data: masukkan deskripsi gambar di sini

Akhirnya, merencanakan varians sampel dari probabilitas keberhasilan estimasi yang ditransformasikan untuk tiga perawatan tidak menunjukkan perbedaan antara perawatan, yang merupakan cara saya mensimulasikan data: masukkan deskripsi gambar di sini

Apakah ada bentuk model linier umum yang dengannya saya dapat mengukur efek kelompok dan pengobatan pada varians probabilitas keberhasilan?

Mungkin model linear umum heteroskedastik atau beberapa bentuk model varians loglinear?

Sesuatu dalam garis model yang memodelkan Variance (y) = Zλ di samping E (y) = Xβ, di mana Z dan X adalah masing-masing regresi dari mean dan varians, yang dalam kasus saya akan identik dan termasuk pengobatan (level t.1, t.2, dan t.3) dan kelompok (level a dan b), dan mungkin N dan R, dan karenanya λ dan β akan memperkirakan efek masing-masing.

Sebagai alternatif, saya dapat menyesuaikan model dengan varians sampel di seluruh ulangan masing-masing gen dari masing-masing kelompok dalam setiap perlakuan, menggunakan glm yang hanya memodelkan nilai respons yang diharapkan. Satu-satunya pertanyaan di sini adalah bagaimana menjelaskan fakta bahwa gen yang berbeda memiliki jumlah ulangan yang berbeda. Saya pikir bobot dalam glm dapat menjelaskan hal itu (varian sampel yang didasarkan pada lebih banyak ulangan harus memiliki bobot lebih tinggi) tetapi sebenarnya bobot mana yang harus ditetapkan?

Catatan: Saya sudah mencoba menggunakan dglmpaket R:

library(dglm)
dglm.fit = dglm(formula = K/N ~ 1, dformula = ~ group + treatment, family = quasibinomial, weights = N, data = data.df)
summary(dglm.fit)
Call: dglm(formula = K/N ~ 1, dformula = ~group + treatment, family = quasibinomial, 
    data = data.df, weights = N)

Mean Coefficients:
               Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -0.09735366 0.01648905 -5.904138 3.873478e-09
(Dispersion Parameters for quasibinomial family estimated as below )

    Scaled Null Deviance: 3600 on 3599 degrees of freedom
Scaled Residual Deviance: 3600 on 3599 degrees of freedom

Dispersion Coefficients:
                Estimate Std. Error      z value  Pr(>|z|)
(Intercept)  9.140517930 0.04409586 207.28746254 0.0000000
group       -0.071009599 0.04714045  -1.50634107 0.1319796
treatment   -0.001469108 0.02886751  -0.05089138 0.9594121
(Dispersion parameter for Gamma family taken to be 2 )

    Scaled Null Deviance: 3561.3 on 3599 degrees of freedom
Scaled Residual Deviance: 3559.028 on 3597 degrees of freedom

Minus Twice the Log-Likelihood: 29.44568 
Number of Alternating Iterations: 5 

Efek grup menurut dglm.fit cukup lemah. Saya bertanya-tanya apakah modelnya diatur dengan benar atau kekuatan yang dimiliki model ini.

pengguna1701545
sumber
1
Ini agak buram. Bisakah Anda mengatakan secara lebih spesifik apa situasi Anda, mungkin secara konkret?
gung - Reinstate Monica
1
Saya mengedit posting saya. Tapi secara umum, saya mencari glm yang selain E (y) = Xβ, model Variance (y) = Zλ, di mana kolom X adalah regresi untuk rata-rata respons dan kolom Z adalah regresi untuk varians dari respon, dan β dan λ adalah parameter yang memperkirakan efeknya pada mean dan varians, masing-masing.
user1701545
2
Saya mengumpulkan sebanyak itu, tetapi sulit untuk memikirkan hal ini dengan cara yang samar-samar. Anda tidak mendiskusikan variabel X atau Z, mis. Tentang apa studi Anda? Apa data Anda? Dll
gung - Reinstate Monica
Tidak masalah - lihat posting terbaru.
user1701545

Jawaban:

1

Mungkin yang Anda cari adalah sesuatu yang disebut model linier umum ganda di mana parameter mean dan dispersi dimodelkan. Bahkan ada paket R dglm yang dirancang agar sesuai dengan model tersebut.

pengguna3897081
sumber