Encontre dados e elipses de confiança (regiões?) Para uma mediana bivariada?

Estou pensando em maneiras de calcular elipses de dados e confiança em torno de uma mediana bivariada. Por exemplo, eu posso calcular facilmente uma elipse de dados ou uma elipse de confiança para a média bivariada dos seguintes dados (aqui apenas mostrando uma elipse de dados)

library("car")
set.seed(1)
df <- data.frame(x = rnorm(200, mean = 4, sd = 1.5),
                 y = rnorm(200, mean = 1.4, sd = 2.5))
plot(df)
with(df, dataEllipse(x, y, level = 0.68, add = TRUE))

Mas estou lutando com como eu faria isso por uma mediana bivariada? No caso univariado, eu poderia simplesmente inicializar a amostra novamente para gerar o intervalo necessário, mas não sei como traduzir isso no caso bivariado.

Conforme apontado por @Andy W, a mediana não está definida de forma exclusiva. Nesse caso, usamos a mediana espacial , encontrando um ponto que minimiza a norma L1 das distâncias entre as observações naquele ponto. Uma otimização foi usada para calcular a mediana espacial a partir dos pontos de dados observados.

Além disso, os pares de dados x, y no caso de uso real são dois vetores próprios de uma análise de coordenadas principais de uma matriz de dissimilaridade; portanto, xey devem ser ortogonais, se isso fornecer uma via de ataque específica.

No caso de uso real, queremos calcular a elipse de dados / confiança para grupos de pontos no espaço euclidiano. Por exemplo:

A análise é um análogo multivariado do teste de Levene de homogeneidade de variâncias entre grupos. Utilizamos medianas espaciais ou centróides de grupo padrão como a medida da tendência central multivariada e desejamos adicionar o equivalente da elipse de dados na figura acima para o caso mediano espacial.

Esta é uma boa pergunta.

Vou seguir a sugestão de @ amoeba e inicializar as medianas espaciais, usando depth::med()with method="Spatial". No entanto, há uma ligeira complicação: mednão gosta quando há pontos de dados duplicados; portanto, não podemos fazer uma inicialização direta. Em vez disso, vou desenhar uma amostra de bootstrap e, em seguida, tremer cada ponto em uma quantidade minúscula - menor que as distâncias mínimas em cada um dos$x$ e $y$ dimensões na amostra de dados original - antes de calcular a mediana espacial.

Por fim, vou calcular a menor elipse cobrindo uma proporção especificada (95%) de medianas e plotagem inicializadas .

library(depth)      # for med()
library(MASS)           # for cov.rob()
library(cluster)    # for ellipsoidhull()

# create data
set.seed(1)
df <- data.frame(x = rnorm(200, mean = 4, sd = 1.5),
                 y = rnorm(200, mean = 1.4, sd = 2.5))

# find minimum distances in each dimension for later jittering
foo <- outer(X=df$x,Y=df$x,FUN=function(xx,yy)abs(xx-yy))
delta.x <- min(foo[upper.tri(foo)])/2
foo <- outer(X=df$y,Y=df$y,FUN=function(xx,yy)abs(xx-yy))
delta.y <- min(foo[upper.tri(foo)])/2

# bootstrap spatial medians, using jittering
n.boot <- 1000
pb <- winProgressBar(max=n.boot)
boot.med <- matrix(NA,nrow=n.boot,ncol=2)
for ( ii in 1:n.boot ) {
    setWinProgressBar(pb,ii,paste(ii,"of",n.boot))
    index <- sample(1:nrow(df),nrow(df),replace=TRUE)
    bar <- df[index,] + 
      data.frame(x=runif(nrow(df),-delta.x,delta.x),
                 y=runif(nrow(df),-delta.y,delta.y))
    boot.med[ii,] <- med(bar,method="Spatial")$median
}
close(pb)

# specify confidence level
pp <- 0.95

# find smallest ellipse containing the specified proportion of bootstrapped medians
fit <- cov.rob(boot.med, quantile.used = ceiling(pp*n.boot), method = "mve")
best_ellipse <- ellipsoidhull( boot.med[fit$best,] )

plot(df)
points(boot.med,pch=19,col="grey",cex=0.5)
points(df)
lines(predict(best_ellipse), col="red")
legend("bottomright",bg="white",pch=c(21,19,NA),
    col=c("black","grey","red"),pt.bg=c("white",NA,NA),lwd=c(0,0,1),
    legend=c("Observations","Bootstrapped medians","Confidence ellipse"))

Por fim, observe que a mediana espacial bivariada é normalmente assintoticamente distribuída (Brown, 1983, JRSS, Série B ) , para que também pudéssemos dispensar o "jittered bootstrap" acima e calcular diretamente a elipse, confiando que$n=200$é "suficientemente assintótico". Posso editar este post para incluir essa elipse de confiança paramétrica, se encontrar o tempo nos próximos dias.