0
$\begingroup$

I have a dataset that is split into two different behavior, a first part of the data follow a linear model and the second part of the data follows a log model. My problem is that I need to get the breakpoint where the data stop following a linear model and start following a log model. I am doing this in R. Here I give a code line to create a dataset with these characteristics:

data <- c(1:30*0.3+rnorm(30,sd=0.4), log10(0:50+10)*10+rnorm(51,sd=0.4))

This dataset has a first part that follows a linear model in a regression against its index, and the second part follows a log model in a regression against its index. My first thought was to do a linear fit until a point j and a log fit from the point j, and compared the error for the fits for different j values. But I get stacked.

$\endgroup$
  • 1
    $\begingroup$ Since these are univariate data--you just have a sequence of numbers--exactly what are you regressing and how?? $\endgroup$ – whuber May 23 '19 at 15:52
  • 1
    $\begingroup$ I am regressing the data with it index, like: lm(data~1:length(data)) $\endgroup$ – Santiago I. Hurtado May 23 '19 at 16:00
  • $\begingroup$ Thank you: I added a change-point tag to reflect that. It is likely you can find answers to your question by examining the highest-voted posts with that tag: see stats.stackexchange.com/questions/tagged/…. $\endgroup$ – whuber May 23 '19 at 16:04
2
$\begingroup$

I'm assuming that you don't know the parameters before or after the breakpoint, that you don't require continuity at the breakpoint (in your example you do have continuity at the breakpoint, but you mentioned doing two separate fits pre- and post-breakpoint, which would contradict it being a requirement), and that the standard deviation of the errors is the same before and after the breakpoint.

So I'm assuming that your observations are drawn from the following model: $$X_t \sim \begin{cases} \beta_0 + \beta_1 t + N(0, \sigma^2) & \text{if } t < \tau,\\ \beta_2 + \beta_3 \log_{10}(10 + t - \tau) + N(0, \sigma^2) & \text{if } t \geq \tau. \end{cases}$$

This is a bit different from a usual fit because you have different types of models before and after the change, but you can still calculate a best estimate for the breakpoint $\tau$. You can do this by calculating the likelihood for every value of $\tau$ and then choosing the one that maximises the likelihood.

calculateLogLikelihood = function(x, tau) {
    n = length(x)
    t = 1:n

    # define a sequence that is 1 if we are pre-breakpoint and 0 otherwise
    # (when we do the regression this will give us beta_0)
    preTauIntercept = c(rep(1, tau - 1), rep(0, n - tau + 1))
    # define a sequence that is 1, 2, 3, ... if we are pre-breakpoint and 0 otherwise
    preTauTime = c(1:(tau - 1), rep(0, n - tau + 1))
    # define a sequence that is 0 if we are pre-breakpoint and 1 otherwise
    postTauIntercept = c(rep(0, tau - 1), rep(1, n - tau + 1))
    # define a sequence that is 0 if we are pre-breakpoint and log10(10 + t - tau otherwise)
    postTauLogTime = c(rep(0, tau - 1), log10(10 + 0:(n - tau)))

    # fit a model using these 4 variables and with no default intercept
    lmTau = lm(x ~ preTauIntercept + preTauTime + postTauIntercept + postTauLogTime + 0)
    return(logLik(lmTau))
}

set.seed(1)
x = c(1:30*0.3+rnorm(30,sd=0.4), log10(0:50+10)*10+rnorm(51,sd=0.4))
# cycle through all possible values of tau, requiring at least two observations in each segment
logLikelihoods = numeric(77)
for (tau in 3:(81 - 2)) {
    logLikelihoods[tau - 2] = calculateLogLikelihood(x, tau)
}
bestTau = 2 + which.max(logLikelihoods)

# Plot the data and the best estimated breakpoint
plot(1:81, x, xlab = "t", ylab = "x[t]", pch = 16, cex = 0.5)
abline(v = bestTau)

# Plot the likelihood as a function of tau
plot(3:(81 - 2), logLikelihoods, type = "l", xlab = "Tau", ylab = "Log Likelihood")
abline(v = bestTau)

For this random seed, the estimate of $\tau$ is actually exactly correct, I get an estimate of $31$. The plots are as follows: x against time

Log likelihood against tau


Postscript Thinking about making a confidence interval for the breakpoint, I've implemented a version of the bootstraping method described in Bootstrapping confidence intervals for the change-point of time series by Marie Hušková, Claudia Kirch. The bootstrapping method is as follows:

  1. Use the data $x_1, \dots, x_n$ to calculate MLE estimates for $\tau, \beta_0, \beta_1, \beta_2, \beta_3, \sigma$.
  2. Generate $B$ data sets $x^1, x^2, \dots, x^B$ from the distribution, using MLE parameters.
  3. For each $x^b$ calculate the maximum likelihood estimate of $\tau$, $\tau^b$.
  4. Estimate the distribution and quantiles of $\tau$ from the bootstrap estimates $\tau^1, \tau^2, \dots, \tau^B$.

Looking at the output from the original fit, we have

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
preTauIntercept  0.105215   0.136669   0.770    0.444    
preTauTime       0.295340   0.007698  38.364   <2e-16 ***
postTauIntercept 0.188212   0.361496   0.521    0.604    
postTauLogTime   9.902841   0.239013  41.432   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.365 on 77 degrees of freedom

So the following code generates the bootstrap samples, which are shown in the plot below:

calculateBestTau = function(x) {
    n = length(x)
    logLikelihoods = numeric(n - 4)
    for (tau in 3:(n - 2)) {
        logLikelihoods[tau - 2] = calculateLogLikelihood(x, tau)
    }
    bestTau = 2 + which.max(logLikelihoods)
    return(bestTau)
}
# multiple sims based on MLE estimates for the parameters
B = 1000
bestTaus = numeric(B)
set.seed(1000)
for (b in 1:B) {
    xb = c(0.105215 + 1:30 * 0.295340 + rnorm(30, sd = 0.365), 
           0.188212 + log10(0:50 + 10) * 9.902841 + rnorm(51, sd = 0.365))
    bestTaus[b] = calculateBestTau(xb)
}
plot(table(bestTaus)/B, xlab = "Tau", ylab = "Probability")

rough distn of tau

And estimating a 95% confidence interval using

quantile(bestTaus, c(0.025, 0.975))

outputs $[29, 33]$.

$\endgroup$
  • $\begingroup$ +1. However, the challenge in this question is to find a way to assess the answer: that is, to be able to test hypotheses about the break point and compute confidence intervals for it. It would be nice to see that also explained. BTW, the example given by the OP is easy to fit because there is a distinct, sudden jump at time 31 of several $\sigma:$ that's why the likelihood profile drops off so quickly. $\endgroup$ – whuber May 25 '19 at 21:39
  • $\begingroup$ Following the thread to be able to test hypothesis !.Since this is time series data .....what did you do to identify and incorporate any arima structure or outliers or level shifts or local time trends in the model residuals which have been shown to impact parametric tests of significance. ? The Gaussian violations that I mentioned are not remedied by power transformations. $\endgroup$ – IrishStat May 31 '19 at 21:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.