如何从ggplot2中的汽车包中制作椭圆？

Question

我需要构建一个图来显示普通最小二乘估计器和岭估计器之间的方差差异。为此，我想绘制方差的轮廓椭圆。

我正在使用的数据是这样的：

set.seed(10)
sigma = 10
U = rnorm(50,0,sigma) #The errors
X = scale(matrix(rnorm(50*2),ncol=2))
Y = scale(U + X%*%c(5,-2))

然后我找到岭（r）和OLS（m）及其协方差矩阵： （编辑：我已经意识到岭估计器的协方差矩阵是错误的，但这并不重要。它只会改变椭圆，而不是如何绘制。）

lbda = 100
r = solve(t(X)%*%X + lbda * diag(2))%*%t(X)%*%Y
m = solve(t(X)%*%X)%*%t(X)%*%Y
covmatr = sigma^2 * solve(t(X)%*%X + lbda * diag(2))
covmatm = sigma^2 * solve(t(X)%*%X)

然后我想构建情节。在标准 R 中，我会使用

ellipse

模块中的

car

函数来做类似的事情：

library(car)
plot(r[1],r[2],xlim=c(-5,5),ylim=c(-5,5), xlab = TeX("$\\beta_1$"), col="red")
points(m[1],m[2], col="blue")
abline(h=0)
abline(v=0)
ellipse(center=c(r),shape=covmatr,center.cex=0,radius=sqrt(qchisq(0.5,2)), col = "red")
ellipse(center=c(r),shape=covmatr,center.cex=0,radius=sqrt(qchisq(0.9,2)), col = "red")
ellipse(center=c(m),shape=covmatm,center.cex=0,radius=sqrt(qchisq(0.5,2)))
ellipse(center=c(m),shape=covmatm,center.cex=0,radius=sqrt(qchisq(0.9,2)))

我意识到我忘了重命名 y 标签

我的问题是我不知道如何在

ggplot2

中做同样的事情。我试过了：

#Building the plot with dots
rd = data.frame(r[1],r[2])
plot1 = ggplot(rd, aes(x=rd$r.1, y = rd$r.2, color = 'r')) +
    geom_point(size = 3) +
    geom_point(aes(x=m[1], y=m[2]), colour="blue", size = 3) + 
    geom_vline(xintercept = 0) +
    geom_hline(yintercept = 0) +
    xlim(-5,5) +
    ylim(-5,5) +
    ggtitle(TeX("Varians for ridge"))+
    xlab(TeX("$\\beta_1$")) + ylab("$\\beta_2$")+
    theme_bw()
plot1

#Trying to add ellipses:
ellipse(center=c(r),shape=covmatr,center.cex=0,radius=sqrt(qchisq(0.5,2)), col = "red")
ellipse(center=c(r),shape=covmatr,center.cex=0,radius=sqrt(qchisq(0.9,2)), col = "red")
ellipse(center=c(m),shape=covmatm,center.cex=0,radius=sqrt(qchisq(0.5,2)))
ellipse(center=c(m),shape=covmatm,center.cex=0,radius=sqrt(qchisq(0.9,2)))

但是这给了我以下错误

Fejl i plot.xy(xy.coords(x, y), type = type, ...) : 
  plot.new has not been called yet

我还尝试使用

car

包中的椭圆函数简单地在 ggplot 中添加省略号：

plot1 = ... +
    ellipse(center=c(m),shape=covmatm,center.cex=0,radius=sqrt(qchisq(0.9,2)))

其中“...”表示缺少行。我知道我的两种尝试都很幼稚，但我希望其中一种能够奏效。我找到了命令

stat_ellipse

但这对我来说没有用，因为据我所知它无法计算岭估计器的协方差矩阵。

Answer 1

我最近不得不做这样的事情。这有点痛苦，但是处理这个问题的方法（无需破解

stat_geom()

，或让其他人来做，从长远来看，这将是最好的解决方案）是使用

car::ellipse

生成 x 和椭圆的 y 坐标，然后用

geom_path

绘制它们，如下所示：

设置基础图

（使用OP中的数据设置，删除一些细节）

dat <- rbind(data.frame(method = "ridge", t(r)),
            data.frame(method = "OLS", t(m))) |> 
    setNames(c("method", "x", "y"))

plot1 <- ggplot(dat, aes(x, y, colour = method)) +
    geom_point(size = 3) +
    scale_colour_manual(values = c("red", "blue")) +
    xlim(-5,5) +
    ylim(-5,5)

计算椭圆

我们希望最终将所有内容都集中在一个数据框中，并通过方法和置信度进行识别

efun <- function(ctr, shp, lvl, nm) {
         e <- car::ellipse(c(ctr), shp, radius = sqrt(qchisq(lvl,2)), draw = FALSE)
         dplyr::tibble(method = nm, lvl, x = e[,"x"], y = e[,"y"])
     }
elist <- purrr::pmap(list(
        list(r, r, m, m),
        list(covmatr, covmatr, covmatm, covmatm),
        list(0.5, 0.9, 0.5, 0.9),
        list("ridge", "ridge", "OLS", "OLS")),
     efun
    )
e_df <- dplyr::bind_rows(elist)

情节

plot1 + geom_path(data = e_df, aes(lty = factor(lvl)))