如何获取ggplot中stat_ellipse椭圆最远点坐标并绘制过该点直线？

我有一段用stat_ellipse绘制椭圆的代码：

df %>% 
  select(starts_with("WindSpeed"), starts_with("Humidity")) %>%
  na.omit() %>%
  mutate(WindSpeed = (WindSpeed9am+WindSpeed3pm)/2,
         Humidity = (Humidity9am+Humidity3pm)/2) %>%
  ggplot(aes(x = WindSpeed, y = Humidity)) +
  geom_point(alpha = .4) +
  stat_ellipse(size = 2, color = "blue") + 
  geom_abline(intercept = 83, slope = -1.4, size = 2, color = "blue")

运行后得到的图像：

我想绘制一条经过该椭圆上最远两点的直线，该怎么实现？

编辑：尝试过的方法但结果不符合预期 #

我尝试用car::dataEllipse获取椭圆点，代码如下：

k <- df %>% 
  select(starts_with("WindSpeed"), starts_with("Humidity")) %>%
  na.omit() %>%
  mutate(WindSpeed = (WindSpeed9am+WindSpeed3pm)/2,
         Humidity = (Humidity9am+Humidity3pm)/2) %>%
  select(WindSpeed, Humidity) %>%
  as.matrix() %>%
  car::dataEllipse(draw = F)

k <- k$`0.5` %>% data.frame()
distances <- dist(k)
max_distance_index <- which.max(distances)
distance_matrix <- as.matrix(distances)
row_index <- row(distance_matrix)[max_distance_index]
col_index <- col(distance_matrix)[max_distance_index]
farthest_points <- k[c(row_index, col_index), ]

a <- farthest_points[1,]
b <- farthest_points[2,]
m <- (b[2]-a[2])/(b[1]-a[1])
int <- a[2]-m*a[1]

把计算出的斜率m和截距int加到ggplot后，结果不符合预期：

问题出在两个椭圆的参数不匹配：stat_ellipse默认绘制的是置信度为95%（即level=0.95）的椭圆，而你用car::dataEllipse取的是k$0.5``，也就是置信度50%的椭圆，两者完全不是同一个椭圆，自然计算出的直线不对。

另外，直接计算椭圆点对的距离效率较低，也可以通过椭圆的几何性质直接找到长轴端点——这是椭圆上距离最远的两个点。

正确步骤 #

方法1：对齐car::dataEllipse的置信度

修改代码，把k$0.5换成`k$`0.95，确保和stat_ellipse的置信度一致：

# 数据预处理和椭圆点提取
k <- df %>% 
  select(starts_with("WindSpeed"), starts_with("Humidity")) %>%
  na.omit() %>%
  mutate(WindSpeed = (WindSpeed9am+WindSpeed3pm)/2,
         Humidity = (Humidity9am+Humidity3pm)/2) %>%
  select(WindSpeed, Humidity) %>%
  as.matrix() %>%
  car::dataEllipse(draw = F, levels=0.95)  # 指定与stat_ellipse一致的置信度

k <- k$`0.95` %>% data.frame()
distances <- dist(k)
max_distance_index <- which.max(distances)
distance_matrix <- as.matrix(distances)
row_index <- row(distance_matrix)[max_distance_index]
col_index <- col(distance_matrix)[max_distance_index]
farthest_points <- k[c(row_index, col_index), ]

# 计算直线参数
a <- farthest_points[1,]
b <- farthest_points[2,]
m <- (b[2]-a[2])/(b[1]-a[1])
int <- a[2]-m*a[1]

# 绘图
df %>% 
  select(starts_with("WindSpeed"), starts_with("Humidity")) %>%
  na.omit() %>%
  mutate(WindSpeed = (WindSpeed9am+WindSpeed3pm)/2,
         Humidity = (Humidity9am+Humidity3pm)/2) %>%
  ggplot(aes(x = WindSpeed, y = Humidity)) +
  geom_point(alpha = .4) +
  stat_ellipse(size = 2, color = "blue") + 
  geom_abline(intercept = int, slope = m, size = 2, color = "red")

方法2：通过几何性质直接计算（更高效）

椭圆的长轴是数据协方差矩阵最大特征值对应的特征向量方向，长轴端点就是均值加减（最大特征值的平方根 * 对应特征向量）乘以椭圆的缩放系数（对应置信度的卡方分位数）。

library(ggplot2)
library(dplyr)

# 预处理数据
data_processed <- df %>% 
  select(starts_with("WindSpeed"), starts_with("Humidity")) %>%
  na.omit() %>%
  mutate(WindSpeed = (WindSpeed9am+WindSpeed3pm)/2,
         Humidity = (Humidity9am+Humidity3pm)/2) %>%
  select(WindSpeed, Humidity)

# 计算均值和协方差矩阵
mu <- colMeans(data_processed)
cov_mat <- cov(data_processed)

# 特征值分解，找到长轴方向
eig <- eigen(cov_mat)
max_eig_val <- max(eig$values)
max_eig_vec <- eig$vectors[, which.max(eig$values)]

# 对应95%置信度的缩放系数（自由度2的卡方分布0.95分位数的平方根）
scale_factor <- sqrt(qchisq(0.95, df=2))

# 计算长轴的两个端点（最远两点）
point1 <- mu + scale_factor * sqrt(max_eig_val) * max_eig_vec
point2 <- mu - scale_factor * sqrt(max_eig_val) * max_eig_vec

# 计算直线参数
m <- (point2[2] - point1[2])/(point2[1] - point1[1])
int <- point1[2] - m * point1[1]

# 绘图
ggplot(data_processed, aes(x = WindSpeed, y = Humidity)) +
  geom_point(alpha = .4) +
  stat_ellipse(size = 2, color = "blue") + 
  geom_abline(intercept = int, slope = m, size = 2, color = "red") +
  # 可选：标记最远两点
  geom_point(aes(x=point1[1], y=point1[2]), color="red", size=3) +
  geom_point(aes(x=point2[1], y=point2[2]), color="red", size=3)

错误原因总结 #

stat_ellipse默认level=0.95，而car::dataEllipse默认生成0.5、0.9、0.95三个置信度的椭圆，你取了0.5的那组点，相当于找的是内部小椭圆的最远点，和图中蓝色大椭圆完全无关，所以直线位置不符合预期。

内容的提问来源于stack exchange，提问作者Rootsyl