这是为 JAGS 脚本实现基本宏替换系统的一个示例。
系统说明
- 定义一个函数，该函数将脚本的任何可选元素作为参数。
- 对于因参数值而异的脚本的任何方面，请记录宏标记。这应该是一些独特的文本。以一些符号开始和结束可能有助于明确这一点。
- 对于每个宏标记，包括函数参数的所有可能组合下的替换字符串应该是什么的代码。
- 根据参数将宏标记替换为适当的放置文本。
- 下面的代码提供了一个如何指定宏以及如何将宏应用于原始脚本的示例（欢迎提出更优雅地执行此操作的建议）。
- 该函数返回一个替换的模型字符串，如有必要，可以将其传递给
textConnection函数以与 rjags 一起使用。
我喜欢这个系统有几个原因：
- 原始脚本易于阅读
- 生成的脚本有适当的缩进
- 您不必担心与出现在同一行的命令相关的错误
例子
具体来说，此示例旨在允许用户拟合特定类型的多级非线性模型。它被设计为允许三种函数形式：二参数幂、三参数幂和三参数指数。
宏部分将宏构造为列表列表。顶级列表包含每个宏标记的一个元素。对于每个宏标记，都有宏标记文本和条件替换文本。
最后，for 循环将所有宏替换应用到原始脚本。
见下文（需要滚动）：
jags_model <- function (f=c('power2', 'power3', 'exp3')) {
f <- match.arg(f)
# raw script
script <-
"model {
# Model
for (i in 1:length(y)) {
$FUNCTION
y[i] ~ dnorm(mu[i], tau[subject[i]])
}
# Random coefficients
for (i in 1:N) {
theta1[i] ~ dnorm(theta1.mu, theta1.tau)T(0, 1000)
theta2[i] ~ dnorm(theta2.mu, theta2.tau)T(-10, 0)
$THETA3DISTRIBUTION
sigma[i] ~ dnorm(sigma.mu, sigma.tau)T(0, 100);
tau[i] <- 1/(sigma[i]^2)
}
theta1.mu ~ dunif(0, 100)
theta2.mu ~ dunif(-2, 0)
$THETA3PRIOR.MU
sigma.mu ~ dunif(0, 20)
theta1.sigma ~ dunif(0, 100)
theta2.sigma ~ dunif(0, 2)
$THETA3PRIOR.SIGMA
sigma.sigma ~ dunif(0, 10)
# Transformations
theta1.tau <- 1/(theta1.sigma^2)
theta2.tau <- 1/(theta2.sigma^2)
$THETA3.TAU
sigma.tau <- 1/(sigma.sigma^2)
}"
# define macros
macros <- list(list("$FUNCTION",
switch(f,
power2="mu[i] <- theta1[subject[i]] * pow(trial[i], theta2[subject[i]])",
power3="mu[i] <- theta1[subject[i]] * pow(trial[i], theta2[subject[i]]) + theta3[subject[i]];",
exp3="mu[i] <- theta1[subject[i]] * exp(theta2[subject[i]] * (trial[i] - 1)) + theta3[subject[i]];") ),
list("$THETA3DISTRIBUTION",
switch(f,
power3=, exp3= "theta3[i] ~ dnorm(theta3.mu, theta3.tau)T(0, 1000)",
power2="") ),
list("$THETA3PRIOR.MU",
switch(f,
power3=, exp3= "theta3.mu ~ dunif(0, 100)",
power2="") ),
list("$THETA3PRIOR.SIGMA",
switch(f,
power3=, exp3= "theta3.sigma ~ dunif(0, 100)",
power2="") ),
list("$THETA3.TAU",
switch(f,
power3=, exp3= "theta3.tau <- 1/(theta3.sigma^2)",
power2="") )
)
# apply macros
for (m in seq(macros)) {
script <- gsub(macros[[m]][1], macros[[m]][2], script, fixed=TRUE)
}
script
}
示范
因此，我们可以生成处理后的 JAGS 模型
x <- jags_model(f='power3')
如果我们想查看生成的模型，我们可以这样做
cat(x)
这导致
model {
# Model
for (i in 1:length(y)) {
mu[i] <- theta1[subject[i]] * pow(trial[i], theta2[subject[i]]) + theta3[subject[i]];
y[i] ~ dnorm(mu[i], tau[subject[i]])
}
# Random coefficients
for (i in 1:N) {
theta1[i] ~ dnorm(theta1.mu, theta1.tau)T(0, 1000)
theta2[i] ~ dnorm(theta2.mu, theta2.tau)T(-10, 0)
theta3[i] ~ dnorm(theta3.mu, theta3.tau)T(0, 1000)
sigma[i] ~ dnorm(sigma.mu, sigma.tau)T(0, 100);
tau[i] <- 1/(sigma[i]^2)
}
theta1.mu ~ dunif(0, 100)
theta2.mu ~ dunif(-2, 0)
theta3.mu ~ dunif(0, 100)
sigma.mu ~ dunif(0, 20)
theta1.sigma ~ dunif(0, 100)
theta2.sigma ~ dunif(0, 2)
theta3.sigma ~ dunif(0, 100)
sigma.sigma ~ dunif(0, 10)
# Transformations
theta1.tau <- 1/(theta1.sigma^2)
theta2.tau <- 1/(theta2.sigma^2)
theta3.tau <- 1/(theta3.sigma^2)
sigma.tau <- 1/(sigma.sigma^2)
}