¿Cómo condicionalmente ejecutar el elemento de script JAGS basado en la variable suministrada por el usuario?

Question

Antecedentes: a menudo me encuentro que cuando se trabaja con EXIGENCIAS de que termino con un montón de EXIGENCIAS de secuencias de comandos, como voy a explorar diferentes modelos. Después de un tiempo me puedan asentarse en un conjunto de modelos que te voy a denunciar, pero luego está el tema de mantener la configuración coherente a través de estos modelos que están destinados a ser consistente. Por ejemplo, los priores de los parámetros que son comunes deben ser coherentes, o los nombres de las variables para los parámetros deben ser consistentes. Mientras que la coherencia no es tan difícil de lograr, es a menudo un poco molesto. Por ejemplo, puedo decidir que quiero llamar a algo theta en lugar de la beta, y ahora todos los scripts de actualización. Así que estoy tratando de averiguar cómo puedo tener una secuencia de comandos que pueden servir para múltiples propósitos.

Ejemplo: He aquí un ejemplo simple. Supongamos que tenemos dos modelos, uno con el tiempo como un predictor y uno sin. Yo de suministro a exigencias de una variable lógica timevariable a indicar si quiero el modelo de tiempo o no. Por lo tanto, me gustaría ser capaz de escribir algo como esto.

for (i in 1:length(y)) {
    if (timevariable) {
        mu[i] <- beta1 + beta2 * time[i]
    } else {
        mu[i] <- beta1
    }
    y[i]  ~ dnorm(mu[i], tau)
}

Sin embargo, si-entonces no es parte de la lanugage. Tenga en cuenta que no estoy tratando de incorporar un "si-entonces" basado en los valores estimados de los parámetros. Más bien, lo que quiero es tener un script que puede tomar flexible de entrada.

Pregunta

Lo que es una buena manera de ejecutar condicionalmente partes de un ENTRECORTADO script basado en un usuario proporcionado variable?

Reflexiones iniciales

• Poner un nombre de macro en el texto y el uso de un texto de reemplazo de la función en R como gsub a reemplazar el nombre de la macro con el texto deseado. Supongo que esto iba a funcionar, pero me parece un poco de un hack. Sería bueno tener todo dentro de la secuencia de comandos.

Answer 1

Aquí es un ejemplo de la implementación de un macro básico de sustitución del sistema de EXIGENCIAS de secuencias de comandos.

Explicación del sistema

• Definir una función que toma como argumentos opcionales de elementos de la secuencia de comandos.
• Para cualquiera de los aspectos de la secuencia de comandos que varían a través de los valores de los argumentos, grabar una macro token. Esto debe ser algo único texto. Empezando y terminando con algunos símbolos pueden ayudar a hacer esto sin ambigüedades.
• Para cada macro token de incluir el código de lo que la cadena de reemplazo debe ser bajo todas las combinaciones posibles de los argumentos de la función.
• Reemplace la macro fichas con la colocación adecuada de texto basado en argumentos.
• El código a continuación se proporciona un ejemplo de cómo las macros se pueden especificar y cómo aplicar las macros para el raw de secuencia de comandos (sugerencias para hacerlo más elegante son bienvenidos).
• La función devuelve un sustituye el modelo de cadena que puede, si es necesario pasarse a una textConnection función para su uso con rjags.

Me gusta este sistema para un par de razones:

• la cruda script es fácil de leer
• la secuencia de comandos resultante tiene sangría adecuada
• usted no tiene que preocuparse acerca de los errores relacionados con los comandos que aparecen en la misma línea

Ejemplo

Específicamente, en este ejemplo, tiene por objetivo permitir que el usuario ajuste a un tipo particular de multinivel no lineales del modelo. Está diseñado para permitir tres formas funcionales: una de dos parámetros de potencia, tres parámetros de potencia, y un de tres parámetro exponencial.

La macro estructuras de sección de macros como una lista de listas. El nivel superior de la lista contiene un elemento por cada macro token. Para cada macro token, no es la macro token de texto, y el condicional texto de reemplazo.

Finalmente, un bucle for se aplica a todos los macro reemplazos para el raw de secuencia de comandos.

Vea a continuación (el desplazamiento es necesario):

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
}

Demostración

Por lo tanto, podemos producir el procesado ENTRECORTADO modelo con

x <- jags_model(f='power3')

Y si queremos ver el modelo resultante, que podemos hacer

cat(x)

que se traduce en

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

Answer 2

En este caso específico se podría establecer la beta2 plazo a cero

for (i in 1:length(y)) {
    mu[i] <- beta1 + indicator[i] * beta2 * time[i]
    y[i]  ~ dnorm(mu[i], tau)
}

donde indicator[] es un vector que es uno de los puntos de datos que se desea modelar con beta2 y 0 en caso contrario. También podría utilizar un escalar para cambiar el modelo como un todo. Sin embargo, este enfoque es probablemente menos eficaz porque todavía muestras de un no deseado del parámetro de la distribución. He utilizado el gsub enfoque de antes, pero estoy de acuerdo en que no es bonita.

Si el modelo está escrito y ejecutado en un R secuencia de comandos puede ser manipulado como cualquier otra cadena:

modelstring <- paste("
model {  
  for ( i in 1:N ) {
",ifelse(timevariable,"
   mu[i] <- beta1 + beta2 * time[i]
","mu[i] <- beta1"),"
  }
  y[i]  ~ dnorm(mu[i], tau) 
}
")
writeLines(modelstring,con=filename)