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

Question

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

如何将这种广义非线性模型与 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

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 组件。