Glmer binomial () tunggal meskipun "banyak" data
Saya telah membaca berbagai pertanyaan tentang betapa mudahnya menemukan singularitas saat menyesuaikan model efek campuran glmer(). Secara umum, idenya adalah bahwa singularitas mungkin muncul dari struktur acak yang sangat kompleks. Jika struktur acaknya sederhana, mungkin juga terjadi ketika datanya tidak cukup untuk menghitung matriks varians-kovarian ... lihat misalnya halaman ini oleh Ben Bolker , jawaban Robert Long untuk posting ini atau halaman bantuanisSingular() .
Namun, model yang saya coba paskan sangat sederhana:
mod.detection_rand <- glmer(reaction ~ Pedra + (1|Channel), family="binomial", data = garotes)
boundary (singular) fit: see ?isSingular
... dan ternyata saya memiliki cukup data untuk berbagai kombinasi variabel prediktor (tetap dan acak):
library(tidyverse)
garotes %>%
group_by(Channel, Pedra) %>%
summarise(n = n())
# A tibble: 16 x 3
# Groups: Channel [8]
Channel Pedra n
<int> <fct> <int>
1 1 No 13
2 1 Yes 13
3 2 No 14
4 2 Yes 12
5 3 No 12
6 3 Yes 14
7 4 No 13
8 4 Yes 13
9 5 No 13
10 5 Yes 13
11 6 No 14
12 6 Yes 12
13 7 No 13
14 7 Yes 13
15 8 No 14
16 8 Yes 12
Bagaimana menurut anda?
EDIT: Berikut ringkasan modelnya, summary(mod.detection_rand)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: reaction ~ Pedra + (1 | Channel)
Data: garotes
AIC BIC logLik deviance df.resid
261.5 271.5 -127.7 255.5 205
Scaled residuals:
Min 1Q Median 3Q Max
-1.8533 -0.9449 0.5396 0.5396 1.0583
Random effects:
Groups Name Variance Std.Dev.
Channel (Intercept) 0 0
Number of obs: 208, groups: Channel, 8
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.1133 0.1946 -0.582 0.56
PedraYes 1.3473 0.3066 4.394 1.11e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
PedraYes -0.635
convergence code: 0
boundary (singular) fit: see ?isSingular
EDIT2: Mengikuti komentar Billy:
bobyqa : boundary (singular) fit: see ?isSingular
[OK]
Nelder_Mead : boundary (singular) fit: see ?isSingular
[OK]
nlminbwrap : boundary (singular) fit: see ?isSingular
[OK]
nmkbw : boundary (singular) fit: see ?isSingular
[OK]
optimx.L-BFGS-B : boundary (singular) fit: see ?isSingular
[OK]
nloptwrap.NLOPT_LN_NELDERMEAD : boundary (singular) fit: see ?isSingular
[OK]
nloptwrap.NLOPT_LN_BOBYQA : boundary (singular) fit: see ?isSingular
[OK]
EDIT3: Mengikuti jawaban Isabella:
Saya memeriksa struktur variabel hasil ( reaction). Berikut tabel hasilnya:
library(tidyverse)
garotes %>%
group_by(Channel, Pedra, reaction) %>%
summarise(n = n()) %>%
print(n = Inf)
# A tibble: 32 x 4
# Groups: Channel, Pedra [16]
Channel Pedra reaction n
<int> <fct> <int> <int>
1 1 No 0 6
2 1 No 1 7
3 1 Yes 0 3
4 1 Yes 1 10
5 2 No 0 7
6 2 No 1 7
7 2 Yes 0 2
8 2 Yes 1 10
9 3 No 0 8
10 3 No 1 4
11 3 Yes 0 6
12 3 Yes 1 8
13 4 No 0 7
14 4 No 1 6
15 4 Yes 0 3
16 4 Yes 1 10
17 5 No 0 8
18 5 No 1 5
19 5 Yes 0 1
20 5 Yes 1 12
21 6 No 0 6
22 6 No 1 8
23 6 Yes 0 2
24 6 Yes 1 10
25 7 No 0 6
26 7 No 1 7
27 7 Yes 0 2
28 7 Yes 1 11
29 8 No 0 8
30 8 No 1 6
31 8 Yes 0 4
32 8 Yes 1 8
Rupanya, ada kedua jenis hasil untuk semua Channelsdan semua Pedraperawatan ... jadi tidak seperti contoh yang disajikan Isabella ... selanjutnya, saya mencoba memodelkan GLMM ini dengan library(GLMMadaptive)dan tidak menyatu juga.
EDIT4: Kumpulan data yang saya gunakan, jika ada yang penasaran.
Channel Pedra reaction
1 No 1
2 No 0
3 No 0
4 No 0
5 No 0
6 No 1
7 No 0
8 No 0
1 No 1
2 No 1
3 No 1
4 No 1
5 No 0
6 No 0
7 No 0
8 No 0
1 No 0
2 No 1
3 No 0
4 No 0
5 No 0
6 No 0
7 No 0
8 No 1
1 No 0
2 No 1
3 Yes 0
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 0
1 Yes 1
2 Yes 1
3 Yes 0
4 Yes 0
5 No 0
6 No 1
7 Yes 1
8 Yes 1
1 Yes 0
2 Yes 1
3 Yes 1
4 Yes 1
5 Yes 1
6 Yes 0
7 No 1
8 No 1
1 Yes 1
2 Yes 1
3 Yes 1
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 1
1 Yes 1
2 Yes 1
3 Yes 1
4 Yes 1
5 Yes 0
6 Yes 1
7 Yes 1
8 Yes 1
1 Yes 1
2 Yes 1
3 Yes 0
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 0
8 Yes 0
1 Yes 1
2 Yes 1
3 Yes 0
4 Yes 0
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 0
1 Yes 1
2 Yes 1
3 Yes 0
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 0
8 Yes 0
1 Yes 1
2 Yes 0
3 Yes 1
4 Yes 0
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 1
1 Yes 1
2 Yes 1
3 Yes 0
4 Yes 1
5 Yes 1
6 Yes 0
7 Yes 1
8 Yes 1
1 Yes 1
2 Yes 1
3 Yes 1
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 1
1 Yes 0
2 Yes 0
3 Yes 1
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 1
1 Yes 1
2 No 0
3 Yes 1
4 No 1
5 Yes 1
6 No 1
7 Yes 1
8 No 1
1 No 0
2 Yes 1
3 No 0
4 Yes 1
5 No 1
6 Yes 1
7 No 1
8 Yes 1
1 Yes 0
2 No 1
3 Yes 1
4 No 0
5 Yes 1
6 No 1
7 Yes 1
8 No 0
1 No 0
2 No 1
3 No 1
4 No 0
5 No 1
6 No 0
7 No 0
8 No 0
1 No 1
5 No 0
3 No 1
4 No 1
2 No 1
6 No 0
7 No 1
8 No 0
1 No 0
5 No 0
3 No 0
4 No 0
2 No 1
6 No 0
7 No 0
8 No 0
1 No 1
5 No 1
3 No 1
4 No 0
2 No 0
6 No 1
7 No 1
8 No 0
1 No 1
5 No 0
3 No 0
4 No 1
2 No 0
6 No 1
7 No 1
8 No 1
1 No 1
5 No 1
3 No 0
4 No 1
2 No 0
6 No 1
7 No 1
8 No 1
1 No 1
5 No 1
3 No 0
4 No 0
2 No 0
6 No 1
7 No 0
8 No 0
1 No 0
5 No 0
3 No 0
4 No 1
2 No 0
6 No 0
7 No 1
8 No 1
Terima kasih banyak atas semua tanggapan Anda! Belajar banyak dari mereka!
Jawaban
Komentar lebih lanjut: Saya telah memeriksa data Anda, dan jelas sekali lagi, bahwa tidak ada bukti perbedaan sistematis antara saluran yang berbeda. Inilah sebabnya mengapa model campuran memperkirakan varians antara saluran menjadi$0$, membuat modelnya tunggal.
Anda dapat melihat ini pada gambar di bawah, di mana kesalahan standar untuk hampir setiap saluran tumpang tindih ...
... dan dapat memastikannya dengan dekomposisi ANOVA dari GLM efek tetap, yang menunjukkan bahwa tidak ada efek utama yang signifikan dari Channel (p = .986).
m_fixed_effects = glm(cbind(n, total) ~ Pedra + factor(Channel),
data=positive, family=binomial)
car::Anova(m_fixed_effects)
# Analysis of Deviance Table (Type II tests)
#
# Response: cbind(n, total)
# LR Chisq Df Pr(>Chisq)
# Pedra 4.9148 1 0.02663 *
# factor(Channel) 1.3859 7 0.98600
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Kode
library(tidyverse)
df = read.csv('/path/to/reaction.csv')
head(df)
# Channel Pedra reaction n
# 1 1 No 0 6
# 2 1 No 1 7
# 3 1 Yes 0 3
# 4 1 Yes 1 10
# 5 2 No 0 7
# 6 2 No 1 7
df = df %>%
group_by(Channel, Pedra) %>%
mutate(total = sum(n),
prop = n / total,
se = sqrt((prop * (1-prop)) / n)) %>%
ungroup()
positive = filter(df, reaction==1)
ggplot(positive, aes(Pedra, prop, group=Channel, color=factor(Channel))) +
geom_path(position = position_dodge(width=.1)) +
geom_point(position = position_dodge(width=.1)) +
stat_summary(fun.data=mean_se, group=1, color='black',
position = position_nudge(x=c(-.2, .2))) +
geom_linerange(mapping=aes(ymin=prop-se, ymax=prop+se),
position = position_dodge(width=.1)) +
geom_hline(linetype='dashed', yintercept=.5) +
coord_cartesian(ylim=c(0, 1)) +
labs(color='Channel', y='Proportion positive reactions',
caption='Error bars show SEM')
m_fixed_effects = glm(cbind(n, total) ~ Pedra + factor(Channel),
data=positive, family=binomial)
car::Anova(m_fixed_effects)
# Analysis of Deviance Table (Type II tests)
#
# Response: cbind(n, total)
# LR Chisq Df Pr(>Chisq)
# Pedra 4.9148 1 0.02663 *
# factor(Channel) 1.3859 7 0.98600
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Karena ini adalah model regresi logistik biner efek campuran, model ini mengasumsikan bahwa variabel hasil Anda adalah biner dengan nilai yang dikodekan sebagai 0 atau 1.
Apa yang perlu Anda selidiki adalah apakah Anda memiliki cukup 1 hadir dalam variabel respons Anda untuk jumlah 'subjek' yang cukup. (Dalam kasus Anda, subjek adalah singkatan dari saluran.)
Berikut adalah contoh buatan yang menghasilkan peringatan yang sama seperti yang Anda dapatkan:
SubjectID <- rep(1:5, each = 3)
SubjectID
Outcome <- rep(0, 15)
Outcome[1] <- 1
Data <- data.frame(Outcome, SubjectID)
str(Data)
Data
library(lme4)
glmer(Outcome ~ 1 + (1|SubjectID), family="binomial", data = Data)
Dalam contoh ini, terdapat 5 mata pelajaran sehingga 4 di antaranya hanya memiliki 0 nilai hasil dan salah satunya memiliki nilai hasil yang menyertakan satu nilai 1. (Setiap mata pelajaran memiliki total 3 nilai hasil.)
Meskipun Anda memberikan setiap subjek dalam contoh yang dibuat-buat ini nilai 1 untuk nilai hasil pertamanya, Anda masih akan mendapatkan pesan kesalahan yang sama saat menyesuaikan model:
Outcome <- rep(0, 15)
Outcome[c(1, 4, 7, 10, 13)] <- 1
Namun, jika keempat subjek yang awalnya hanya memiliki 0 nilai diizinkan untuk menyimpan nilai-nilai ini dan subjek pertama menerima dua nilai 1, pesan kesalahan menghilang:
Outcome <- rep(0, 15)
Outcome[c(1,2)] <- 1
Setelah Anda lebih memahami pola nilai 0 dan 1 untuk variabel hasil di antara subjek studi Anda, hal lain yang dapat Anda coba adalah menyesuaikan model Anda dengan fungsi mixed_model () dari paket GLMMadaptive di R.
Untuk contoh kecil yang diberikan di sini, fungsi ini akan digunakan seperti ini:
library(GLMMadaptive)
m <- mixed_model(fixed = Outcome ~ 1,
random = ~ 1 | SubjectID,
data = Data,
family = binomial())
summary(m)
Isabella membuat beberapa poin bagus. Hal ini juga dapat terjadi jika variasi yang sangat sedikit di tingkat saluran. Mungkin saluran sangat mirip satu sama lain sehingga variansnya benar-benar mendekati nol dan oleh karena itu tidak diperlukan dalam model. Anda dapat menilai ini dengan menyesuaikan glm dan melihat apakah kesimpulannya serupa.
Karena penasaran, apakah kesalahan muncul saat Anda menggunakan penduga alternatif? Bisa jadi penaksir karena alasan tertentu terjebak pada singularitas. Anda hanya dapat mencoba berikut: mod.alt_est <- allFit(mod.detection_rand). Sebagai alternatif, Anda mungkin memerlukan solusi Bayesian untuk membantu mengatur estimasi dan memaksanya menjauh dari singularitas (coba blmepaket jika allFitfungsi tidak menghasilkan estimator yang berfungsi).