title: changepoint detection in R
body: <p>Anybody an idea how to solve this R issue? It is to find a changepoint in a relation, like x=5 in data below.</p>
  
  <pre><code>fitDat &lt;- data.frame(x=1:10, y=c(rnorm(5, sd=0.2),(1:5)+rnorm(5, sd=0.2)))
  plot(fitDat$x, fitDat$y)
  </code></pre>
  
  <p>SAS works beautifully well;</p>
  
  <pre><code>/*simulated example*/
  data test;
  INPUT id $ x y;
  cards;
  1   1 -0.22711769
  2   2 -0.08712909
  3   3  0.06922072
  4   4 -0.12940913
  5   5 -0.43152927
  6   6  1.17685016
  7   7  1.83410448
  8   8  2.88528795
  9   9  4.30078012
  10 10  4.84517101
  ;
  run;
  
  proc nlmixed data=test;
      parms B0 = 0 B1=1 t0=5 sigma_r=1;
      mu1 = (x &lt;= t0)*B0;
      mu2 = (x &gt; t0)*(B0 + B1*(x - t0));
      mu=mu1+mu2;
      model y~normal(mu,sigma_r);
  run;
  </code></pre>
  
  <p>R fails;</p>
  
  <pre><code>changePointOptim &lt;- function(x, int, slope, delta){ # code adapted from https://stats.stackexchange.com/questions/7527/change-point-analysis-using-rs-nls
         int + (pmax(delta-x, 0) * slope)
  }
  
  # nls
  nls(y ~ changePointOptim(x, b0, b1, delta), 
                data = fitDat, 
                start = c(b0 = 0, b1 = 1, delta = 5))
  
  # optimization
  sqerror &lt;- function (par, x, y)
         sum((y - changePointOptim(x, par[1], par[2], par[3]))^2)
  minObj &lt;- optim(par = c(0, 1, 3), fn = sqerror,
                x = fitDat$x, y = fitDat$y, method = "BFGS")
  minObj$par
  
  # nlme
  library(nlme)
  nlmeChgpt &lt;- nlme(y ~ changePointOptim(b0, b1, delta,x), 
                data = fitDat, 
                fixed = b0 + b1 + delta, 
                start = c(b0=0, b1=1, delta=5)) 
  summary(nlmeChgpt)
  </code></pre>
  
  <p>As it is a simple case with a clear signal, R should work. I wonder what I’m doing wrong in R (I used some code from <a href="https://stats.stackexchange.com/questions/7527/change-point-analysis-using-rs-nls">https://stats.stackexchange.com/questions/7527/change-point-analysis-using-rs-nls</a>). Anybody a suggestion/solution?</p>
  
  <p>Thanks!</p>
  
  <p>Willem</p>