r - 用分类变量拟合 nls 模型

Question

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

这是我目前所拥有的:

# Create data
df1 <- cbind.data.frame (height = c (0.5, 0.6, 0.9, 1.3, 1.5, 1.6, 1.6,
                                     0.6, 0.6, 0.8, 1.3, 1.5, 1.6, 1.5,
                                     0.6, 0.8, 1.0, 1.4, 1.6, 1.6, 1.6,
                                     0.5, 0.8, 1.0, 1.3, 1.6, 1.7, 1.6),
                         age = c (0.5, 0.9, 3.0, 7.3, 12.2, 15.5, 20.0,
                                  0.4, 0.8, 2.3, 8.5, 11.5, 14.8, 21.3,
                                  0.5, 1.0, 5.1, 11.1, 12.3, 16.0, 19.8,
                                  0.5, 1.1, 5.5, 10.2, 12.2, 15.4, 20.5),
                         region = as.factor (c (rep ("A", 7),
                                                rep ("B", 7),
                                                rep ("C", 7),
                                                rep ("D", 7))))

> head (df1)
  height  age region
1    0.5  0.5      A
2    0.6  0.9      A
3    0.9  3.0      A
4    1.3  7.3      A
5    1.5 12.2      A
6    1.6 15.5      A

# Create linear-plateau function
lp <- function(x, a, b, c){
  ifelse (x < c, a + b * x, a + b * c)
  } # Where 'a' is the intercept, 'b' the slope and 'c' the breakpoint

# Fit the model ignoring region
m1 <- nls (height ~ lp (x = age, a, b, c),
           data = df1,
           start = list (a = 0.5, b = 0.1, c = 13))

> summary (m1)

Formula: height ~ lp(x = age, a, b, c)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
a  0.582632   0.025355   22.98   <2e-16 ***
b  0.079957   0.003569   22.40   <2e-16 ***
c 12.723995   0.511067   24.90   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07468 on 25 degrees of freedom

Number of iterations to convergence: 2 
Achieved convergence tolerance: 5.255e-09

我想拟合相同的模型，但考虑到 region，并测试 a、b 和 c > 不同地区的估计不同。

我相信this post可能有用，但我不知道如何将它应用于此数据/函数。

这是数据的样子:

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

Answer 1

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

我们对 fm1 使用 plinear 算法，因为它不需要线性参数的起始值。在这种情况下，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)