3

我想拟合以下广义非线性模型：Probit(G)=K+1/Sigma*(Temp-T0)*Time. 作为天真的模型，我创建Probits(G)并qnorm(G)安装了Nonlinear Model. 但我想适应类似于功能的Nonlinear Model链接。logitglmR

我怎样才能用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
4

1 回答 1

9

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

实际上可以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组件。

于 2014-10-16T14:35:21.367 回答