R函数代码语法错误排查求助：封装为函数前可正常运行

我把一段原本能正常运行的代码封装后出现了语法错误，却找不到问题所在。作为R语言新手，这是我写完测试小函数后第一个实用函数，代码如下：

#model bzip
{
  #Likelihood:
  for(i in 1:n)
  {
    Y1[i] ~ dpois(mu1[i])
    Y2[i] ~ dpois(mu2[i])
    mu1[i] <- (1-u[i, 1])(1-u[i, 3])
    (lambda0 + lambda1)
    mu2[i] <- (1-u[i, 1])(1-u[i, 2])
    (lambda0+lambda2)
    u[i, 1:4] ~ dmulti(p[], 1)
  }
  for(i in 1:3)
  {
    log(q[i]) <-alpha[i]
    p[i] <- q[i]*p[4]
  }
  p[4] <- 1/(q[1]+q[2]+q[3]+1)
  log(lambda0) <-beta[1]
  log(lambda1) <-beta[2]
  log(lambda2) <-beta[3]
  zdp <- p[1]+p[2]*(exp(-lambda0 -lambda1))
  +p[3]*(exp(-lambda0 -lambda2)) +
    p[4]*(exp(-lambda0 -lambda1 -lambda2))
  #P(Y1=Y2=0)
  zdp1 <- p[1]+p[3]+(p[2]+p[4])*(exp(-lambda0
                                     -lambda1)) #P(Y1=0)
  zdp2 <- p[1]+p[2]+(p[3]+p[4])*(exp(-lambda0
                                     -lambda2))#P(Y2=0)
  e1 <- (p[2] +p[4])*(lambda0 + lambda1) # E(Y1)
  e2 <- (p[3] +p[4])*(lambda0 + lambda2) # E(Y2)
  #Priors:
  for(j in 1:3)
  {
    beta[j]~ dnorm(0, 0.001)
    alpha[j]~dnorm(0, 0.001)
  }
}

错误点分析 #

• 乘法运算符缺失：mu1[i]和mu2[i]的表达式被换行拆分，没有用*连接各因子，导致语法断裂，需将所有因子用*串联。
• 行首加号语法错误：计算zdp时，+号放在行首会被识别为新表达式，应将+留在上一行末尾，保持表达式连续性。
• 未声明数组长度：q数组使用前未定义长度，Stan无法识别索引，需提前声明vector[3] q;。
• 分布符号空格不规范：~前后缺少空格，虽部分解析器兼容，但统一加空格可避免语法解析问题。

修正后的代码 #

#model bzip
{
  #Likelihood:
  for(i in 1:n)
  {
    Y1[i] ~ dpois(mu1[i])
    Y2[i] ~ dpois(mu2[i])
    mu1[i] <- (1-u[i, 1])*(1-u[i, 3])*(lambda0 + lambda1)
    mu2[i] <- (1-u[i, 1])*(1-u[i, 2])*(lambda0 + lambda2)
    u[i, 1:4] ~ dmulti(p[], 1)
  }
  
  # 声明q数组长度
  vector[3] q;
  for(i in 1:3)
  {
    log(q[i]) <- alpha[i]
    p[i] <- q[i]*p[4]
  }
  p[4] <- 1/(q[1]+q[2]+q[3]+1)
  
  log(lambda0) <- beta[1]
  log(lambda1) <- beta[2]
  log(lambda2) <- beta[3]
  
  # 修正加号位置，保持表达式连续
  zdp <- p[1] + p[2]*exp(-lambda0 - lambda1) + 
         p[3]*exp(-lambda0 - lambda2) + 
         p[4]*exp(-lambda0 - lambda1 - lambda2)
  #P(Y1=Y2=0)
  zdp1 <- p[1] + p[3] + (p[2]+p[4])*exp(-lambda0 - lambda1) #P(Y1=0)
  zdp2 <- p[1] + p[2] + (p[3]+p[4])*exp(-lambda0 - lambda2) #P(Y2=0)
  
  e1 <- (p[2] + p[4])*(lambda0 + lambda1) # E(Y1)
  e2 <- (p[3] + p[4])*(lambda0 + lambda2) # E(Y2)
  
  #Priors:
  for(j in 1:3)
  {
    beta[j] ~ dnorm(0, 0.001)
    alpha[j] ~ dnorm(0, 0.001)
  }
}

R中调用的封装示例 #

注意这段是Stan模型代码，不是纯R函数。若要在R中调用，需用rstan包封装为可执行函数：

library(rstan)

# 封装为可调用的R函数
fit_bzip_model <- function(data) {
  model_code <- "
  #model bzip
  {
    #Likelihood:
    for(i in 1:n)
    {
      Y1[i] ~ dpois(mu1[i])
      Y2[i] ~ dpois(mu2[i])
      mu1[i] <- (1-u[i, 1])*(1-u[i, 3])*(lambda0 + lambda1)
      mu2[i] <- (1-u[i, 1])*(1-u[i, 2])*(lambda0 + lambda2)
      u[i, 1:4] ~ dmulti(p[], 1)
    }
    
    vector[3] q;
    for(i in 1:3)
    {
      log(q[i]) <- alpha[i]
      p[i] <- q[i]*p[4]
    }
    p[4] <- 1/(q[1]+q[2]+q[3]+1)
    
    log(lambda0) <- beta[1]
    log(lambda1) <- beta[2]
    log(lambda2) <- beta[3]
    
    zdp <- p[1] + p[2]*exp(-lambda0 - lambda1) + 
           p[3]*exp(-lambda0 - lambda2) + 
           p[4]*exp(-lambda0 - lambda1 - lambda2)
    zdp1 <- p[1] + p[3] + (p[2]+p[4])*exp(-lambda0 - lambda1)
    zdp2 <- p[1] + p[2] + (p[3]+p[4])*exp(-lambda0 - lambda2)
    
    e1 <- (p[2] + p[4])*(lambda0 + lambda1)
    e2 <- (p[3] + p[4])*(lambda0 + lambda2)
    
    for(j in 1:3)
    {
      beta[j] ~ dnorm(0, 0.001)
      alpha[j] ~ dnorm(0, 0.001)
    }
  }
  "
  # 拟合模型
  fit <- stan(model_code = model_code, data = data)
  return(fit)
}

# 使用示例（需传入包含n、Y1、Y2、u等的数据集）
# fit <- fit_bzip_model(data = list(n = 100, Y1 = ..., Y2 = ..., u = ...))

内容的提问来源于stack exchange，提问作者Fakhri Albarz