在R中拟合广义非线性模型

Question

我希望拟合以下广义非线性模型：Probit(G)=K+1/Sigma*(Temp-T0)*Time。作为Naive Model，我通过qnorm(G)创建了Probits(G)，之后拟合了非线性模型。但是我想用logit链接来拟合R中的glm函数类似于Nonlinear Model。

如何在R中使用logit链接拟合这样的广义非线性模型？

Data <-
  structure(list(Temp = c(23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L,
  27L, 27L, 27L, 27L, 27L, 27L, 33L, 33L, 33L, 33L, 33L, 33L, 33L,
  35L, 35L, 35L, 35L, 35L), Time = c(144L, 168L, 192L, 216L, 240L,
  264L, 288L, 312L, 120L, 144L, 168L, 192L, 216L, 240L, 72L, 96L,
  120L, 144L, 168L, 192L, 216L, 96L, 120L, 144L, 168L, 192L), G = c(15,
  25.5, 27, 28, 28.5, 39.5, 41.5, 43, 13, 21.5, 29.5, 30.5, 32.5,
  35, 13.5, 28, 32.5, 33.5, 35, 39.5, 42, 6.5, 30, 39.5, 57, 58.5
  )), .Names = c("Temp", "Time", "G"), class = "data.frame", row.names = c(NA,
  -26L))

Data$GermRate <- 1/Data$Time
Data$Probits <- qnorm(p=Data$G/100) # Get Probits


fm1 <-
  nls(
      formula= Probits ~ K+1/Sigma*(Temp-T0)*Time
    , data=Data
    , start=list(K=1, Sigma=2, T0=2)
    #, algorithm= "port"
    )
fm1Summary <- summary(fm1)
fm1Coef <- summary(fm1)$coef

- MYaseen208

---

如果您想要一个logit链接，那么您的响应值应该在0-1之间，但是您的“Probits”值包含负值。nls没有链接函数属性，您只需要在模型公式中包含转换即可（假设您使用了有意义的转换）。或者也许我误解了您在这里想要使用链接函数的意思。您能否更明确地指定您想要拟合的模型？ - MrFlick

感谢 @MrFlick 对我的问题表示关注。 实际上，响应变量是 Data $ G / 100，我从中创建了 Probits。 作为统计学的学生，我知道可以通过使用 linking 来避免这种 transformation 。 - MYaseen208

只有当拟合算法允许时，链接才起作用。nls（非线性最小二乘）通常不允许这样做。链接通常与glm相关联。在R中，您模型的响应是公式中~左侧的内容。您是想要进行probit还是logit链接？听起来您可能需要更多的统计建议而不是编程建议。如果是这样，您可以考虑发布到[stats.se]。 - MrFlick

你应该能够通过在 RHS 上应用反向链接函数来进行类似链接的等价操作，例如 GermRate ~ pnorm(K+1/Sigma*(Temp-T0)*Time) -- 第一次尝试没有成功，但这是我会采取的方向。 - Ben Bolker

我认为你应该在公式的右侧加上I，就像这样：Probits ~ I(K+1/Sigma*(Temp-T0)*Time)。否则，它将尝试使用公式语法解析它（将*解释为交互作用），而不是使用正常的数学语法。 - shadowtalker

另外，您应该发布出现了什么问题。您是否遇到错误？您的结果是否毫无意义？ - shadowtalker

Answer 1

您可以使用广义非线性模型的gnm包来适配这种类型的模型。这需要一些工作，因为gnm使用类"nonlin"的预定义函数来指定模型中的非线性项，而该软件包提供的函数通常不足以指定任意非线性函数。但是，可以定义自定义的"nonlin"函数与gnm一起使用。

在您的模型中，k是一个线性参数，因此我们只需要担心第二个术语。可以通过以下"nonlin"函数来指定它：

customNonlin <- function(Temp, Time){
    list(predictors = list(sigma = 1, t0 = 1),
         variables = list(substitute(Temp), substitute(Time)),
         term = function(predLabels, varLabels) {
             sprintf("1/%s * (%s - %s) * %s",
                     predLabels[1], varLabels[1], 
                     predLabels[2], varLabels[2])
         })
}
class(customNonlin) <- "nonlin"

在返回的列表中，

• predictors 指定 sigma 和 t0 是具有单个截距项的预测变量（即为独立参数）。
• variables 指定有两个变量，将通过 Temp 和 Time 参数由用户提供。
• term 指定一个函数，用于创建给定预测变量和变量标签的解析数学表达式。

关于 "nonlin" 函数的更多细节可以在 gnm vignette 的第 3.5 节中找到。

现在我们可以尝试按以下方式拟合您的模型。

mod1 <- gnm(cbind(G, 100 - G) ~ customNonlin(Temp, Time), family = binomial, 
            data  = Data, start = c(1, 2, 2))

请注意，与glm一样，默认情况下会向公式中添加截距项，这里将代表k。虽然起始值与解决方案相差甚远，但在此时满足gnm的收敛准则，因此算法不执行任何迭代。在这种情况下，需要更好的起始sigma估计值，以便gnm收敛到更合理的结果。

mod2 <- gnm(cbind(G, 100 - G) ~ customNonlin(Temp, Time), family = binomial, 
            data  = Data, start = c(1, 2000, 2))

mod2

Call:
gnm(formula = cbind(G, 100 - G) ~ customNonlin(Temp, Time), family = binomial, 
    data = Data, start = c(1, 2000, 2))

Coefficients:
(Intercept)        sigma           t0  
     -2.589     1915.602        8.815  

Deviance:            53.53157 
Pearson chi-squared: 49.91347 
Residual df:         23 

实际上，使用由 gnm 提供的 Mult 函数来指定此模型是可以的，只要您不介意重新参数化模型：

mod3 <- gnm(cbind(G, 100 - G) ~ Mult(1, 1 + offset(Temp), offset(Time)), 
            family = binomial, data  = Data,
            start = c(1, 1/2000, -2))

mod3

Call:

gnm(formula = cbind(G, 100 - G) ~ Mult(1, offset(Temp) + 1, offset(Time)), 
    family = binomial, data = Data, start = c(1, 1/2000, -2))

Coefficients:
                             (Intercept)  
                               -2.588874  
Mult(., 1 + offset(Temp), offset(Time)).  
                                0.000522  
Mult(1, . + offset(Temp), offset(Time)).  
                               -8.815152  

Deviance:            53.53157 
Pearson chi-squared: 49.91347 
Residual df:         23 

编辑

参数的功能在由customNonlin返回的列表的term组件中指定，您可以通过以下方式查看：

customNonlin(Temp, Time)$term(c("sigma", "t0"), c("Temp", "Time"))
"1/sigma * (Temp - t0) * Time"

如果你只想改变功能形式，你需要修改term函数。如果你想添加/删除参数，你也需要修改predictors组件中的列表。同样地，如果新术语需要你添加/删除变量，你将修改variables组件。

- Heather Turner

---

感谢@HeatherTurner的回答。我想知道我们在哪里指定了模型项。不知道如果公式改变，那么customNonlin函数会怎样改变？ - MYaseen208

我已经使用 sprintf 来创建字符串，每个 %s 都被调用 sprintf 的参数替换。更多详情请查看 ?sprintf。你也可以使用 paste 来创建字符串，不一定非要用 sprintf。 - Heather Turner

对于 nonlin 函数的调用只是模型中的一个项 - 因此您可以像在 glm 中一样添加到公式中。具体来说，要删除截距，您可以将 0 或 -1 添加到模型公式中。 - Heather Turner

请问，何时使用gnm，何时使用glm或nls？比如说我想拟合一个逻辑回归模型。 - skan

glm用于使用指数族分布（例如高斯，二项式，泊松）对响应进行建模，并将期望值链接到线性预测器（其中每个项的形式为coef * variable）。nls用于对具有非线性预测器的连续无界响应进行建模（项是参数的非线性函数）。gnm将这两种方法结合起来，即类似于glm，但允许预测器中的非线性项。逻辑回归使用logit链接对二项式响应进行建模-因此根据预测器使用glm或gnm，还可以参见gnm演示文稿的第2节以获取glm增强功能。 - Heather Turner