Data extracted using WebPlotDigitizer.

In [94]:
mw <- read.csv("~/Downloads/nalexoneMidwest.csv", header = F)
names(mw) <- c("months", "residuals")
cutpoints <- as.data.frame(t(sapply(mw$months, function(x) x > seq(-21, 9, by = 3))))
names(cutpoints) <- (seq(-21, 9, by = 3))   
mw <- cbind.data.frame(mw, cutpoints)
In [182]:
par(mfrow=c(3,3))

## Create split plots
for (cut in seq(-18,6,by=3)) {
    plot(residuals ~ months, data = mw, col = "cadetblue3", lwd =2,
        xaxt = 'n', main = "Midwest")
    axis(1, at = seq(-24,12,by=3))
    with(mw[mw$months < cut, ], {
        pred <- predict(loess(residuals ~ months, span = 2), se=T)
        polygon(x = c(months, rev(months)), 
                y = c(pred$fit + 1.96*pred$se, 
                      rev(pred$fit - 1.96*pred$se)),
                col = adjustcolor("grey", alpha.f = .2),
                border = NA)
        lines(months, pred$fit)
    })
    
    with(mw[mw$months > cut, ], {
        pred <- predict(loess(residuals ~ months, span = 2), se=T)
        polygon(x = c(months, rev(months)), 
                y = c(pred$fit + 1.96*pred$se, 
                      rev(pred$fit - 1.96*pred$se)),
                col = adjustcolor("grey", alpha.f = .2),
                border = NA)
        lines(months, pred$fit)
    })
    
    abline(v = cut, col = "red", lty = 2, lwd = 1)
}
In [190]:
par(mfrow=c(3,4))

## Create split plots
for (cut in seq(-21,9,by=3)) {
    plot(residuals ~ months, data = mw, col = "cadetblue3", lwd =2,
        xaxt = 'n', main = "Midwest")
    axis(1, at = seq(-24,12,by=3))
    with(mw[mw$months < cut, ], {
        pred <- predict(lm(residuals ~ poly(months,1)), se=T)
        polygon(x = c(months, rev(months)), 
                y = c(pred$fit + 1.96*pred$se, 
                      rev(pred$fit - 1.96*pred$se)),
                col = adjustcolor("grey", alpha.f = .2),
                border = NA)
        lines(months, pred$fit)
    })
    
    with(mw[mw$months > cut, ], {
        pred <- predict(lm(residuals ~ poly(months,1)), se=T)
        polygon(x = c(months, rev(months)), 
                y = c(pred$fit + 1.96*pred$se, 
                      rev(pred$fit - 1.96*pred$se)),
                col = adjustcolor("grey", alpha.f = .2),
                border = NA)
        lines(months, pred$fit)
    })
    
    abline(v = cut, col = "red", lty = 2, lwd = 1)
}