如何在R中正确绘制双对数(log-log)尺度的线性拟合模型

问题 #

使用以下R代码拟合线性模型并绘图：

benchmark <- read.csv("data.csv")
model <- lm(formula = y ~ x, data = benchmark)
summary(model)
ggplot(benchmark, aes(x=x, y=y)) +
  geom_point() +
  geom_smooth(method='lm', formula = y ~ x)

绘图后发现x轴为2的幂次，数据点集中在左侧，尝试添加scale_x_log10()和scale_y_log10()后结果不符合预期。原始数据如下：

1,          3.50225e-07
2,          4.75979e-07
4,          2.65479e-07
8,          3.55327e-07
16,         2.98107e-07
32,         3.11732e-07
64,         2.40219e-07
128,        2.78199e-07
256,        3.27337e-07
512,        2.79105e-07
1024,       7.27129e-07
2048,       4.59182e-07
4096,       7.8944e-07
8192,       1.1326e-06
16384,      2.11749e-06
32768,      3.53295e-06
65536,      4.822e-06
131072,     7.3658e-06
262144,     1.2835e-05
524288,     2.34737e-05
1048576,    4.86131e-05
2097152,    9.87735e-05
4194304,    0.000215783
8388608,    0.000510538
16777216,   0.00108616
33554432,   0.00226483
67108864,   0.00480124
134217728,  0.00982463
268435456,  0.0198695
536870912,  0.0402285
1073741824, 0.0818823

寻求正确的双对数尺度线性模型绘制方法。

解决方案 #

问题根源 #

直接添加scale_x_log10()和scale_y_log10()仅会变换坐标轴显示，但geom_smooth仍基于原始数据拟合y~x模型，再将拟合线做对数变换，这并非我们需要的在对数尺度下拟合线性模型。正确的做法是拟合log(y) ~ log(x)的线性模型，再对应绘制双对数坐标轴。

具体实现 #

1. 正确读取数据

原始数据无表头，需指定列名：

benchmark <- read.csv("data.csv", header = FALSE, col.names = c("x", "y"))

2. 拟合双对数线性模型

# 拟合log(y) ~ log(x)模型
log_model <- lm(log(y) ~ log(x), data = benchmark)
summary(log_model)

3. 双对数尺度绘图（两种方式）

方式一：直接在ggplot中指定对数拟合公式

无需修改原始数据，直接在geom_smooth中指定对数变换后的公式，同时设置双对数坐标轴：

library(ggplot2)

ggplot(benchmark, aes(x = x, y = y)) +
  geom_point(size = 1.5) +
  # 指定拟合log(y) ~ log(x)的线性模型
  geom_smooth(method = "lm", formula = log(y) ~ log(x), color = "#E63946", se = TRUE) +
  # 设置x轴为2的幂次刻度，匹配数据规律
  scale_x_log10(breaks = 2^(0:30), labels = scales::comma) +
  scale_y_log10(labels = scales::scientific) +
  labs(x = "x", y = "y") +
  theme_bw()

方式二：先对数据做对数变换再绘图

先创建对数变换后的列，再拟合线性模型，绘图时可将刻度转回原始值：

# 添加对数变换列
benchmark$log_x <- log10(benchmark$x)
benchmark$log_y <- log10(benchmark$y)

# 拟合模型（与log_model等价）
log_model2 <- lm(log_y ~ log_x, data = benchmark)

ggplot(benchmark, aes(x = log_x, y = log_y)) +
  geom_point(size = 1.5) +
  geom_smooth(method = "lm", color = "#E63946", se = TRUE) +
  # 将刻度标签转回原始数值
  scale_x_continuous(labels = function(x) scales::comma(10^x)) +
  scale_y_continuous(labels = function(x) scales::scientific(10^x)) +
  labs(x = "x", y = "y") +
  theme_bw()

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