两条 geom_smooth() 行之间的区别
Difference between two geom_smooth() lines
我为我的数据绘制了一个图，现在我想得到 geom_smooth() 估计的每个 x 的 y 差值。不幸的是，有一个 similiar question 没有答案。例如，如何获得以下图（数据如下）的差异：
编辑
提了两个建议还是不知道怎么算差值
第一个建议 是从 ggplot 对象访问数据。我这样做了
pb <- ggplot_build(p)
pb[["data"]][[1]]
这种方法有点管用，但数据并未对各组使用相同的 x 值。例如，第一组的第一个x值是-3.21318853，但是第二组没有-3.21318853的x，因此，我无法计算两组之间-3.21318853的y的差异
第二个建议是看看geom_smooth()中使用了什么公式。包描述说“loess() 用于少于 1,000 个观察；否则 mgcv::gam() 与公式 = y ~ s(x, bs = "cs")”。我的N超过60000，所以默认使用gam。我不熟悉 gam； 任何人都可以提供一个简短的答案如何计算两条线之间的差异 考虑到刚才描述的事情？
R代码
library("ggplot2") # library ggplot
set.seed(1) # make example reproducible
n <- 5000 # set sample size
df <- data.frame(x= rnorm(n), g= factor(rep(c(0,1), n/2))) # generate data
df$y <- NA # include y in df
df$y[df$g== 0] <- df$x[df$g== 0]**2 + rnorm(sum(df$g== 0))*5 # y for group g= 0
df$y[df$g== 1] <-2 + df$x[df$g== 1]**2 + rnorm(sum(df$g== 1))*5 # y for g= 1 (with intercept 2)
ggplot(df, aes(x, y, col= g)) + geom_smooth() + geom_point(alpha= .1) # make a plot
您好，欢迎来到 Stack Overflow，
第一个建议很好。要使 x-sequences 匹配，您可以使用 approx 函数（在 stats 中）在两者之间插入值。
library("ggplot2") # library ggplot
set.seed(1) # make example reproducible
n <- 5000 # set sample size
df <- data.frame(x= rnorm(n), g= factor(rep(c(0,1), n/2))) # generate data
df$y <- NA # include y in df
df$y[df$g== 0] <- df$x[df$g== 0]**2 + rnorm(sum(df$g== 0))*5 # y for group g= 0
df$y[df$g== 1] <-2 + df$x[df$g== 1]**2 + rnorm(sum(df$g== 1))*5 # y for g= 1 (with intercept 2)
p <- ggplot(df, aes(x, y, col= g)) + geom_smooth() + geom_point(alpha= .1) # make a plot
pb <- ggplot_build(p) # Get computed data
data.of.g1 <- pb[['data']][[1]][pb[['data']][[1]]$group == 1, ] # Extract info for group 1
data.of.g2 <- pb[['data']][[1]][pb[['data']][[1]]$group == 2, ] # Extract info for group 2
xlimit.inf <- max(min(data.of.g1$x), min(data.of.g2$x)) # Get the minimum X the two smoothed data have in common
xlimit.sup <- min(max(data.of.g1$x), max(data.of.g2$x)) # Get the maximum X
xseq <- seq(xlimit.inf, xlimit.sup, 0.01) # Sequence of X value (you can use bigger/smaller step size)
# Based on data from group 1 and group 2, interpolates linearly for all the values in `xseq`
y.g1 <- approx(x = data.of.g1$x, y = data.of.g1$y, xout = xseq)
y.g2 <- approx(x = data.of.g2$x, y = data.of.g2$y, xout = xseq)
difference <- data.frame(x = xseq, dy = abs(y.g1$y - y.g2$y)) # Compute the difference
ggplot(difference, aes(x = x, y = dy)) + geom_line() # Make the plot
输出：
正如我在上面的评论中提到的，你真的最好在 ggplot 之外做这件事，而是用你可以计算的两个平滑的完整模型来做差异等方面的不确定性
这基本上是我大约一年前写的 blog post 的简短版本。
OP 的示例数据
set.seed(1) # make example reproducible
n <- 5000 # set sample size
df <- data.frame(x= rnorm(n), g= factor(rep(c(0,1), n/2))) # generate data
df$y <- NA # include y in df
df$y[df$g== 0] <- df$x[df$g== 0]**2 + rnorm(sum(df$g== 0))*5 # y for group g= 0
df$y[df$g== 1] <-2 + df$x[df$g== 1]**2 + rnorm(sum(df$g== 1))*5 # y for g= 1 (with intercept 2)
首先为示例数据拟合模型：
library("mgcv")
m <- gam(y ~ g + s(x, by = g), data = df, method = "REML")
这里我用 factor-smooth 交互（by 位）拟合 GAM，对于这个模型，我们还需要将 g 作为参数效果包含在 group-specific 平滑都以 0 为中心，因此我们需要在模型的参数部分包含组均值。
接下来我们需要沿着 x 变量的数据网格，我们将在其中估计两个估计平滑度之间的差异：
pdat <- with(df, expand.grid(x = seq(min(x), max(x), length = 200),
g = c(0,1)))
pdat <- transform(pdat, g = factor(g))
然后我们使用此预测数据生成 Xp 矩阵，该矩阵将协变量的值映射到平滑的基础扩展的值；我们可以操纵这个矩阵来平滑我们想要的差异：
xp <- predict(m, newdata = pdat, type = "lpmatrix")
接下来是一些代码，用于识别 xp 中的哪些行和列属于 g 各个级别的平滑；由于模型中只有两个级别并且只有一个平滑项，这完全是微不足道的，但对于更复杂的模型来说，这是必需的，并且重要的是要使 grep() 位工作的平滑组件名称正确。
## which cols of xp relate to splines of interest?
c1 <- grepl('g0', colnames(xp))
c2 <- grepl('g1', colnames(xp))
## which rows of xp relate to sites of interest?
r1 <- with(pdat, g == 0)
r2 <- with(pdat, g == 1)
现在我们可以区分 xp 的行来比较我们正在比较的水平对
## difference rows of xp for data from comparison
X <- xp[r1, ] - xp[r2, ]
当我们关注差异时，我们需要将与所选平滑对无关的所有列清零，其中包括任何参数项。
## zero out cols of X related to splines for other lochs
X[, ! (c1 | c2)] <- 0
## zero out the parametric cols
X[, !grepl('^s\(', colnames(xp))] <- 0
（在这个例子中，这两行做完全一样的事情，但在更复杂的例子中，这两行都需要。）
现在我们有一个矩阵 X，它包含我们感兴趣的一对平滑的两个基础扩展之间的差异，但是要根据响应的拟合值得到它 y 我们需要将该矩阵乘以系数向量：
## difference between smooths
dif <- X %*% coef(m)
现在 dif 包含两个平滑之间的差异。
我们可以再次使用 X 和模型系数的协方差矩阵来计算此差异的标准误差，从而计算出估计差异的 95%（在本例中）置信区间。
## se of difference
se <- sqrt(rowSums((X %*% vcov(m)) * X))
## confidence interval on difference
crit <- qt(.975, df.residual(m))
upr <- dif + (crit * se)
lwr <- dif - (crit * se)
请注意，在 vcov() 调用中，我们使用的是经验贝叶斯协方差矩阵，而不是因选择平滑度参数而校正的协方差矩阵。我很快展示的函数允许您通过参数 unconditional = TRUE.
来解释这种额外的不确定性
最后我们收集结果并绘制：
res <- data.frame(x = with(df, seq(min(x), max(x), length = 200)),
dif = dif, upr = upr, lwr = lwr)
ggplot(res, aes(x = x, y = dif)) +
geom_ribbon(aes(ymin = lwr, ymax = upr, x = x), alpha = 0.2) +
geom_line()
这会产生
这与评估结果一致，该评估显示具有 group-level 平滑度的模型并未提供比具有不同组均值但仅在 x 中仅具有单一公共平滑度的模型更好的拟合：
r$> m0 <- gam(y ~ g + s(x), data = df, method = "REML")
r$> AIC(m0, m)
df AIC
m0 9.68355 30277.93
m 14.70675 30285.02
r$> anova(m0, m, test = 'F')
Analysis of Deviance Table
Model 1: y ~ g + s(x)
Model 2: y ~ g + s(x, by = g)
Resid. Df Resid. Dev Df Deviance F Pr(>F)
1 4990.1 124372
2 4983.9 124298 6.1762 73.591 0.4781 0.8301
总结
我提到的博客 post 有一个函数将上面的步骤包装成一个简单的函数，smooth_diff():
smooth_diff <- function(model, newdata, f1, f2, var, alpha = 0.05,
unconditional = FALSE) {
xp <- predict(model, newdata = newdata, type = 'lpmatrix')
c1 <- grepl(f1, colnames(xp))
c2 <- grepl(f2, colnames(xp))
r1 <- newdata[[var]] == f1
r2 <- newdata[[var]] == f2
## difference rows of xp for data from comparison
X <- xp[r1, ] - xp[r2, ]
## zero out cols of X related to splines for other lochs
X[, ! (c1 | c2)] <- 0
## zero out the parametric cols
X[, !grepl('^s\(', colnames(xp))] <- 0
dif <- X %*% coef(model)
se <- sqrt(rowSums((X %*% vcov(model, unconditional = unconditional)) * X))
crit <- qt(alpha/2, df.residual(model), lower.tail = FALSE)
upr <- dif + (crit * se)
lwr <- dif - (crit * se)
data.frame(pair = paste(f1, f2, sep = '-'),
diff = dif,
se = se,
upper = upr,
lower = lwr)
}
使用这个函数我们可以重复整个分析并绘制差异：
out <- smooth_diff(m, pdat, '0', '1', 'g')
out <- cbind(x = with(df, seq(min(x), max(x), length = 200)),
out)
ggplot(out, aes(x = x, y = diff)) +
geom_ribbon(aes(ymin = lower, ymax = upper, x = x), alpha = 0.2) +
geom_line()
我不会在这里显示绘图，因为它与上面显示的相同，除了轴标签。
我为我的数据绘制了一个图，现在我想得到 geom_smooth() 估计的每个 x 的 y 差值。不幸的是，有一个 similiar question 没有答案。例如，如何获得以下图（数据如下）的差异：
编辑
提了两个建议还是不知道怎么算差值
第一个建议 是从 ggplot 对象访问数据。我这样做了
pb <- ggplot_build(p)
pb[["data"]][[1]]
这种方法有点管用，但数据并未对各组使用相同的 x 值。例如，第一组的第一个x值是-3.21318853，但是第二组没有-3.21318853的x，因此，我无法计算两组之间-3.21318853的y的差异
第二个建议是看看geom_smooth()中使用了什么公式。包描述说“loess() 用于少于 1,000 个观察；否则 mgcv::gam() 与公式 = y ~ s(x, bs = "cs")”。我的N超过60000，所以默认使用gam。我不熟悉 gam； 任何人都可以提供一个简短的答案如何计算两条线之间的差异 考虑到刚才描述的事情？
R代码
library("ggplot2") # library ggplot
set.seed(1) # make example reproducible
n <- 5000 # set sample size
df <- data.frame(x= rnorm(n), g= factor(rep(c(0,1), n/2))) # generate data
df$y <- NA # include y in df
df$y[df$g== 0] <- df$x[df$g== 0]**2 + rnorm(sum(df$g== 0))*5 # y for group g= 0
df$y[df$g== 1] <-2 + df$x[df$g== 1]**2 + rnorm(sum(df$g== 1))*5 # y for g= 1 (with intercept 2)
ggplot(df, aes(x, y, col= g)) + geom_smooth() + geom_point(alpha= .1) # make a plot
您好，欢迎来到 Stack Overflow，
第一个建议很好。要使 x-sequences 匹配，您可以使用 approx 函数（在 stats 中）在两者之间插入值。
library("ggplot2") # library ggplot
set.seed(1) # make example reproducible
n <- 5000 # set sample size
df <- data.frame(x= rnorm(n), g= factor(rep(c(0,1), n/2))) # generate data
df$y <- NA # include y in df
df$y[df$g== 0] <- df$x[df$g== 0]**2 + rnorm(sum(df$g== 0))*5 # y for group g= 0
df$y[df$g== 1] <-2 + df$x[df$g== 1]**2 + rnorm(sum(df$g== 1))*5 # y for g= 1 (with intercept 2)
p <- ggplot(df, aes(x, y, col= g)) + geom_smooth() + geom_point(alpha= .1) # make a plot
pb <- ggplot_build(p) # Get computed data
data.of.g1 <- pb[['data']][[1]][pb[['data']][[1]]$group == 1, ] # Extract info for group 1
data.of.g2 <- pb[['data']][[1]][pb[['data']][[1]]$group == 2, ] # Extract info for group 2
xlimit.inf <- max(min(data.of.g1$x), min(data.of.g2$x)) # Get the minimum X the two smoothed data have in common
xlimit.sup <- min(max(data.of.g1$x), max(data.of.g2$x)) # Get the maximum X
xseq <- seq(xlimit.inf, xlimit.sup, 0.01) # Sequence of X value (you can use bigger/smaller step size)
# Based on data from group 1 and group 2, interpolates linearly for all the values in `xseq`
y.g1 <- approx(x = data.of.g1$x, y = data.of.g1$y, xout = xseq)
y.g2 <- approx(x = data.of.g2$x, y = data.of.g2$y, xout = xseq)
difference <- data.frame(x = xseq, dy = abs(y.g1$y - y.g2$y)) # Compute the difference
ggplot(difference, aes(x = x, y = dy)) + geom_line() # Make the plot
输出：
正如我在上面的评论中提到的，你真的最好在 ggplot 之外做这件事，而是用你可以计算的两个平滑的完整模型来做差异等方面的不确定性
这基本上是我大约一年前写的 blog post 的简短版本。
OP 的示例数据
set.seed(1) # make example reproducible
n <- 5000 # set sample size
df <- data.frame(x= rnorm(n), g= factor(rep(c(0,1), n/2))) # generate data
df$y <- NA # include y in df
df$y[df$g== 0] <- df$x[df$g== 0]**2 + rnorm(sum(df$g== 0))*5 # y for group g= 0
df$y[df$g== 1] <-2 + df$x[df$g== 1]**2 + rnorm(sum(df$g== 1))*5 # y for g= 1 (with intercept 2)
首先为示例数据拟合模型：
library("mgcv")
m <- gam(y ~ g + s(x, by = g), data = df, method = "REML")
这里我用 factor-smooth 交互（by 位）拟合 GAM，对于这个模型，我们还需要将 g 作为参数效果包含在 group-specific 平滑都以 0 为中心，因此我们需要在模型的参数部分包含组均值。
接下来我们需要沿着 x 变量的数据网格，我们将在其中估计两个估计平滑度之间的差异：
pdat <- with(df, expand.grid(x = seq(min(x), max(x), length = 200),
g = c(0,1)))
pdat <- transform(pdat, g = factor(g))
然后我们使用此预测数据生成 Xp 矩阵，该矩阵将协变量的值映射到平滑的基础扩展的值；我们可以操纵这个矩阵来平滑我们想要的差异：
xp <- predict(m, newdata = pdat, type = "lpmatrix")
接下来是一些代码，用于识别 xp 中的哪些行和列属于 g 各个级别的平滑；由于模型中只有两个级别并且只有一个平滑项，这完全是微不足道的，但对于更复杂的模型来说，这是必需的，并且重要的是要使 grep() 位工作的平滑组件名称正确。
## which cols of xp relate to splines of interest?
c1 <- grepl('g0', colnames(xp))
c2 <- grepl('g1', colnames(xp))
## which rows of xp relate to sites of interest?
r1 <- with(pdat, g == 0)
r2 <- with(pdat, g == 1)
现在我们可以区分 xp 的行来比较我们正在比较的水平对
## difference rows of xp for data from comparison
X <- xp[r1, ] - xp[r2, ]
当我们关注差异时，我们需要将与所选平滑对无关的所有列清零，其中包括任何参数项。
## zero out cols of X related to splines for other lochs
X[, ! (c1 | c2)] <- 0
## zero out the parametric cols
X[, !grepl('^s\(', colnames(xp))] <- 0
（在这个例子中，这两行做完全一样的事情，但在更复杂的例子中，这两行都需要。）
现在我们有一个矩阵 X，它包含我们感兴趣的一对平滑的两个基础扩展之间的差异，但是要根据响应的拟合值得到它 y 我们需要将该矩阵乘以系数向量：
## difference between smooths
dif <- X %*% coef(m)
现在 dif 包含两个平滑之间的差异。
我们可以再次使用 X 和模型系数的协方差矩阵来计算此差异的标准误差，从而计算出估计差异的 95%（在本例中）置信区间。
## se of difference
se <- sqrt(rowSums((X %*% vcov(m)) * X))
## confidence interval on difference
crit <- qt(.975, df.residual(m))
upr <- dif + (crit * se)
lwr <- dif - (crit * se)
请注意，在 vcov() 调用中，我们使用的是经验贝叶斯协方差矩阵，而不是因选择平滑度参数而校正的协方差矩阵。我很快展示的函数允许您通过参数 unconditional = TRUE.
最后我们收集结果并绘制：
res <- data.frame(x = with(df, seq(min(x), max(x), length = 200)),
dif = dif, upr = upr, lwr = lwr)
ggplot(res, aes(x = x, y = dif)) +
geom_ribbon(aes(ymin = lwr, ymax = upr, x = x), alpha = 0.2) +
geom_line()
这会产生
这与评估结果一致，该评估显示具有 group-level 平滑度的模型并未提供比具有不同组均值但仅在 x 中仅具有单一公共平滑度的模型更好的拟合：
r$> m0 <- gam(y ~ g + s(x), data = df, method = "REML")
r$> AIC(m0, m)
df AIC
m0 9.68355 30277.93
m 14.70675 30285.02
r$> anova(m0, m, test = 'F')
Analysis of Deviance Table
Model 1: y ~ g + s(x)
Model 2: y ~ g + s(x, by = g)
Resid. Df Resid. Dev Df Deviance F Pr(>F)
1 4990.1 124372
2 4983.9 124298 6.1762 73.591 0.4781 0.8301
总结
我提到的博客 post 有一个函数将上面的步骤包装成一个简单的函数，smooth_diff():
smooth_diff <- function(model, newdata, f1, f2, var, alpha = 0.05,
unconditional = FALSE) {
xp <- predict(model, newdata = newdata, type = 'lpmatrix')
c1 <- grepl(f1, colnames(xp))
c2 <- grepl(f2, colnames(xp))
r1 <- newdata[[var]] == f1
r2 <- newdata[[var]] == f2
## difference rows of xp for data from comparison
X <- xp[r1, ] - xp[r2, ]
## zero out cols of X related to splines for other lochs
X[, ! (c1 | c2)] <- 0
## zero out the parametric cols
X[, !grepl('^s\(', colnames(xp))] <- 0
dif <- X %*% coef(model)
se <- sqrt(rowSums((X %*% vcov(model, unconditional = unconditional)) * X))
crit <- qt(alpha/2, df.residual(model), lower.tail = FALSE)
upr <- dif + (crit * se)
lwr <- dif - (crit * se)
data.frame(pair = paste(f1, f2, sep = '-'),
diff = dif,
se = se,
upper = upr,
lower = lwr)
}
使用这个函数我们可以重复整个分析并绘制差异：
out <- smooth_diff(m, pdat, '0', '1', 'g')
out <- cbind(x = with(df, seq(min(x), max(x), length = 200)),
out)
ggplot(out, aes(x = x, y = diff)) +
geom_ribbon(aes(ymin = lower, ymax = upper, x = x), alpha = 0.2) +
geom_line()
我不会在这里显示绘图，因为它与上面显示的相同，除了轴标签。