Credibility interval at x% when the points do not fit an ellipse

Imagine a logistic model with these data:

t <- c(25, 25, 26, 28, 30, 32, 33, 33, 34)
n <- c(10, 12, 13, 15, 10, 16, 10, 12, 15)
x <- c(10, 12, 13, 15, 0, 0, 0, 0, 0)

S <- -0.5
P <- 29

Lbinom <- function(data, x) -sum(dbinom(x=data$x, size=data$size, prob=1/(1+exp((1/x["S"])*(x["P"]-data$t))), log=TRUE))

Lbinom(data=list(x=x, size=n, t=t), x=c(P=P, S=S))

library("HelpersMG")

parameters_mcmc <- data.frame(Density=c('dunif', 'dunif'), 
            Prior1=c(25, -2), Prior2=c(32, 2), SDProp=c(0.35, 0.2), 
            Min=c(25, -2), Max=c(32, 2), Init=c(29, -0.5), stringsAsFactors = FALSE, 
            row.names=c('P', 'S'))

mcmc_run <- MHalgoGen(n.iter=50000, parameters=parameters_mcmc, data=list(x=x, size=n, t=t), 
                      likelihood=Lbinom, n.chains=1, n.adapt=100, thin=1, trace=1, adaptive = TRUE)

plot(mcmc_run, parameters = 1, xlim=c(20, 35))
plot(mcmc_run, parameters = 2, xlim=c(-2.5, 2.5))

raftery.diag(as.mcmc(mcmc_run))
heidel.diag(as.mcmc(mcmc_run))
1-rejectionRate(as.mcmc(mcmc_run))

xMCMC <- as.numeric(mcmc_run$resultMCMC[[1]][, "P"])
yMCMC <- as.numeric(mcmc_run$resultMCMC[[1]][, "S"])

plot(x=xMCMC, 
     y=yMCMC, 
     pch=".", col="grey", bty="n", las=1, xlab="P", ylab="S", 
     xlim=c(25, 33), ylim=c(-1, 0), cex=0.1)

Clearly the distribution of values for the parameters do not fit with an ellipse.
library("car")

ellipse95 <- dataEllipse(x=xMCMC, 
                                       y=yMCMC, 
                                       levels=c(0.95), 
                                       draw=FALSE)
polygon(ellipse95[, 1], ellipse95[, 2], border = NA, col=rgb(0, 0.5, 0, 0.5))

Let do using Maximun Convex Polygon but with a special strategy that will work even when data have special shape.

The classical strategy in such a case is to remove the points based on their distance to the gravity center. However, automatically the shape will be an oval.

The original strategy I propose here is to remove the points which are the closest to the edges of the minimum convex polygon. Then the shape of the remaining points is not constrained to be a simple shape.

plot(x=xMCMC, 

     y=yMCMC, 

     pch=".", col="grey", bty="n", las=1, xlab="P", ylab="S", 

     xlim=c(25, 33), ylim=c(-1, 0.1), cex=0.1)

sx <- xMCMC

sy <- yMCMC

mx0 <- mean(xMCMC)

my0 <- mean(yMCMC)

minx <- min(xMCMC); maxx <- max(xMCMC)

sx <- (xMCMC-mx0)/(maxx-minx)

miny <- min(yMCMC); maxy <- max(yMCMC)

sy <- (yMCMC-my0)/(maxy-miny)

mx <- mean(sx)

my <- mean(sy)

# je calcule les coordonnées du MCP qui inclut tous les points

mcp_bypnts0 <- convexhull.xy(x=sx, y=sy)

bx <- mcp_bypnts0$bdry[[1]]$x

by <- mcp_bypnts0$bdry[[1]]$y

polygon(x=bx*(maxx-minx)+mx0, y=by*(maxy-miny)+my0, border=NA, col=rgb(0, 0.8, 0, 0.2))

bx_c <- NULL

by_c <- NULL

by_x <- 0.0001

for (i in 1:length(bx)) {

  bxec <- ifelse(i==1, tail(bx, n=1), bx[i-1])

  byec <- ifelse(i==1, tail(by, n=1), by[i-1])

  m <- lm(y ~ x, data=data.frame(x=c(bxec, bx[i]), y=c(byec, by[i])))

  bx_c <- c(bx_c, seq(from=bxec, to=bx[i], by=sign(bx[i]-bxec)*by_x))

  by_c <- c(by_c, predict(m, newdata = data.frame(x=seq(from=bxec, to=bx[i], by=sign(bx[i]-bxec)*by_x))))

}

points(x=bx_c*(maxx-minx)+mx0, y=by_c*(maxy-miny)+my0, pch=19, cex=0.01, col="red")

points(x=bx*(maxx-minx)+mx0, y=by*(maxy-miny)+my0, pch=19, cex=0.5)

Now make the polygon being smaller to retain only a fraction x% of the original points.

# For each point, I estimate the closest distance to the border

dx <- sapply(sx, FUN=function(x) min((x-bx_c)^2))

dy <- sapply(sy, FUN=function(y) min((y-by_c)^2))

dst <- sqrt(dx+dy)

# I delete the 5% the closest of the border. So it will be the 95% MCP

ldst <- quantile(dst, 1-0.95)

mcp_bypnts <- convexhull.xy(x=xMCMC[dst>=ldst], 

                            y=yMCMC[dst>=ldst])

polygon(x=mcp_bypnts$bdry[[1]]$x  ,y=mcp_bypnts$bdry[[1]]$y, border=NA, col=rgb(0.2, 0, 0, 0.2))

# points(x=sx0[dst>=ldst], 

#        y=sy0[dst>=ldst], pch=".", col="blue")

# I delete the 15% the closest of the border. So it will be the 85% MCP

ldst <- quantile(dst, 1-0.85)

mcp_bypnts <- convexhull.xy(x=xMCMC[dst>=ldst], 

                            y=yMCMC[dst>=ldst])

# points(x=sx0[dst>=ldst], 

#        y=sy0[dst>=ldst], pch=".", col="blue")

polygon(x=mcp_bypnts$bdry[[1]]$x, y=mcp_bypnts$bdry[[1]]$y, border=NA, col=rgb(0.2, 0, 0, 0.2))

ldst <- quantile(dst, 1-0.5)

mcp_bypnts <- convexhull.xy(x=xMCMC[dst>=ldst], 

                            y=yMCMC[dst>=ldst])

polygon(x=mcp_bypnts$bdry[[1]]$x, y=mcp_bypnts$bdry[[1]]$y, border=NA, col=rgb(0.2, 0, 0, 0.2))

ldst <- quantile(dst, 1-0.1)

mcp_bypnts <- convexhull.xy(x=xMCMC[dst>=ldst], 

                            y=yMCMC[dst>=ldst])

# points(x=sx0[dst>=ldst], 

#        y=sy0[dst>=ldst], pch=".", col="green")

polygon(x=mcp_bypnts$bdry[[1]]$x, y=mcp_bypnts$bdry[[1]]$y, border=NA, col=rgb(0.2, 0., 0, 0.2))