проблема подгонки модели dirlichet к моделируемым данным в JAGS, реализованной в R (бета-версия) - PullRequest
Новая
       91

проблема подгонки модели dirlichet к моделируемым данным в JAGS, реализованной в R (бета-версия)

0 голосов
/ 15 мая 2018

Допустим, у меня есть некоторые данные об относительной численности 3 видов.Вид 1 имеет отрицательную связь со среднегодовой температурой (mat).Данные относительной численности хранятся в матрице, r.spp.y: 

#Simulate some species count data.
y1 <- round(rnorm(100, 100, 5)) 
mat <- rnorm(length(y1), 10, 3)
y1 <- y1 + round(mat*-3)
y2 <- round(rnorm(100, 100, 5)) 
y3 <- round(rnorm(100, 100, 5)) 
spp.y <- as.matrix(data.frame(y1,y2,y3))
r.spp.y <- spp.y / rowSums(spp.y)

Я могу подогнать бета-регрессию для каждого вида по одному в JAGS, чтобы показать, что существует отрицательная корреляция между mat иОтносительная численность вида 1, но положительная корреляция между mat и относительной численностью вида 2 и 3 с использованием JAGS: 

library(runjags)
beta.model = "
model{
  #priors
  m0 ~ dnorm(0, 1.0E-3) #intercept prior
  m1 ~ dnorm(0, 1.0E-3) #      mat prior
  t0 ~ dnorm(0, .01)
  tau <- exp(t0) 

  #drop beta
  for (i in 1:N){
  y[i] ~ dbeta(p[i], q[i])
  p[i] <- mu[i] * tau
  q[i] <- (1 - mu[i]) * tau
  logit(mu[i]) <- m0 + m1 * mat[i]
  }

}#close model loop.
"
output <- list()
for(i in 1:ncol(r.spp.y)){
  beta.data <- list(y = r.spp.y[,i], mat = mat, N = nrow(r.spp.y))
  jags.beta <- run.jags(beta.model,
                        data=beta.data,
                        adapt = 200,
                        burnin = 1000,
                        sample = 1000,
                        n.chains = 3,
                        monitor = c('m0','m1'))
  output[[i]] <- summary(jags.beta)
}
names(output) <- c('species 1','species 2','species 3')

Это показывает, что у каждого есть значение перехвата и некоторая связь с mat.В частности, вид 1 имеет отрицательную связь с mat, тогда как виды 2 и 3 демонстрируют положительную связь с mat, что отражено в значениях параметра m1. 

$`species 1`
       Lower95      Median     Upper95        Mean          SD        Mode        MCerr MC%ofSD SSeff     AC.10     psrf
m0 -0.71012252 -0.66164521 -0.61324227 -0.66175203 0.025504290 -0.65828829 0.0025395302    10.0   101 0.4951512 1.005070
m1 -0.04608964 -0.04139645 -0.03671938 -0.04138209 0.002453942 -0.04148652 0.0002420001     9.9   103 0.5013018 1.006634

$`species 2`
       Lower95      Median     Upper95        Mean          SD        Mode        MCerr MC%ofSD SSeff     AC.10     psrf
m0 -0.74993189 -0.70794571 -0.67146076 -0.70773078 0.019593399 -0.70704999 0.0018992840     9.7   106 0.4913926 1.006231
m1  0.01479929  0.01849042  0.02212619  0.01844102 0.001845884  0.01825645 0.0001812411     9.8   104 0.4964613 1.006281

$`species 3`
       Lower95      Median     Upper95        Mean          SD        Mode       MCerr MC%ofSD SSeff     AC.10     psrf
m0 -0.71786559 -0.67713045 -0.63342213 -0.67726218 0.021507927 -0.67767494 0.001908073     8.9   127 0.4398105 1.024932
m1  0.01160666  0.01560139  0.01958145  0.01561193 0.002017519  0.01576153 0.000187649     9.3   116 0.4513887 1.026890

Я бы хотелчтобы приспособить все изобилие одновременно, используя многомерную модель.Это можно сделать с помощью распределения dirlichet, которое является многомерным обобщением бета-распределения.Для этого я устанавливаю модель dirlichet в JAGS, аналогичную бета-версии выше: 

dirlichet.model = "
model {
  #setup priors for each species
  for(j in 1:N.spp){
    m0[j] ~ dnorm(0, 1.0E-3) #intercept prior
    m1[j] ~ dnorm(0, 1.0E-3) #      mat prior
  }

  #implement dirlichet
  for(i in 1:N){
    for(j in 1:N.spp){
      logit(a0[i,j]) <- m0[j] + m1[j] * mat[i]
  }
  y[i,1:N.spp] ~ ddirch(a0[i,1:N.spp]) 
  }

} #close model loop.
"

Затем я настраиваю немного другой объект данных и запускаю модель: 

jags.data <- list(y = r.spp.y,mat = mat,N = nrow(r.spp.y), N.spp = ncol(r.spp.y))
jags.out <- run.jags(dirlichet.model,
                     data=jags.data,
                     adapt = 200,
                     burnin = 2000,
                     sample = 2000,
                     n.chains=3,
                     monitor=c('m0','m1'))
summary(jags.out)

Однакопараметры m0 и m1, возвращаемые никоим образом, не похожи на параметры, возвращаемые циклом бета-модели.Кажется, что все параметры отражают предоставленные (плоские) предыдущие распределения, и данные никак не ограничивают параметры.Эта модель не может уловить отрицательную связь между mat и относительной численностью вида 1. Она даже не слишком сильно ограничивает данные: 

> summary(jags.out)
         Lower95    Median  Upper95      Mean       SD      Mode     MCerr MC%ofSD SSeff         AC.10     psrf
m0[1] -49.005761  6.181091 63.01256  6.393293 28.43734  6.628375 0.5359853     1.9  2815 -0.0018363868 1.000590
m0[2] -45.469227  6.723323 67.85340  6.978617 28.52249  5.924171 0.5417058     1.9  2772  0.0129357542 1.000718
m0[3] -49.227870  6.233448 61.60683  6.599494 28.43137  4.742497 0.5279870     1.9  2900 -0.0057574414 1.000546
m1[1]  -1.869337 24.268826 64.16340 27.391725 19.13417 20.322548 0.4572891     2.4  1751  0.0109859218 1.006568
m1[2]  -1.465451 23.725515 64.21979 27.071232 18.98391 15.823296 0.4294849     2.3  1954 -0.0005437107 1.002363
m1[3]  -1.907332 24.031360 62.18261 27.017599 18.57743 18.597688 0.4290431     2.3  1875 -0.0086604393 1.001171

Стоит отметить, что если яв соответствии с моделью перехвата только для перехвата (без ковариации mat), она может хорошо фиксировать относительные содержания, однако, только если я использую лог-ссылку, а не логит-ссылку, использованную в бета-случае (пример здесь).Расширение регистрационной ссылки и добавление ковариаты mat не решает проблему.Что мне здесь не хватает?

...