Как указать более сильные приоры в MCMCglmm?

Я несколько недель работал с R-пакетом MCMCglmm. Работаю с ним впервые. Я прочитал много статей и руководств для лучшего понимания, но я не могу решить проблему, которая у меня есть:

Это часть моих данных (только для одного человека):

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 и Thin, но в конце всегда пишет одну и ту же ошибку. Однако, если я уменьшу эти значения, я не получу никакой ошибки, но мои графики не такие, какими должны быть.


r mcmc
person Luz García-Longoria Batanete    schedule 29.11.2017    source источник
comment
Это ваш код? Кажется, что каждая строка следует другому соглашению, из-за чего создается впечатление, что вы частично позаимствовали ее у кого-то. Может быть, что кто-то еще может помочь вам? Если это не удастся, в институте обязательно найдется кто-нибудь, кто может помочь вам. Это кажется довольно сложным исследовательским вопросом, и кому-то без конкретных знаний будет сложно помочь.   -  person AkselA    schedule 29.11.2017

Ответы (1)


arrow_upward
0
arrow_downward

Во-первых, в спецификации модели вы не указали «предыдущий» как ваш предыдущий список, поэтому спецификация вашей модели должна выглядеть так:

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

Надеюсь, это поможет!

person Z. Radai    schedule 04.06.2018