我们需要更好的起始值。拟合一个非高原模型，model0，并使用其中的参数来拟合给出模型的所有数据点，然后使用其中的 a 和 b 以及 clx 的值网格（由于其有问题的性质）给出 model.Ab 和模型.La。（请注意，它无法根据网格的某些起始值生成拟合，从而导致错误消息，但 nls2 将继续处理进一步的起始值，因此可以忽略这些错误。）

library(nls2)

# ensure data is sorted for plotting
o <- with(healing, order(Type, Days))
h <- healing[o, ]

# last argument specifies whether there is or is not a plateau
quadplat = function(x, a, b, clx, plat = TRUE) {
  if (plat) x <- pmin(x, clx)
  a + b * x + (-0.5*b/clx) * x * x
}

# fit no plateau model with all data
st <- c(a = 1, b = 1, clx = 1)
model0 <- nls(Area ~ quadplat(Days, a, b, clx, FALSE), h, start = st)

# fit all data model
model <- nls(Area ~ quadplat(Days, a, b, clx), h, start = coef(model0))
co <- coef(model)

我们现在可以使用上面在起始值中计算的值来拟合和绘制子集模型。

if (exists("model.Ab")) rm(model.Ab)
model.Ab <- nls2(Area ~ quadplat(Days, a, b, clx), h, subset = h$Type == "Abrasion",
  start = data.frame(a = co[[1]], b = co[[2]], clx = 0:140))

if (exists("model.La")) rm(model.La)
model.La <- nls2(Area ~ quadplat(Days, a, b, clx), h, subset = h$Type == "Laceration",
  start = data.frame(a = co[[1]], b = co[[2]], clx = 0:140))
  
cols <- c(Abrasion = "red", Laceration = "blue")
plot(Area ~ Days, h, col = cols[Type], pch = 20, cex = 1.5)
lines(fitted(model.Ab) ~ Days, subset(h, Type == "Abrasion"), 
  col = cols["Abrasion"])
lines(fitted(model.La) ~ Days, subset(h, Type == "Laceration"), 
  col = cols["Laceration"])

（图后续）

替代模型

如果可以考虑其他模型，则该模型只有两个参数，更容易拟合，尽管参数较少，但残差平方和也较低。

model.Ab2 <- nls(Area ~ a * (1 - exp(- b * Days)), h, 
   subset = Type == "Abrasion", start = c(a = 100, b = .1))
 
model.La2 <- nls(Area ~ a * (1 - exp(- b * Days)), h, 
   subset = Type == "Laceration", start = c(a = 100, b = .1))

# plot
cols <- c(Abrasion = "red", Laceration = "blue")
plot(Area ~ Days, h, col = cols[Type], pch = 20, cex = 1.5)
lines(fitted(model.Ab2) ~ Days, subset(h, Type == "Abrasion"), 
  col = cols["Abrasion"])
lines(fitted(model.La2) ~ Days, subset(h, Type == "Laceration"), 
  col = cols["Laceration"])

（图后续）

一参数模型

如果我们在上一节的 2 参数模型中固定 a = 100，我们会得到一个 1 参数模型，它在统计上与 2 参数模型没有区别。从方差分析中显示的大于 0.05 的 p 值可以看出这一点，这表明我们不能拒绝零假设，即 1 和 2 参数模型对两个子集中的每一个都同样地描述了数据。

model.Ab3 <- nls(Area ~ 100 * (1 - exp(- b * Days)), h, 
   subset = Type == "Abrasion", start = c(b = .1))
 
model.La3 <- nls(Area ~ 100 * (1 - exp(- b * Days)), h, 
   subset = Type == "Laceration", start = c(b = .1))

anova(model.Ab3, model.Ab2)

anova(model.La3, model.La2)

还要注意它达到 y = 95 的点，即接近高原，是-log(1 - 95/100)/b（基于反转模型方程）。分子大约是 3，所以它大约在 0 处达到 95 3/b。

其他

如果m <- nls(...)thensummary(m)将给出系数和其他信息的标准误差。

另外，如果有人擅长解释公式，你能帮我把这段代码写成一个可读的公式吗？

function(x, a, b, clx) {
  ifelse(x  < clx, a + b * x + (-0.5*b/clx) * x * x, 
         a + b * clx + (-0.5*b/clx) * clx * clx)}

$$f(x, a, b, x_{cl}) = \begin{cases} a + bx + (\frac{-0.5b}{x_{cl}}) \times x^2 , & \text{if}\ x < x_{cl} \\ a + bx_{cl} + (\frac{-0.5b}{x_{cl}}) \times {x_{cl}}^2 , & \text{otherwise} \end{cases}$$

这简化为：

$$f(x, a, b, x_{cl}) = \begin{cases} a + bx \left( 1 - \frac{x}{2x_{cl}} \right) , & \text{if}\ x < x_{cl} \\ a + \frac{bx_{cl}}{2} , & \text{otherwise} \end{cases}$$

我替换的地方$x_{cl}$为了clx使其更具可读性。