解析出错（文本 = 粘贴（“~”，粘贴（nVal，折叠 =“/”）））：<文本>：2：0：在 R 中运行 nlme 包时输入意外结束

Question

我正在尝试使用 nlme 包将第二类分布的广义 beta 拟合到模拟的健康成本数据。

在测试数据集上运行以下代码：

包安装（如有必要）

install.packages("withr", dependencies = T)
library(withr)
with_makevars(c(PKG_CFLAGS ="-std=gnu99"), 
       install.packages("cubature"), assignment="+=") 
install.packages("GB2", dependencies = T)
    install.packages("nlme", dependencies = T)

---

# load packages
library(cubature)
library(GB2)
library(nlme)

# Binary independent variables
age <- rbinom(n=1000, size=1, prob=.3)
sex <- rbinom(n=1000, size=1, prob=.5)
trmt <- rbinom(n=1000, size=1, prob=.5)

# GB2 parameter equations
shape1 <- exp(rnorm(n=1000, mean=.1 + age/100 - sex/10 + trmt/10, sd=.3))
scale <- exp(rnorm(n=1000, mean=7 + age/50 + sex - trmt, sd=.5))
shape2 <- exp(rnorm(n=1000, mean=1.5 + age/100 + sex/10 - trmt/10, sd=.3))
shape3 <- exp(rnorm(n=1000, mean=.5 + age/100 - sex/10 - trmt/10, sd=.3))

# Outcome
y <- rgb2(1000, shape1, scale, shape2, shape3)

# Create test dataset
df <- data.frame(cbind(y,age,sex,trmt,shape1,scale,shape2,shape3))

# Fit GB2 distribution to data
gb2_fit <- nlme(y ~ scale*beta(shape2 + 1/shape1, shape3 - 1/shape1)/beta(shape2, shape3),
           # data = list(y=df_gb2_test[,1]),
           data = df,
           fixed = list(shape1 ~ age + sex + trmt, 
                        scale ~ age + sex + trmt, 
                        shape2 ~ age + sex + trmt, 
                        shape3 ~ age + sex + trmt),
           start = list(fixed = c(shape1 = 1.00, scale = 100, shape2 = 1.00, shape3 = 1.00)))

我得到错误：

Error in parse(text = paste("~", paste(nVal, collapse = "/"))) : 
  <text>:2:0: unexpected end of input
1: ~ 
   ^

任何想法我做错了什么？我似乎正确地使用了波浪号运算符。

Answer 1

我认为nlme 没有做你认为它做的事情。它做非线性最小二乘 混合模型；即，假设响应是高斯的，并且假设存在随机效应（也许您将此与更通用的 SAS PROC NLMIXED 混淆了？

library(bbmle)

## we need a version of the density function that takes a 'log' argument
dgb2B <- function(..., log=FALSE) {
  r <- GB2::dgb2(...)
  if (!log) r else log(r)
}

## don't include shape1, scale shape2, shape3 in the data, that confuses things
df2 <- df[,c("y","age","sex", "trmt")]


## fit homogeneous model
m1 <- mle2(y ~ dgb2B(shape1, scale, shape2, shape3),
     method="Nelder-Mead",
     trace=TRUE,
     data=df2,
     start = list(shape1 = 1.00, scale = 100, shape2 = 1.00, shape3 = 1.00))

## allow parameters to vary by group 
mle2(y ~ dgb2B(shape1, scale, shape2, shape3),
     ## parameters need to be in the same order!
     parameters=list(shape1 ~ age + sex + trmt,
                     scale ~ age + sex + trmt,
                     shape2 ~ age + sex + trmt,
                     shape3 ~ age + sex + trmt),
     method="Nelder-Mead",
     control=list(maxit=10000,
                  ## set parameter scales equal to magnitude
                  ## of starting values; each top-level parameter
                  ## has 4 associated values (intercept, + 3 cov effects)
                  parscale=rep(abs(coef(m1)), each=4)),
     trace=TRUE,
     data=df2,
     start = as.list(coef(m1))
)

对于它的价值，对于这个例子，你可以通过将八个不同的模型拟合到所有年龄 × 性别 × 治疗组来实现相同的目标（但我可以理解你的实际应用可能更复杂，即你可能只希望参数的子集在组之间变化，或者可能希望允许参数根据连续协变量变化。

如果您要尝试更难的问题，您可能希望在对数刻度上拟合参数。

Answer 2

之前还发生了一个错误：

y <- rgb2(1, shape1, scale, shape2, shape3)
Error in rgb2(1, shape1, scale, shape2, shape3) : 
  could not find function "rgb2"

您可能需要为此加载所需的包：

https://www.rdocumentation.org/packages/gamlss.dist/versions/5.3-2/topics/GB2

它似乎在library(gamlss.dist)