1

下面的 R 代码来自 Coursera 课程，但它似乎不起作用。我收到一个错误，表明它y[1]与它的父母不兼容，但鉴于它的值为 1 并且它的父母形成了一个 n 等于 4 的二项分布，我看不出它是如何不兼容的。

以下命令显示y[1] = 1and n[1] = 4：

head(str(data_jags))
List of 6
 $ Age            : num [1:712] 30 30 30 30 30 30 30 30 30 30 ...
 $ OMElow         : num [1:712] 1 1 1 1 1 1 1 1 1 1 ...
 $ Loud           : num [1:712] 35 35 40 40 45 45 50 50 55 55 ...
 $ Noiseincoherent: num [1:712] 0 1 0 1 0 1 0 1 0 1 ...
 $ y              : int [1:712] 1 4 0 1 2 2 3 4 3 2 ...
 $ n              : int [1:712] 4 5 3 1 4 2 3 4 3 2 ...

完整的错误是：

Error in jags.model(textConnection(mod_string), data = data_jags, inits = inits,  : 
  Error in node y[1]
Node inconsistent with parents

完整的代码是：

library("rjags")
library("coda")
library("MASS")
data("OME")
?OME # background on the data
head(OME)

any(is.na(OME)) # check for missing values
dat = subset(OME, OME != "N/A") # manually remove OME missing values identified with "N/A"
dat$OME = factor(dat$OME)
str(dat)

plot(dat$Age, dat$Correct / dat$Trials )
plot(dat$OME, dat$Correct / dat$Trials )
plot(dat$Loud, dat$Correct / dat$Trials )
plot(dat$Noise, dat$Correct / dat$Trials )

mod_glm = glm(Correct/Trials ~ Age + OME + Loud + Noise, data=dat, weights=Trials, family="binomial")
summary(mod_glm)

plot(residuals(mod_glm, type="deviance"))
plot(fitted(mod_glm), dat$Correct/dat$Trials)


X = model.matrix(mod_glm)[,-1] # -1 removes the column of 1s for the intercept
head(X)

mod_string = " model {
    for (i in 1:length(y)) {
y[i] ~ dbin(phi[i], n[i])
logit(phi[i]) = b0 + b[1]*Age[i] + b[2]*OMElow[i] + b[3]*Loud[i] + b[4]*Noiseincoherent[i]
}

b0 ~ dnorm(0.0, 1.0/5.0^2)
for (j in 1:4) {
b[j] ~ dnorm(0.0, 1.0/4.0^2)
}

} "

data_jags = as.list(as.data.frame(X))
data_jags$y = dat$Correct # this will not work if there are missing values in dat (because they would be ignored by model.matrix). Always make sure that the data are accurately pre-processed for JAGS.
data_jags$n = dat$Trials
str(data_jags) # make sure that all variables have the same number of observations (712).

set.seed(54)

params = c("b0", "b")

inits = function() {
  inits = list("b0" = rnorm(1, 0.0, 100.0), "b"=rnorm(4, 0.0, 100.0))
}

mod = jags.model(textConnection(mod_string), data=data_jags, inits=inits, n.chains=3)

update(mod, 1000)

mod.sim = coda.samples(model=mod, variable.names=params, n.iter=5e3)

mod.c.sim = do.call(rbind, mod.sim)
4

0 回答 0