R中拟合广义非线性模型

Question

我想拟合以下广义非线性模型:Probit(G)=K+1/Sigma*(Temp-T0)*Time.作为天真的模特,我创造Probits(G)了qnorm(G)然后适合了Nonlinear Model.但我想,以适应Nonlinear Model与logit链接类似glm的功能R.如果有人帮我配合带logit链接的广义非线性模型,我将非常感激R.在此先感谢您的帮助.

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插图的第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 

其实这是可以使用指定这个模型Mult所提供的功能gnm,只要你不介意重新参数化模型:

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 

编辑

参数的功能在term返回的列表的组件中指定customNonlin,您可以通过该组件查看

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

因此,如果您只想更改功能表单,则需要修改该term功能.如果要添加/删除参数,还需要修改predictors组件中的列表.同样,如果新术语要求您添加/删除变量,则可以修改该variables组件.