估计断点直到收敛

Question

我有以下数据：

> str(sample)
'data.frame':   500 obs. of  3 variables:
 $ X        : int  1 2 3 4 5 6 7 8 9 10 ...
 $ predictor: num  0.446 0.395 0.484 0.919 0.844 ...
 $ response : num  6.65 6.22 6.54 7.85 7.34 ...

我需要估计逐步回归中使用的断点，我正在遵循 Muggeo (2003) 的算法。这是预测变量与响应的关系图：

我创建了以下函数：

breakpoint<-function(psi,data){
repeat{
  psi<-psi
  predictor<-data[,1]
  response<-data[,2]
  
  U<-ifelse(predictor>=psi,predictor-psi,0)
  V<-ifelse(predictor>=psi,-1,0)
  
  model<-lm(response~predictor+U+V,data=data)
  psi_new<-(model$coefficients[4]/model$coefficients[3])+psi
  psi<-psi_new
  if (abs(psi-psi_new)<0.00001 && model$coefficients[4]<0.00001) break
}
  result<-c(psi,model$coefficients[4])
  return(result)
}

但是这个函数对psi的初始值非常敏感，我不确定为什么会这样。我认为这可能与我的收敛标准有关。当 psi 和 psi_new 之间的“差距”可以忽略不计并且 V 的系数基本为零时，函数应该收敛。

包分段确定的正确断点是0.27，当我输入0.25的psi时我得到正确的值，但是当我输入0.4时，估计值变为0.32。

Answer 1

也许是这样的：

## geneate some test data
set.seed(132)
predictor <- seq(0, 0.999, length.out = 500)
response  <- approx(c(0, 0.27, 1), c(6.5, 6, 8.5), predictor)$y + rnorm(predictor, sd = 0.2)
sample <- data.frame(predictor, response)

breakpoint <- function(psi, data, tol=1e-5, maxit=50){
  iterations <- 0
  predictor  <- data[,1]

  repeat{
    U <- ifelse(predictor >= psi, predictor - psi, 0)
    V <- ifelse(predictor >= psi,              -1, 0)

    model <- lm(response ~ predictor + U + V, data = data)
    p <- coef(model)
    psi_new <- (p[4]/p[3]) + psi
    psi <- psi_new
    if ((abs(psi - psi_new) < tol) & abs(p[4]) < tol) break
    if (iterations > maxit) {
      warning("too many iterations")
      break
    }
    iterations <- iterations + 1
  }
  result <- c(psi = unname(psi), p[4])
  return(result)
}

psi <- 0.1
breakpoint(psi, sample, tol=1e-5) # warning, but psi is ok
breakpoint(psi, sample, tol=1e-2) # no warning

plot(predictor, response, cex=0.5, pch=16)
abline(v=breakpoint(psi, sample, tol=1e-2)["psi"], lty="dashed", col="red")

## for comparison
library("segmented")
segmented(lm(response ~ predictor, data=sample))

在这个例子中，断点被很好地识别，但优化在理论上可以陷入局部最小值。因此，尝试不同的起始值可能是明智的。

编辑：

再次编辑，使用abs(p4)
避免无限循环的最大迭代次数

估计断点直到收敛

Estimate breakpoint until convergence

estimation

regression

r

breakpoints

repeat