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
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
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
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
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] ##  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.
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.