楼主: 我的素质低
2226 1

[R] 〖移花接木〗分段回归的拐点连续性 [推广有奖]

学术权威

83%

还不是VIP/贵宾

-

TA的文库  其他...

〖素质文库〗

结构方程模型

考研资料库

威望
8
论坛币
23388 个
通用积分
28302.5637
学术水平
2705 点
热心指数
2881 点
信用等级
2398 点
经验
224313 点
帖子
2977
精华
52
在线时间
2175 小时
注册时间
2012-11-24
最后登录
2024-1-13

一级伯乐勋章 初级学术勋章 初级热心勋章 初级信用勋章 中级热心勋章 中级学术勋章 中级信用勋章 高级学术勋章 高级热心勋章 高级信用勋章 特级学术勋章

+2 论坛币
k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群
赵安豆老师微信:zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

感谢您参与论坛问题回答

经管之家送您两个论坛币!

+2 论坛币

——来自谢益辉老师




      我翻了一下,上一次写跟统计有关的话题已经是2010年的事情了,多数时候,这里写的都是(计算机)技术类与生活类的话题,于是我又想起来很久以前有客官说,你还不如改成“Keep on eating”算了。废话少说,言归正传。





       最近有客官问:

     在分段线性回归中,转折点的连续性是必须的吗?

     意思是这样,下图中,两段直线应该接着呢(实线),还是可以各走各的(虚线)?

     思是这样,下图中,两段直线应该接着呢(实线),还是可以各走各的(虚线)?


  1. set.seed(123)
  2. # 真实模型
  3. x = sort(runif(100))
  4. y = 2 + 1 * x + 4 * (x > 0.5) + 3 * (x - 0.5) * (x > 0.5) + rnorm(100)
  5. par(mar = c(4, 4, 0, 0), family = "serif", mgp = c(2, 1, 0))
  6. plot(x, y, pch = 20, col = "darkgray")
  7. # 斜率不同,截距限制
  8. fit1 = lm(y ~ 1 + x + I((x - 0.5) * (x > 0.5)))
  9. lines(x, fitted(fit1), lwd = 2, col = 2)
  10. # 斜率不同,截距也不同
  11. fit2 = lm(y ~ 1 + x + I(x > 0.5) + I((x - 0.5) * (x > 0.5)))
  12. lines(x, fitted(fit2), lwd = 2, lty = 2)
复制代码








       虚线线段的中间一段不应该被画出来,我偷了个懒画上去了,本来应该是间断的。连与不连,归根结底,只不过是回归设计阵要不要多一列的问题。这是线性回归中经典的模型比较问题,上面fit2是全模型,fit1是取消一个限制的模型,自由度之差为1,残差平方和之差除以1,再除以全模型的残差平方和除以相应自由度(88/96),就是F统计量,P值也就出来了。R里面所谓的anova()函数,就是干这事儿的:



  1. # 两个嵌套模型做F检验
  2. anova(fit1, fit2)
  3. ## Analysis of Variance Table
  4. ##
  5. ## Model 1: y ~ 1 + x + I((x - 0.5) * (x > 0.5))
  6. ## Model 2: y ~ 1 + x + I(x > 0.5) + I((x - 0.5) * (x > 0.5))
  7. ## Res.Df RSS Df Sum of Sq F Pr(>F)
  8. ## 1 97 155
  9. ## 2 96 88 1 66.6 72.6 2.2e-13 ***
  10. ## ---
  11. ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
复制代码


      这里因为我们设计的真实模型就是两部分直线斜率和截距在0.5前后都不同,所以这里F检验高度著,说明加上截距限制条件不合理。道理很简单:除非简约模型和全模型没有显著差异,我们只能认为全模型相对“好”一些。




      然而,线性模型本身就是对现实的简化。按照上面的想法,其实我们还需要看更复杂的全模型,否则不能认为线性模型就是合适的模型,比如可以考虑非线性模型。复杂模型是没有尽头的,如果要一个个检查下去,那大家都不用干活了。2008年我在COS写了一篇简单的LOWESS文章,如今我仍然可以把它搬出来,因为LOWESS实在是回归模型中的战斗机,它是模型复杂与简约的很好平衡。前面的直线回归可以看作LOWESS的特例,如下fit4。我们可以构造一个相对复杂的LOWESS模型(span参数取小一些),然后和一个简单的模型比较,如:



.
  1. fit3 = loess(y ~ x, span = 0.2)
  2. fit4 = loess(y ~ 1 + x + I(x > 0.5) + I((x - 0.5) * (x > 0.5)), span = 1,
  3. degree = 1)
  4. par(mar = c(4, 4, 0, 0), family = "serif", mgp = c(2, 1, 0))
  5. plot(x, y, pch = 20, col = "darkgray")
  6. lines(x, fitted(fit3), lwd = 2, col = 2)
  7. lines(x, fitted(fit4), lwd = 2, lty = 2)
复制代码






  1. anova(fit3, fit4)
  2. ## Model 1: loess(formula = y ~ x, span = 0.2)
  3. ## Model 2: loess(formula = y ~ 1 + x + I(x > 0.5) + I((x - 0.5) * (x > 0.5)), span = 1, degree = 1)
  4. ##
  5. ## Analysis of Variance:   denominator df 82.38
  6. ## ##        ENP  RSS F-value Pr(>F)## [1,] 15.16 83.2               
  7. ## [2,]  4.02 87.7   0.321   0.99
复制代码






      P值很大,这毫不奇怪,因为真实模型就是按照两段直线构造的。弯弯曲曲的复杂模型无法打败简单模型,此时我们可以有点底气说两个线段的直线回归可能是对数据的一个恰当描述。








客官的第二个问题是:

转折点怎么确定?




      好问题。我对学术研究其实没多大兴趣,所以懒得去文献堆里翻找。我总觉得很多研究是把问题复杂化了,这种文献海让我觉得窒息。如果是我,第一步肯定是画平滑曲线,然后再去结合具体事件背景去决定。如前面的数据用ggplot2很容易分组画图:




  1. library(ggplot2)qplot(x, y) + geom_smooth()  # 总趋势
复制代码


  1. qplot(x, y, group = x > 0.5) + geom_smooth()  # 按0.5前后分组
复制代码




      显然所有数据放一条线上建模不合适。作为普通青年的我,到此就止步了(所以我离学术还十万八千里)。不过作为二青年想起一件事,如果不知道转折点,而靠不断推移转折点、新建模型、比较模型得来的转折点,是否有多重比较的嫌疑?比如先试0.2为转折点,再试0.3,……,一直试到0.9,看哪个模型跟全模型比起来没有显著差异。或者说,如果本来真实模型中不存在这样一个转折点,而通过这种方法找出转折点的概率有多大?我估计我八成问了一个很二的问题,应该早有人研究过了。做个模拟应该也不难。









library(knitr)


render_jekyll()

opts_chunk$set(cache = TRUE)

opts_knit$set(upload.fun = imgur_upload)



local({

  all.md = list.files(pattern = '^_.*\\.Rmd$')



  for (f in all.md) {

    b = gsub('^_|\\.Rmd$', '', f)

    unlink(paste0(b, '.html'))

    unlink(list.files('_figure/', '^(^\\d)+', full.names = TRUE))

    opts_chunk$set(cache.path = paste('_cache', b, sep = '/'),

                   fig.path = paste('_figure', b, sep = '/'))

    message('processing ', f)

    knit(f, output = paste0(b, '.md'))

  }

})




二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

关键词:移花接木 连续性 Processing variance Analysis 计算机 family eating 转折点 谢益辉

已有 1 人评分经验 论坛币 学术水平 热心指数 信用等级 收起 理由
niuniuyiwan + 60 + 60 + 5 + 5 + 5 精彩帖子

总评分: 经验 + 60  论坛币 + 60  学术水平 + 5  热心指数 + 5  信用等级 + 5   查看全部评分

本帖被以下文库推荐

心晴的时候,雨也是晴;心雨的时候,晴也是雨!
扣扣:407117636,欢迎一块儿吐槽!!
沙发
xddlovejiao1314 学生认证  发表于 2015-3-27 14:00:26 |只看作者 |坛友微信交流群
好贴,赞一个。谢谢楼主分享。有空我也做一个Stata做分段回归的列子。
已有 1 人评分论坛币 热心指数 收起 理由
niuniuyiwan + 10 + 2 精彩帖子

总评分: 论坛币 + 10  热心指数 + 2   查看全部评分

使用道具

您需要登录后才可以回帖 登录 | 我要注册

本版微信群
加好友,备注jltj
拉您入交流群

京ICP备16021002-2号 京B2-20170662号 京公网安备 11010802022788号 论坛法律顾问:王进律师 知识产权保护声明   免责及隐私声明

GMT+8, 2024-6-16 23:46