Как найти p-значение гладкой регрессии сплайна / лёсса?

10

У меня есть некоторые переменные, и мне интересно найти нелинейные отношения между ними. Поэтому я решил добавить несколько сплайнов или лессов и напечатать красивые графики (см. Код ниже). Но я также хочу иметь некоторую статистику, которая дает мне представление о том, насколько вероятно, что отношение является вопросом случайности ... т.е. мне нужно некоторое общее значение p, как, например, для линейной регрессии. Другими словами, мне нужно знать, имеет ли подобранная кривая какой-либо смысл, поскольку мой код подгонит кривую под любые данные.

x <- rnorm(1000)
y <- sin(x) + rnorm(1000, 0, 0.5)

cor.test(x,y)
plot(x, y, xlab = xlab, ylab = ylab)
spl1 <- smooth.spline(x, y, tol = 1e-6, df = 8)
lines(spl1, col = "green", lwd = 2)

spl2 <- loess(y ~ x)
x.pr <- seq(min(x), max(x), length.out = 100)
lines(x.pr, predict(spl2, x.pr), col = "blue", lwd = 2)
любознательный
источник

Ответы:

8

Библиотека шлицы имеет функции bsи nsкоторые будут создавать сплайн базиса для использования с lmфункцией, то вы можете разместить линейную модель и модель , включая шлицы и использовать anovaфункцию , чтобы сделать полный и уменьшенный тест модели , чтобы увидеть , если сплайн модель подходит значительно лучше чем линейная модель.

Вот пример кода:

x <- rnorm(1000)
y <- sin(x) + rnorm(1000, 0, 0.5)

library(splines)

fit1 <- lm(y~x)
fit0 <- lm(y~1)
fit2 <- lm(y~bs(x,5))

anova(fit1,fit2)
anova(fit0,fit2)

plot(x,y, pch='.')
abline(fit1, col='red')
xx <- seq(min(x),max(x), length.out=250)
yy <- predict(fit2, data.frame(x=xx))
lines(xx,yy, col='blue')

Вы также можете использовать polyфункцию для подбора полинома и проверки нелинейных членов как критерий кривизны.

R2

Существуют методы для вычисления и построения доверительного интервала для подбора лесса (я думаю, что в пакете ggplot2 может быть встроенный способ), вы можете построить доверительный интервал и посмотреть, подойдет ли прямая линия в пределах диапазона (это не является р-значением, но все равно дает да / нет.

Вы можете подогнать линейную модель и взять остатки и подогнать модель лесса к остаткам в качестве ответа (и представляющей интерес переменной в качестве предиктора), если истинная модель является линейной, тогда это соответствие должно быть близко к плоской линии и переупорядочивать точки по отношению к предиктору не должно иметь никакого значения. Вы можете использовать это для создания теста перестановки. Подгоните лесс, найдите предсказанное значение, самое дальнее от 0, теперь случайным образом переставьте точки и установите новый лесс и найдите самую дальнюю предсказанную точку от 0, повторите несколько раз, значение p - это доля переставленных значений, которые далее от 0, чем первоначальное значение.

Вы также можете рассмотреть перекрестную проверку как метод выбора полосы пропускания лесса. Это не дает p-значения, но бесконечная полоса пропускания соответствует идеальной линейной модели, если перекрестная проверка предполагает очень большую полосу пропускания, то это предполагает, что линейная модель может быть разумной, если более высокие полосы пропускания явно уступают некоторым из чем меньше пропускная способность, то это означает, что определенной кривизны и линейного недостаточно.

Грег Сноу
источник
Спасибо, Грег! Я думаю, что первый абзац звучит как путь, за исключением того, что я не заинтересован в сравнении с линейной моделью, просто чтобы увидеть, объясняет это сплайн или нет. Не могли бы вы предоставить код или более конкретные указатели о том, как протестировать сплайн с помощью anova? Я смотрел на функции bs и ns, но я не настолько хорош в статистике, чтобы иметь возможность сам изобрести его.
Любопытно
R2R2
anovaR2R21R2R21R2
Грег, спасибо! 1) Не могли бы вы объяснить, что lm(y~bs(x,5))делает и почему нет lm(y~I(bs(x,5)))? Я очень смущен этим вызовом, потому что результат bs (x, 5) не является переменной ... 2) Правильно ли я понимаю, что искомое значение p является результатом anova(fit0,fit2)?
Любопытно
1
xx2x3bsxlm