<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/rupture-detection/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://freakonometrics.hypotheses.org/49848"> R-english – Freakonometrics</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>There are some graphs that you cannot forget. One graph that I found puzzling was mentioned on Andrew Gelman’s blog, <a target="_blank" rel="nofollow" href="http://andrewgelman.com/2013/08/05/evidence-on-the-impact-of-sustained-use-of-polynomial-regression-on-causal-inference-a-claim-that-coal-heating-is-reducing-lifespan-by-5-years-for-half-a-billion-people/">a few years back</a>, and was related to rupture detection</p>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2016/12/Screen-Shot-2013-08-03-at-4.23.29-PM-1024x874.png?resize=456%2C389&ssl=1" data-recalc-dims="1"></p>
<p>What I remember from this graph is that if you want to get a rupture, you can easily find one…</p>
<p><span></span></p>
<p>Recently, I had to review a paper, and Imbens &amp; Lemieux (2008) was mentioned (the paper on <a target="_blank" rel="nofollow" href="http://www.nber.org/papers/w13039">Regression Discontinuity</a>). It is a great reference and a major contribution, undoublty. But now, each time I read about ruptures and discontinuities, I have this graph in mind, and I start having doubt. I there really a discontinuity ? or is it somehow artificial, simply because the author is looking for a discontinuity ? For instance, consider the following (simulated) dataset</p>
<p><code>f=function(x) pnorm(x,mean=3)-x/4<br>
n=100<br>
set.seed(1)<br>
X=runif(n)*6<br>
Y=f(X)+rnorm(n)/8<br>
DB=data.frame(X,Y)</code></p>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2016/12/rupture2.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>It is a smooth continuous model. To test for discontinuity, a standard procedure seems to fit models on the left, and on the right. For instance two linear regression models,</p>
<p><code>s=2<br>
i=1<br>
reg1=lm(Y~poly(X,i),data=DB,subset=X&lt;=s)<br>
reg2=lm(Y~poly(X,i),daia=DB,subset=X>=s)<br>
u=seq(0,6,by=.025)<br>
v1=predict(reg1,newdata=data.frame(X=u[u&lt;=s]))<br>
v2=predict(reg2,newdata=data.frame(X=u[u>=s]))<br>
lines(u[u&lt;=s],v1,lwd=2)<br>
lines(u[u>=s],v2,lwd=2)</code></p>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2016/12/rupt12.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>or two quadratic ones,</p>
<p><code>i=2<br>
reg1=lm(Y~poly(X,i),data=DB,subset=X&lt;=s)<br>
reg2=lm(Y~poly(X,i),daia=DB,subset=X>=s)<br>
u=seq(0,6,by=.025)<br>
v1=predict(reg1,newdata=data.frame(X=u[u&lt;=s]))<br>
v2=predict(reg2,newdata=data.frame(X=u[u>=s]))<br>
lines(u[u&lt;=s],v1,lwd=2)<br>
lines(u[u>=s],v2,lwd=2)</code></p>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2016/12/rupt11.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>We observe here a discontinuity… but is it significant ? Consider here simulations of the dataset, and try several breakpoints</p>
<p><code>SIMU=function(ns=1e3){<br>
MAT=matrix(NA,ns,21)<br>
n=100<br>
for(j in 1:ns){<br>
X=runif(n)*6<br>
Y=f(X)+rnorm(n)/8<br>
DB=data.frame(X,Y)<br>
saut=function(s){<br>
reg1=lm(Y~poly(X,2),data=DB,subset=X&lt;=s)<br>
reg2=lm(Y~poly(X,2),data=DB,subset=X>=s)<br>
predict(reg2,newdata=data.frame(X=s))-predict(reg1,newdata=data.frame(X=s))<br>
}<br>
S=seq(.5,5.5,by=.25)<br>
MAT[j,]=Vectorize(saut)(S)<br>
}<br>
return(MAT)}</code></p>
<p>If we plot those differences, between the model on the left and the model on the right of the breakpoint, we get</p>
<p><code>SM=SIMU(1e3)<br>
S=seq(.5,5.5,by=.25)<br>
plot(S,SM[1,],type=&quot;l&quot;,ylim=c(-.7,.5),col=&quot;light blue&quot;)<br>
for(j in 2:100) lines(S,SM[j,],col=&quot;light blue&quot;)<br>
abline(h=0,lty=2)<br>
lines(S,apply(SM,2,mean),col=&quot;red&quot;,lwd=2)<br>
lines(S,apply(SM,2,function(x) quantile(x,.9)),col=&quot;red&quot;)<br>
lines(S,apply(SM,2,function(x) quantile(x,.1)),col=&quot;red&quot;)</code></p>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2016/12/Rupture1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>for various breakpoints. The good thing is that it looks like we should reject the idea of having a breakpoint here, since it might not be “significant”. But what if the change of concavity was sharper, in our simulated dataset</p>
<p><code>f=function(x) pnorm(x,mean=3)-x/4<br>
n=100<br>
set.seed(1)<br>
X=runif(n)*6<br>
Y=f(X)+rnorm(n)/8<br>
DB=data.frame(X,Y)</code></p>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2016/12/rupture3.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>With a simular code, we might now accept the idea of having a discontinuity, that we do not have actually…</p>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2016/12/rupture4.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>But more puzzeling, if we get back to the initial statement in <a target="_blank" rel="nofollow" href="http://andrewgelman.com/2013/08/05/evidence-on-the-impact-of-sustained-use-of-polynomial-regression-on-causal-inference-a-claim-that-coal-heating-is-reducing-lifespan-by-5-years-for-half-a-billion-people/">the post </a>of Andrew Gelman, it might actually be possible to chose the degrees of the polynomial regressions to validate the discontinuity hypothesis. Let us get back on our original dataset.</p>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2016/12/rupture2.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p><code>SIMU=function(ns){<br>
MAT=matrix(NA,ns,21)<br>
n=100<br>
for(j in 1:ns){<br>
X=runif(n)*6<br>
Y=f(X)+rnorm(n)/8<br>
DB=data.frame(X,Y)<br>
saut=function(s){<br>
reg1=reg2=rep(NA,3)<br>
for(i in 1:3){<br>
reg1[i]=predict(lm(Y~poly(X,i),data=DB,subset=X&lt;=s),newdata=data.frame(X=s))<br>
reg2[i]=predict(lm(Y~poly(X,i),data=DB,subset=X>=s),newdata=data.frame(X=s))<br>
}<br>
max(abs(rep(reg1,each=3)-rep(reg2,3)))}<br>
S=seq(.5,5.5,by=.25)<br>
MAT[j,]=Vectorize(saut)(S)<br>
}<br>
return(MAT)}</code></p>
<p>Here we try a linear regression, a quadratic and a cubic one. If we generate one thousand datasets, we get</p>
<p><code>SM=SIMU(1e3)<br>
S=seq(.5,5.5,by=.25)<br>
plot(S,SM[1,],type=&quot;l&quot;,ylim=c(-.1,.9),col=&quot;light blue&quot;)<br>
for(j in 2:100) lines(S,SM[j,],col=&quot;light blue&quot;)<br>
abline(h=0,lty=2)<br>
lines(S,apply(SM,2,mean),col=&quot;red&quot;,lwd=2)<br>
lines(S,apply(SM,2,function(x) quantile(x,.9)),col=&quot;red&quot;)<br>
lines(S,apply(SM,2,function(x) quantile(x,.1)),col=&quot;red&quot;)</code></p>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2016/12/rupt10.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>i.e. wherever we seek a discontinuity, we can actually find one. Even if we did generate a very smooth regression model….</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://freakonometrics.hypotheses.org/49848"> R-english – Freakonometrics</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/4qceF4tcCWo"></a></div>