R2OpenBUGS error in parse: unexpected symbol

Question

I am currently trying to fit a WinBUGS model using R2OpenBUGS. However, I get the following error:

Error in parse(text = x, srcfile = src): <text>:5:44: unexpected symbol
4:         for(j in 1:P){
5:             y[i,j] ~ dnorm(mu[i,j], psi[j])I

Apparently, this has something to do with the inappropriate R syntax on line 5:

y[i, j] ~ dnorm(mu[i, j], psi[j])I(thd[j, z[i, j]], thd[j, z[i, j] + 1])

Changing the line to appropriate R syntax:

y[i, j] ~ dnorm(mu[i, j], psi[j]) * I(thd[j, z[i, j]], thd[j, z[i, j] + 1])

will remove the error, but cause the model run to crash with the following log message:

"invalid or unexpected token scanned error pos 100 (error on line 5)"

I am trying to discern what's causing the problem. Any help is appreciated.

Full model:

model <- function() {
    for(i in 1:N){
        #measurement equation model
        for(j in 1:P){
            y[i,j] ~ dnorm(mu[i,j], psi[j]) * I(thd[j, z[i, j]], thd[j, z[i, j] + 1])
            ephat[i,j] <- y[i,j] - mu[i,j]
        }
        mu[i,1] <- eta[i]
        mu[i,2] <- lam[1] * eta[i]
        mu[i,3] <- xi[i,1]
        mu[i,4] <- lam[2] * xi[i,1]
        mu[i,5] <- lam[3] * xi[i,1]
        mu[i,6] <- lam[4] * xi[i,1]
        mu[i,7] <- lam[5] * xi[i,1]
        mu[i,8] <- lam[6] * xi[i,1]
        mu[i,9] <- lam[7] * xi[i,1]
        mu[i,10] <- xi[i,2]
        mu[i,11] <- lam[8] * xi[i,2]
        mu[i,12] <- lam[9] * xi[i,2]
        mu[i,13] <- lam[10] * xi[i,2]
        mu[i,14] <- lam[11] * xi[i,2]
        mu[i,15] <- lam[12] * xi[i,2]
        mu[i,16] <- xi[i,3]
        mu[i,17] <- lam[13] * xi[i,3]
        mu[i,18] <- lam[14] * xi[i,3]
        mu[i,19] <- xi[i,4]
        mu[i,20] <- lam[15] * xi[i,4]
        mu[i,21] <- lam[16] * xi[i,4]
        mu[i,22] <- lam[17] * xi[i,4]
        mu[i,23] <- lam[18] * xi[i,4]
        mu[i,24] <- lam[19] * xi[i,4]
        mu[i,25] <- lam[20] * xi[i,4]
        mu[i,26] <- lam[21] * xi[i,4]

        #structural equation model
        xi[i,1:4] ~ dmnorm(u[1:4], phi[1:4,1:4])
        eta[i] ~ dnorm(nu[i], psd)
        nu[i] <- gam[1] * xi[i,1] + gam[2] * xi[i,2] + gam[3] * xi[i,3] + gam[4] * xi[i,4]
        dthat[i] <- eta[i] - nu[i]
    }# end of i

    for(i in 1:4){
        u[i] <- 0.0
    }

    #priors on loadings and coefficients
    var.lam[1] <- 4.0 * psi[2]
    var.lam[2] <- 4.0 * psi[4]
    var.lam[3] <- 4.0 * psi[5]
    var.lam[4] <- 4.0 * psi[6]
    var.lam[5] <- 4.0 * psi[7]
    var.lam[6] <- 4.0 * psi[8]
    var.lam[7] <- 4.0 * psi[9]
    var.lam[8] <- 4.0 * psi[11]
    var.lam[9] <- 4.0 * psi[12]
    var.lam[10] <- 4.0 * psi[13]
    var.lam[11] <- 4.0 * psi[14]
    var.lam[12] <- 4.0 * psi[15]
    var.lam[13] <- 4.0 * psi[17]
    var.lam[14] <- 4.0 * psi[18]
    var.lam[15] <- 4.0 * psi[20]
    var.lam[16] <- 4.0 * psi[21]
    var.lam[17] <- 4.0 * psi[22]
    var.lam[18] <- 4.0 * psi[23]
    var.lam[19] <- 4.0 * psi[24]
    var.lam[20] <- 4.0 * psi[25]
    var.lam[21] <- 4.0 * psi[26]
    
    for(i in 1:21){
        lam[i] ~ dnorm(0.8,var.lam[i])
    }

    var.gam <- 4.0 * psd
    gam[1] ~ dnorm(0.6,var.gam)
    gam[2] ~ dnorm(0.6,var.gam)
    gam[3] ~ dnorm(0.4,var.gam)
    gam[4] ~ dnorm(0.4,var.gam)

    #priors on precisions
    for(j in 1:P){
        psi[j] ~ dgamma(10,8)
        sgm[j] <- 1/psi[j]
    }
    psd ~ dgamma(10,8)
    sgd <- 1/psd
    phi[1:4,1:4] ~ dwish(R[1:4,1:4], 30)
    phx[1:4,1:4] <- inverse(phi[1:4, 1:4])

} #end of model