我想要拟合以下广义非线性模型: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
最佳答案
您可以使用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
组件。
关于r - 在 R 中拟合广义非线性模型,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/26344429/