Stack Overflow 上有一个非常相似的帖子:使用二次多项式和在断点处平滑连接的直线进行分段回归 https://stackoverflow.com/q/39418783/4891738。唯一的区别是我们现在考虑:
事实证明,函数est
, choose.c
and pred
定义于我的答案 https://stackoverflow.com/a/40817449/4891738根本不需要改变;我们只需要修改getX
返回分段回归的设计矩阵:
getX <- function (x, c) cbind("beta0" = 1, "beta1" = pmin(x - c, 0))
现在,我们按照以下代码进行操作玩具示例 https://stackoverflow.com/a/40817494/4891738使模型适合您的数据:
x <- c(1, 2, 3, 1, 2, 1, 6, 1, 2, 3, 2, 1, 4, 3, 1)
y <- c(0.041754212, 0.083491254, 0.193129615, 0.104249201, 0.17280516,
0.154342335, 0.303370501, 0.025503008, 0.123934121, 0.191486527,
0.183958737, 0.156707866, 0.31019215, 0.281890206, 0.25414608)
x
范围从1到6,所以我们考虑
c.grid <- seq(1.1, 5.9, 0.05)
fit <- choose.c(x, y, c.grid)
fit$c
# 4.5
最后我们做出预测图:
x.new <- seq(1, 6, by = 0.1)
p <- pred(fit, x.new)
plot(x, y, ylim = c(0, 0.4))
matlines(x.new, p[,-2], col = c(1,2,2), lty = c(1,2,2), lwd = 2)
我们在拟合模型中有丰富的信息:
str(fit)
#List of 12
# $ coefficients : num [1:2] 0.304 0.055
# $ residuals : num [1:15] -0.06981 -0.08307 -0.02844 -0.00731 0.00624 ...
# $ fitted.values: num [1:15] 0.112 0.167 0.222 0.112 0.167 ...
# $ R : num [1:2, 1:2] -3.873 0.258 9.295 -4.37
# $ sig2 : num 0.00401
# $ coef.table : num [1:2, 1:4] 0.3041 0.055 0.0384 0.0145 7.917 ...
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:2] "beta0" "beta1"
# .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
# $ aic : num -34.2
# $ bic : num -39.5
# $ c : num 4.5
# $ RSS : num 0.0521
# $ r.squared : num 0.526
# $ adj.r.squared: num 0.49
例如,我们可以检查系数汇总表:
fit$coef.table
# Estimate Std. Error t value Pr(>|t|)
#beta0 0.30406634 0.03840657 7.917039 2.506043e-06
#beta1 0.05500095 0.01448188 3.797915 2.216095e-03