Different ways to show confidence interval as an ellipse or a minimum convex polygon, with gradient of colors

Let calculate of sample of 100000 points defining credibility interval from MCMC with Metropolis-Hasting algorithm based on a logistic model:

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, 11, 12, 13, 5, 3, 1, 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"])

Let draw the results using different methods.

First draw simply all the points:

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

Now add a 95% 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))

Now a 99% ellipse with gradient color:

plot(x=xMCMC, 
     y=yMCMC, 
     pch=".", col=rgb(red = 0, green = 0, blue=0, alpha=0.01), 
     bty="n", las=1, xlab="P", ylab="S", 
     xlim=c(25, 33), ylim=c(-2, 0), type="p")


n.levels <- 128
col.levels <- grey.colors(n.levels)
ellipse99gradient <- dataEllipse(x=xMCMC, 
                                 y=yMCMC, 
                                       levels=seq(from=0.99, to=0.01, length=n.levels), 
                                       draw=FALSE)

for (ellip in 1:n.levels) {
  polygon(ellipse99gradient[[ellip]][, 1], ellipse99gradient[[ellip]][, 2], col=rev(col.levels)[ellip], border = NA)
}

points(x=xMCMC, 
        y=yMCMC, 
        pch=".", col=rgb(red = 1, green = 0, blue=0, alpha=0.02))

Now do a Minimum Convex Polygon that encompasses all the points:

plot(x=xMCMC, 
     y=yMCMC, 
     pch=".", col=rgb(red = 0, green = 0, blue=0, alpha=0.01), 
     bty="n", las=1, xlab="P", ylab="S", 
     xlim=c(25, 33), ylim=c(-2, 0), type="p")


library("spatstat")
mcp_bypnts <- convexhull.xy(x=xMCMC, 
                            y=yMCMC)
polygon(x=mcp_bypnts$bdry[[1]]$x  , y=mcp_bypnts$bdry[[1]]$y, border=NA, col=rgb(0.5, 0, 0, 0.2))

Let the Minimum Convex Polygon encompasses only 95% of the points:

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)

dst <- sqrt((sx-mx)^2+(sy-my)^2)

ldst <- quantile(dst, 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.5, 0, 0, 0.2))

Note that the shape converge to a circle (in standardized scale, then an oval in non standardized scale).

Let converge to a square (in standardized scale, then an parallelpiped in non standardized scale).

dst <- abs(sx-mx)+abs(sy-my)

ldst <- quantile(dst, 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, 0, 0.5, 0.2))

Now do a gradient color for the Minimum Convex Polygon:

plot(x=xMCMC, 
     y=yMCMC, 
     pch=".", col=rgb(red = 0, green = 0, blue=0, alpha=0.01), 
     bty="n", las=1, xlab="P", ylab="S", 
     xlim=c(25, 33), ylim=c(-2, 0), type="p")


n.levels <- 128
col.levels <- rev(grey.colors(n.levels))

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)

# This formula implies that it is a circle because it is applied on
# standardized values
dst <- sqrt((sx-mx)^2+(sy-my)^2)
probs <- seq(from=1, to=0.01, length=n.levels)

for (k in 1:n.levels) {
  # print(k)
  ellip <- probs[k]
  
  ldst <- quantile(dst, ellip)
  
  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=col.levels[k])
}