Как написать файл модели для бинома JAGS с помощью функции logit - PullRequest
Новая
       26

Как написать файл модели для бинома JAGS с помощью функции logit

2 голосов
/ 25 апреля 2019

Я работаю над назначением, используя JAGS для моделирования биномиального распределения, параметр p которого является функцией другой переменной d.

Вот что я пытаюсь сделать:

  1. генерирует 10000 образцов из апостериорного для двух параметров альфа / бета
  2. производит выборки из заданного прогнозируемого числа успеха, когда dist = 25 для 100 попыток
  3. рассчитать 95 вероятных интервалов для успеха на расстоянии 25 футов

Я написал модель, но выдает ошибку.

Ниже приведен код, который я уже пробовал 

#R-code
distance=seq(from=2,to=20,by=1)
Ntrys=c(1443,694,455,353,272,256,240,217,200,237,202,192,174,167,201,195,191,147,152)
Nsucc=c(1346,577,337,208,149,136,111,69,67,75,52,46,54,28,27,31,33,20,24)

psucc=Nsucc/Ntrys

glm1.data=list(N=19, Nsucc=Nsucc,psucc=psucc,distance=distance)

glm1.model=jags.model("glm1.model",glm1.data,n.chains=2)

glm1.samps=coda.samples(glm1.model, variable.names=c("alpha", "beta"), 1e5)

#model file
model{ 
    for (i in 1:N){
            Nsucc[i] ~ dbern(psucc[i])
            log((psucc[i])/(1-psucc[i])) <- alpha + beta*(distance[i])
    }
    alpha ~ dunif(-10,10)
    beta ~ dunif(-10,10)
}

Я получаю ошибку

Ошибка в jags.model ("glm1.model", glm1.data, n.chains = 2):
ОШИБКА ВРЕМЕНИ:
Ошибка компиляции в строке 4.
pmiss [1] является логическим узлом и его нельзя наблюдать

Я не думаю, что файл модели даже настроен на то, что я пытаюсь сделать.

1 Ответ

0 голосов
/ 25 апреля 2019

Вам не нужно вычислять вероятности за пределами rjags, но вы можете использовать функцию биномиального распределения dbin(p,N), которая принимает аргументы, p, вероятность успеха и N, количество попыток , Кроме того, в качестве функции связи может использоваться функция logit.

Обновленная функция модели будет 

mod <-
"model{ 
    # likelihood
    for (i in 1:N){
            Nsucc[i] ~ dbin(p[i], Ntrys[i])
            logit(p[i]) <- alpha + beta*distance[i]
    }
    # priors
    alpha ~ dunif(-10,10)
    beta ~ dunif(-10,10)

}"

Прогнозы могут быть сгенерированы с учетом некоторого значения предикторов путем добавления значений предикторов к данным и добавления соответствующего числа NA к вектору результатов. Таким образом, данные, переданные в rjags, становятся 

glm1.data <- list(N=20, Nsucc=c(Nsucc, NA), Ntrys=c(Ntrys, 100), distance=c(distance, 25))

Затем скомпилируйте и запустите модель 

# set.seed so sampling is reproducible
library(rjags)
load.module("glm")

glm1.model <- jags.model(textConnection(mod), glm1.data, 
                         n.chains=2,
                         inits=list(.RNG.name="base::Wichmann-Hill",
                                    .RNG.seed=1))
update(glm1.model, n.iter = 1000, progress.bar="none")

# sample: monitor the unknown predictions, Nsucc[20], p[20]
glm1.samps <- coda.samples(glm1.model, variable.names=c("alpha", "beta", "Nsucc[20]", "p[20]"), 1e5)

Затем вы можете генерировать интервалы из квантилей 

s <- summary(glm1.samps)
s$quantiles

или самый высокий интервал плотности 

library(HDInterval)
hdi(glm1.samps)

(просто ради интереса, сравните коэффициенты от glm: summary(glm(cbind(Nsucc, Ntrys-Nsucc) ~ distance, family=binomial)))

...