# Estimate point in metric variable where probability of success of a conditionally Bernoulli distributed variable changes (independent observations)

There is a metric variable X and a conditionally Bernoully distributed variable Y, where the probability of success of Y changes at a threshold x of variable X. The obervations are independent. I want to estimate x together with the variance of the estimated x in order to be able to estimate a confidence interval around the estimated x. What method could be used for this setting?

A simple non-parametric technique for both testing and selecting a changepoint is to use maximally-selected test statistics. Essentially, these look for the point in x where a two-sample contrast of y is maximized. In R, this is implemented in the coin package for conditional inference in the maxstat_test() function.

Alternatively, you can also use a 0/1 dummy coding of the response and rely on structural change methods for OLS regressions. In R, the strucchange package can be leveraged for this. As this was originally developed for time series data, you should order your data by x and then regress on a constant to see whether the (segment-wise) constant mean/proportion of y changes over x.

Consider the following simple artificial example with a uniform x from [-1, 1] and a response y whose success proportion changes from 20% to 80% at x = 0.5.

set.seed(1)
x <- runif(250, -1, 1)
y <- rbinom(250, prob = ifelse(x > 0.5, 0.8, 0.2), size = 1)
plot(factor(y) ~ x)


To apply a two-sample maximally selected test statistic (essentially a type of two-sample t-test), you simply do

library("coin")
maxstat_test(factor(y) ~ x)
##  Asymptotic Generalized Maximally Selected Statistics
## data:  factor(y) by x
## maxT = 9.0955, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
##   "best" cutpoint: <= 0.4933965


For applying strucchange it is easiest to include the data in a data.frame and order it by x to resemble a time-series structure:

d <- data.frame(y = y, x = x)
d <- d[order(d$x), ]  Then you can test for a structural change using a CUSUM-type test: library("strucchange") sc <- gefp(y ~ 1, order.by = ~ x, data = d) plot(sc) sctest(sc) ## M-fluctuation test ## ## data: sc ## f(efp) = 3.8979, p-value = 1.271e-13  The corresponding plot already shows that the test is significant (because the boundary is crossed) and that the changepoint is close to 0.5 (because of the peak in the process). To formally estimate the changepoint and obtain a confidence interval you can use breakpoints(): bp <- breakpoints(y ~ 1, data = d) confint(bp) ## Confidence intervals for breakpoints ## of optimal 2-segment partition: ## ## Call: ## confint.breakpointsfull(object = bp) ## ## Breakpoints at observation number: ## 2.5 % breakpoints 97.5 % ## 1 188 191 196 ## ## Corresponding to breakdates: ## 2.5 % breakpoints 97.5 % ## 1 0.752 0.764 0.784  Here, the breakdates are simply the proportion of the data. To convert the breakpoints (= numbers of observations) to the corresponding x values one has to do this manually: d$x[confint(bp)$confint] ## [1] 0.4635850 0.4933965 0.5326213  Thus, in this simple example, the changepoint estimated by OLS regression and by the maximally-selected two-sample statistic above coincides exactly. • Thank you very much! I prefer the solution with breakpoints() from the strucchange-package in order to obtain a confidence interval (CI). I exchanged the expression confint(bp) by confint(bp, breaks=1), because breakpoints() chose breaks = 0, wherefore the CI had not been estimable by confint(bp). However, now I encountered the problem that the upper limit of the estimated CI is larger than the maximum of the observed values of x (warning message: Confidence intervals outside data time interval). How could one obtain a meaningful upper limit of the CI here? Jan 27, 2016 at 21:13 • This usually indicates that there is no strong enough evidence for a breakpoint. Is the sctest() or the maxstat_test() significant? Jan 27, 2016 at 22:59 • No, in this case they are both not significant. I did this analysis separated by gender - while the breakpoint was significant for the males (according to maxstat_test() but not according to sctest()) it was not for the females. We just want to assume that there is a breakpoint for both sexes and estimate both of them together with their confidence intervals. Interestingly doing the analysis taking males and females together gave a signicant result also for the sctest() and in the separated analysis the breakpoints were very similar. Jan 28, 2016 at 12:09 • The inference in maxstat_test is based on somewhat different principles than that of sctest, typically leading the latter to be somewhat conservative in small to moderate samples. But given your description, it seems that it is a combination of only moderately clear significance comoared with somewhat conservative inference that leads to the large confidence intervals. Jan 28, 2016 at 19:56 • Thank you again! Is there a way to obtain an upper limit of the confidence interval although it is outside the range of the observed values of x? Jan 28, 2016 at 21:04 If you have data for the$Y$, and you assume further that$X$is a beta random variable, the unconditional distribution of$Y $is a beta-binomial distribution. Also, the posterior distribution for$X $is beta as well, if you're doing inference on that. This sounds like what you want. From$p (x|y) $you can make Bayesian credible intervals. If you don't assume anything about the distribution of$X $you can still do something Bayesian. Just keep in mind that you don't think of it as a threshold random variable like$Y $, but rather as a random probability on the unit interval. Edit: ehh, I misunderstood the question. See comments below. • Thank you! Unfortunately X cannot be assumed to be beta distributed, because its values are not confined to [0,1] (X represents 'age'). May you provide a reference or key word about a corresponding Bayesian method which does not assume a specific distribution for X? Jan 25, 2016 at 14:35 • Do you have data for$X$? If you do, and you knew the threshold, you could turn$X$into a categorical variable and use logistic regression to regress$Y$on$X\$. Jan 25, 2016 at 17:24
• Yes, I have the data for X. However, unfortunately the threshold is not known. Actually, the goal is to estimate this threshold together with its standard error. For now I just optimize the common Likelihood of the threshold and the two probabilities and use bootstrap to estimate a percentile confidence-interval for the threshold. Jan 26, 2016 at 12:38