Как определить более сильные приоры в MCMCglmm?
Я работал в течение нескольких недель с пакетом MCMCglmm R. Это первый раз, когда я работаю с ним. Я прочитал много статей и руководств для лучшего понимания, но я не могу решить проблему, которая у меня есть:
Это часть моих данных (только для одного человека):
Species Individual Lineage Prevalence day breeding Year phylo
Aegithalos_caudatus I1 SGS1 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 CARDUEL1 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 CARDUEL2 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 CARDUEL3 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 DELURB1 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 DELURB2 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 DELURB3 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 DELURB5 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 DURB6 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 GRW2 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 GRW4 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 GRW9 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 RTSR1 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 ROBIN1 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 GRW9 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 HIPOL1 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 COLL1 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 GRW11 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 PADOM_5 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 PADOM01 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 PADOM08 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 PADOM22 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 PAHIS_01 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 CCF2 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 SYMEL1 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 SYAT_05 0 125 Yes i2010 Aegithalos_caudatus
Aegithalos_caudatus I1 TURDUS3 0 125 Yes i2010
У меня 815 особей, которые принадлежат к 23 различным видам. У каждого человека есть 24 записи (одна на паразитическую линию). Я также хочу принять во внимание филогенетическое дерево-хозяин (птицу). Основная цель моих анализов состоит в том, чтобы проверить, могут ли распространенность (зараженный =0 или не зараженный =1) быть затронутой видами хозяина или происхождением.
Итак, я запускаю свою модель MCMCglmm и получаю сообщения об ошибках моих приоров:
data= read.table("data.txt", header = T)
phylo<-read.nexus("tree.nex")
inv.phylo<-inverseA(phylo,nodes="TIPS",scale=TRUE)
prior.ex.<- list(G = list(G1 = list(V = 1, nu = 2, alpha.mu = 0, alpha.V = 10000),
G2 = list(V = 1, nu = 2, alpha.mu = 0, alpha.V = 10000),
G3 = list(V = 1, nu = 2, alpha.mu = 0, alpha.V = 10000),
G4 = list(V = 1, nu = 2, alpha.mu = 0, alpha.V = 10000),
R = list(R1 = list(V =1, nu = 0.02, fix = TRUE),
R2 = list(V =1, nu = 0.02, fix = TRUE))))
model <- MCMCglmm(fixed = Prevalence ~ day + breeding,
random = ~ Year + Lineage + Individual + Individual:Species,
data = data,
family = "categorical", ginverse=list(phylo=inv.phyl$Ainv), nitt=300000, burnin=60000, thin=200, verbose = TRUE)
# MCMC iteration = 0
# Acceptance ratio for liability set 1 = 0.000658
# MCMC iteration = 1000
# Acceptance ratio for liability set 1 = 0.443318
# MCMC iteration = 2000
# Acceptance ratio for liability set 1 = 0.428678
# MCMC iteration = 3000
# Acceptance ratio for liability set 1 = 0.438693
# MCMC iteration = 4000
# Acceptance ratio for liability set 1 = 0.437736
# MCMC iteration = 5000
# Acceptance ratio for liability set 1 = 0.440917
# MCMC iteration = 6000
# Acceptance ratio for liability set 1 = 0.430484
# MCMC iteration = 7000
# Acceptance ratio for liability set 1 = 0.434346
# MCMC iteration = 8000
# Acceptance ratio for liability set 1 = 0.428129
# MCMC iteration = 9000
# Acceptance ratio for liability set 1 = 0.435863
# MCMC iteration = 10000
# Acceptance ratio for liability set 1 = 0.437513
# MCMC iteration = 11000
# Acceptance ratio for liability set 1 = 0.437331
# MCMC iteration = 12000
# Acceptance ratio for liability set 1 = 0.431787
# MCMC iteration = 13000
# Acceptance ratio for liability set 1 = 0.429239
# MCMC iteration = 14000
# Acceptance ratio for liability set 1 = 0.435133
# MCMC iteration = 15000
# Acceptance ratio for liability set 1 = 0.436029
# Error in MCMCglmm(fixed = Prevalence ~ day + breeding, random = ~Year + :
# Mixed model equations singular: use a (stronger) prior
Я пытался изменить Burnin, Nitt и Slim, но в конце он всегда говорит ту же ошибку. Однако, если я уменьшу это значение, я не получу никакой ошибки, но мои графики не будут такими, какими они должны быть.
1 ответ
Во-первых, в спецификации модели вы пропустили задание "предыдущий" в качестве предыдущего списка, поэтому спецификация модели должна выглядеть следующим образом:
model <- MCMCglmm(fixed = Prevalence ~ day + breeding,
random = ~ Year + Lineage + Individual + Individual:Species,
prior=prior.ex, data = data,
family = "categorical", ginverse=list(phylo=inv.phyl$Ainv), nitt=300000, burnin=60000, thin=200, verbose = TRUE)
Кроме того, я думаю, что ваш предыдущий не полностью работоспособен, попробуйте это:
prior.ex<- list(G = list(G1 = list(V = 1, nu = 2, alpha.mu = 0, alpha.V = 10000),
G2 = list(V = 1, nu = 2, alpha.mu = 0, alpha.V = 10000),
G3 = list(V = 1, nu = 2, alpha.mu = 0, alpha.V = 10000),
G4 = list(V = 1, nu = 2, alpha.mu = 0, alpha.V = 10000)),
R = list(V=1, fix=1))
Надеюсь, это поможет!