Как я могу получить значение x через заданное значение y при использовании языка R для линейной регрессии?
У меня есть следующая таблица данных. Я визуализировал это с помощью ggplot2. Мне интересно, когда y = 1, каковы соответствующие значения оси X, соответствующие нижней границе стандартной ошибки (se) и кривой регрессии?
library(tidyverse)
set.seed(114)
tableInput <- tibble(
x = sample(60:270,100,replace = T) %>% sort(decreasing = F),
y = rnorm(100,1,1) %>% abs(),
)
ggplot(tableInput, aes(x,y)) +
geom_smooth(method = lm) +
geom_hline(yintercept = 1) +
scale_x_continuous(n.breaks = 10) +
coord_cartesian(ylim = c(.5,2),xlim = c(50,280)) +
theme_classic()
см. этот пост
Ответы 1
Решение в базе R состоит в том, чтобы использовать optimize, чтобы найти значение x, которое дает желаемое значение y.
Начнем с написания функции, которая будет минимальной, когда линия регрессии находится в точке y = 1.
f1 <- function(val) {
mod <- lm(y ~ x, tableInput)
pred <- predict(mod, newdata = data.frame(x = val), se = TRUE)
(1 - pred$fit)^2
}
Мы можем сделать то же самое для нижней границы ленты SE следующим образом:
f2 <- function(val) {
mod <- lm(y ~ x, tableInput)
pred <- predict(mod, newdata = data.frame(x = val), se = TRUE)
(1 - (pred$fit + qnorm(0.025) * pred$se.fit))^2
}
Теперь мы можем использовать optimize, чтобы автоматически найти правильные значения x:
x_regress <- optimize(f1, c(50, 275))$minimum
x_lower <- optimize(f2, c(50, 275))$minimum
x_regress
#> [1] 208.0677
x_lower
#> [1] 157.8854
И мы можем подтвердить их правильность, добавив их в качестве аннотаций к сюжету:
ggplot(tableInput, aes(x,y)) +
geom_smooth(method = lm) +
geom_hline(yintercept = 1) +
scale_x_continuous(n.breaks = 10) +
annotate('point', x = c(x_regress, x_lower), y = 1, col = 'red') +
coord_cartesian(ylim = c(.5,2),xlim = c(50,280)) +
theme_classic()
Инвертирование линейной функции не требует числовой оптимизации. Вам нужно только применить алгебру на уровне средней школы: (y - coef(mod)[[1]]) / coef(mod)[[2]]. В качестве доверительного интервала я бы использовал qt(0.025, df = mod$df.residual) на случай, если он будет использоваться для модели, подходящей для небольшого количества наблюдений.
@Роланд, спасибо - в данном случае это, конечно, правда. Преимущество использования числовой оптимизации заключается в том, что оно применимо к более широкому набору моделей.
Спасибо за ваш ответ. Я хочу знать, в чем заключается эффект qnorm(0.025) * pred$se.fit(1 - (pred$fit + qnorm(0.025) * pred$se.fit))^2 . Почему pred$se.fit нужно умножать на qnorm(0.025) ?
@tenjyosyu на ленте показан доверительный интервал 95% для линии регрессии. 95%-ный доверительный интервал находится между квантилями 2,5% и 97,5% нормального распределения, среднее значение которого соответствует подгонке модели, а стандартное отклонение представляет собой подобранную стандартную ошибку. qnorm(0.025) равен примерно -1,96, другими словами, это количество стандартных ошибок ниже подходящей модели, где мы найдем нижнюю границу ленты.