使用分类变量拟合 nls 模型

Question

我想要拟合一个线性平台 (nls) 模型，将身高描述为年龄的函数，并且我想测试区域之间模型的任何参数是否存在显着差异。

\n

这是我到目前为止所拥有的：

\n

# Create data\ndf1 <- cbind.data.frame (height = c (0.5, 0.6, 0.9, 1.3, 1.5, 1.6, 1.6,\n                                     0.6, 0.6, 0.8, 1.3, 1.5, 1.6, 1.5,\n                                     0.6, 0.8, 1.0, 1.4, 1.6, 1.6, 1.6,\n                                     0.5, 0.8, 1.0, 1.3, 1.6, 1.7, 1.6),\n                         age = c (0.5, 0.9, 3.0, 7.3, 12.2, 15.5, 20.0,\n                                  0.4, 0.8, 2.3, 8.5, 11.5, 14.8, 21.3,\n                                  0.5, 1.0, 5.1, 11.1, 12.3, 16.0, 19.8,\n                                  0.5, 1.1, 5.5, 10.2, 12.2, 15.4, 20.5),\n                         region = as.factor (c (rep ("A", 7),\n                                                rep ("B", 7),\n                                                rep ("C", 7),\n                                                rep ("D", 7))))\n\n> head (df1)\n  height  age region\n1    0.5  0.5      A\n2    0.6  0.9      A\n3    0.9  3.0      A\n4    1.3  7.3      A\n5    1.5 12.2      A\n6    1.6 15.5      A\n\n# Create linear-plateau function\nlp <- function(x, a, b, c){\n  ifelse (x < c, a + b * x, a + b * c)\n  } # Where \'a\' is the intercept, \'b\' the slope and \'c\' the breakpoint\n\n# Fit the model ignoring region\nm1 <- nls (height ~ lp (x = age, a, b, c),\n           data = df1,\n           start = list (a = 0.5, b = 0.1, c = 13))\n\n> summary (m1)\n\nFormula: height ~ lp(x = age, a, b, c)\n\nParameters:\n   Estimate Std. Error t value Pr(>|t|)    \na  0.582632   0.025355   22.98   <2e-16 ***\nb  0.079957   0.003569   22.40   <2e-16 ***\nc 12.723995   0.511067   24.90   <2e-16 ***\n---\nSignif. codes:  0 \xe2\x80\x98***\xe2\x80\x99 0.001 \xe2\x80\x98**\xe2\x80\x99 0.01 \xe2\x80\x98*\xe2\x80\x99 0.05 \xe2\x80\x98.\xe2\x80\x99 0.1 \xe2\x80\x98 \xe2\x80\x99 1\n\nResidual standard error: 0.07468 on 25 degrees of freedom\n\nNumber of iterations to convergence: 2 \nAchieved convergence tolerance: 5.255e-09\n\n

\n

我想拟合相同的模型，但考虑region、 并测试a、b和c估计值在区域之间是否不同。

\n

我相信这篇文章可能有用，但我不知道如何将其应用于此数据/函数。

\n

数据如下所示：

\n

\n

不使用 nls 的解决方案也受欢迎

\n

Answer 1

对每个区域使用相同的参数拟合模型，给出 fm1，并再次使用不同的参数拟合模型，给出 fm2，并使用方差分析来测试差异。

我们使用plinearfm1 算法，因为它不需要线性参数的起始值。在这种情况下，RHS 应该是一个矩阵，其第一列乘以截距，第二列乘以斜率。这两个线性参数将被命名为.lin1和.lin2。我们使用重复 4 次的 fm1 系数作为 fm2 拟合的起始值。

fm1 <- nls(height ~ cbind(1, pmin(age, c)), df1, start = list(c = mean(df1$age)),
  algorithm = "plinear")
co <- as.list(coef(fm1))

fm2 <- nls(height ~ a[region] + b[region] * pmin(age, c[region]), df1, 
  start = list(a = rep(co$.lin1, 4), b = rep(co$.lin2, 4), c = rep(co$c, 4)))

anova(fm1, fm2)

给予：

Analysis of Variance Table

Model 1: height ~ cbind(1, pmin(age, c))
Model 2: height ~ a[region] + b[region] * pmin(age, c[region])
  Res.Df Res.Sum Sq Df   Sum Sq F value Pr(>F)
1     25    0.13944                           
2     16    0.11895  9 0.020483  0.3061 0.9617

因此我们不能拒绝跨区域参数相同的假设。

如果我们希望测试不同的 c 值但常见的截距和斜率，我们可以使用

fm3 <- nls(height ~ cbind(1, pmin(age, c[region])), df1, 
   start = list(c = rep(co$c, 4)), algorithm = "plinear")

anova(fm1, fm3)

虽然我们不能拒绝这样的假设，即 c 的值在视觉上在各个区域是相同的，但我们看到平台值的截止年龄确实看起来有些不同，因此我们可能想要使用 fm3，即使它与 fm1 没有显着差异。我们可能希望以与此处的应用程序相关的其他因素为指导，而不仅仅是适合性。

图形

下面我们展示了 fm2 的单独拟合和 fm1 的整体拟合。

library(ggplot2)

df1$Everything <- "Everything"
ggplot(df1, aes(age, fitted(fm2), col = region)) +
  geom_line() +
  geom_point() +
  geom_line(aes(age, fitted(fm1), col = Everything), lty = 2, lwd = 2)