Сравните статистическую значимость разницы между двумя полиномиальными регрессиями в R

Итак, прежде всего, я провел некоторое исследование на этом форуме, и я знаю, что были заданы чрезвычайно похожие вопросы, но на них обычно не отвечали должным образом, или иногда ответ просто не был достаточно подробным, чтобы я мог понять. Итак, на этот раз мой вопрос: у меня есть два набора данных, на каждом я делаю полиномиальную регрессию следующим образом:

Ratio<-(mydata2[,c(2)])
Time_in_days<-(mydata2[,c(1)])
fit3IRC <- lm( Ratio~(poly(Time_in_days,2)) )

Графики полиномиальной регрессии:

Коэффициенты:

> as.vector(coef(fit3CN))
[1] -0.9751726 -4.0876782  0.6860041
> as.vector(coef(fit3IRC))
[1] -1.1446297 -5.4449486  0.5883757 

И теперь я хочу знать, есть ли способ использовать функцию R для выполнения теста, который бы сообщал мне, есть ли статистическая значимость в разнице между двумя регрессиями полиномов, зная, что соответствующий интервал дней равен [ 1100].

Из того, что я понял, я не могу применить тест anova напрямую, потому что значения получены из двух разных наборов данных или из AIC, который используется для сравнения данных модели / истинных данных.

Я пытался следовать инструкциям, данным @Roland в соответствующем вопросе, но, вероятно, что-то неправильно понял, глядя на мои результаты:

Вот что я сделал:

Я объединил оба моих набора данных в один.

fпеременный фактор, о котором говорил @Roland. Я поставил 1 для первого сета и 0 для другого.

y<-(mydata2[,c(2)])
x<-(mydata2[,c(1)])
f<-(mydata2[,c(3)])

plot(x,y, xlim=c(1,nrow(mydata2)),type='p')

fit3ANOVA <- lm( y~(poly(x,2)) )

fit3ANOVACN <- lm( y~f*(poly(x,2)) )

Мои данные теперь выглядят так:

Красный, fit3ANOVAкоторый все еще работает, но у меня есть проблема с синим, fit3ANOVACNу модели странные результаты. Я не знаю, верна ли подходящая модель, я не понимаю, что именно имел в виду @Roland.

Рассматривая решение @DeltaIV, я полагаю, что в этом случае: модели значительно отличаются, даже если они перекрываются. Правильно ли я так предполагаю?

#Create some example data
mydata1 <- subset(iris, Species == "setosa", select = c(Sepal.Length, Sepal.Width))
mydata2 <- subset(iris, Species == "virginica", select = c(Sepal.Length, Sepal.Width))

#add a grouping variable
mydata1$g <- "a"
mydata2$g <- "b"

#combine the datasets
mydata <- rbind(mydata1, mydata2)

#model without grouping variable
fit0 <- lm(Sepal.Width ~ poly(Sepal.Length, 2), data = mydata)

#model with grouping variable
fit1 <- lm(Sepal.Width ~ poly(Sepal.Length, 2) * g, data = mydata)

#compare models 
anova(fit0, fit1)
#Analysis of Variance Table
#
#Model 1: Sepal.Width ~ poly(Sepal.Length, 2)
#Model 2: Sepal.Width ~ poly(Sepal.Length, 2) * g
#  Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
#1     97 16.4700                                  
#2     94  7.1143  3    9.3557 41.205 < 2.2e-16 ***
#  ---
#  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Как видите, fit1значительно лучше, чем fit0, то есть влияние группирующей переменной является значительным. Так как переменная группировки представляет соответствующие наборы данных, полиномиальные соответствия для двух наборов данных могут считаться существенно различными.

Ответ @Ronald 's является лучшим и широко применим ко многим подобным проблемам (например, есть ли статистически значимая разница между мужчинами и женщинами в отношениях между весом и возрастом?). Тем не менее, я добавлю другое решение, которое, хотя и не так количественно (оно не обеспечивает p- значение), дает хорошее графическое отображение разницы.

РЕДАКТИРОВАТЬ : в соответствии с этим вопросом , похоже predict.lm, что функция, используемая ggplot2для вычисления доверительных интервалов, не вычисляет одновременные доверительные полосы вокруг кривой регрессии, а только точечные доверительные полосы. Эти последние полосы не являются правильными для оценки того, являются ли две подобранные линейные модели статистически различными или, по-другому, могут ли они быть совместимы с одной и той же истинной моделью или нет. Таким образом, они не являются правильными кривыми, чтобы ответить на ваш вопрос. Поскольку, по-видимому, нет встроенной функции R для получения одновременных доверительных интервалов (странно!), Я написал свою собственную функцию. Вот:

simultaneous_CBs <- function(linear_model, newdata, level = 0.95){
    # Working-Hotelling 1 – α confidence bands for the model linear_model
    # at points newdata with α = 1 - level

    # summary of regression model
    lm_summary <- summary(linear_model)
    # degrees of freedom 
    p <- lm_summary$df[1]
    # residual degrees of freedom
    nmp <-lm_summary$df[2]
    # F-distribution
    Fvalue <- qf(level,p,nmp)
    # multiplier
    W <- sqrt(p*Fvalue)
    # confidence intervals for the mean response at the new points
    CI <- predict(linear_model, newdata, se.fit = TRUE, interval = "confidence", 
                  level = level)
    # mean value at new points
    Y_h <- CI$fit[,1]
    # Working-Hotelling 1 – α confidence bands
    LB <- Y_h - W*CI$se.fit
    UB <- Y_h + W*CI$se.fit
    sim_CB <- data.frame(LowerBound = LB, Mean = Y_h, UpperBound = UB)
}

library(dplyr)
# sample datasets
setosa <- iris %>% filter(Species == "setosa") %>% select(Sepal.Length, Sepal.Width, Species)
virginica <- iris %>% filter(Species == "virginica") %>% select(Sepal.Length, Sepal.Width, Species)

# compute simultaneous confidence bands
# 1. compute linear models
Model <- as.formula(Sepal.Width ~ poly(Sepal.Length,2))
fit1  <- lm(Model, data = setosa)
fit2  <- lm(Model, data = virginica)
# 2. compute new prediction points
npoints <- 100
newdata1 <- with(setosa, data.frame(Sepal.Length = 
                                       seq(min(Sepal.Length), max(Sepal.Length), len = npoints )))
newdata2 <- with(virginica, data.frame(Sepal.Length = 
                                          seq(min(Sepal.Length), max(Sepal.Length), len = npoints)))
# 3. simultaneous confidence bands
mylevel = 0.95
cc1 <- simultaneous_CBs(fit1, newdata1, level = mylevel)
cc1 <- cc1 %>% mutate(Species = "setosa", Sepal.Length = newdata1$Sepal.Length)
cc2 <- simultaneous_CBs(fit2, newdata2, level = mylevel)
cc2 <- cc2 %>% mutate(Species = "virginica", Sepal.Length = newdata2$Sepal.Length)

# combine datasets
mydata <- rbind(setosa, virginica)
mycc   <- rbind(cc1, cc2)    
mycc   <- mycc %>% rename(Sepal.Width = Mean) 
# plot both simultaneous confidence bands and pointwise confidence
# bands, to show the difference
library(ggplot2)
# prepare a plot using dataframe mydata, mapping sepal Length to x,
# sepal width to y, and grouping the data by species
p <- ggplot(data = mydata, aes(x = Sepal.Length, y = Sepal.Width, color = Species)) + 
# add data points
geom_point() +
# add quadratic regression with orthogonal polynomials and 95% pointwise
# confidence intervals
geom_smooth(method ="lm", formula = y ~ poly(x,2)) +
# add 95% simultaneous confidence bands
geom_ribbon(data = mycc, aes(x = Sepal.Length, color = NULL, fill = Species, ymin = LowerBound, ymax = UpperBound),alpha = 0.5)
print(p)

Внутренние полосы - это те, которые рассчитываются по умолчанию geom_smooth: это точечные 95% доверительные полосы вокруг кривых регрессии. Внешние, полупрозрачные полосы (спасибо за графический совет, @Roland) вместо этого представляют собой одновременные 95% доверительные полосы. Как вы можете видеть, они больше, чем точечные полосы, как и ожидалось. Тот факт, что одновременные доверительные интервалы на двух кривых не перекрываются, можно считать показателем того, что различие между двумя моделями является статистически значимым.

Конечно, для проверки гипотезы с действительным p-значением следует придерживаться подхода @Roland, но этот графический подход можно рассматривать как анализ поисковых данных. Также сюжет может дать нам дополнительные идеи. Понятно, что модели для двух наборов данных статистически различны. Но похоже, что две модели степени 1 будут соответствовать данным почти так же хорошо, как две квадратичные модели. Мы можем легко проверить эту гипотезу:

fit_deg1 <- lm(data = mydata, Sepal.Width ~ Species*poly(Sepal.Length,1))
fit_deg2 <- lm(data = mydata, Sepal.Width ~ Species*poly(Sepal.Length,2))
anova(fit_deg1, fit_deg2)
# Analysis of Variance Table

# Model 1: Sepal.Width ~ Species * poly(Sepal.Length, 1)
# Model 2: Sepal.Width ~ Species * poly(Sepal.Length, 2)
#  Res.Df    RSS Df Sum of Sq      F Pr(>F)
# 1     96 7.1895                           
# 2     94 7.1143  2  0.075221 0.4969   0.61

Разница между моделью степени 1 и моделью степени 2 незначительна, поэтому мы также можем использовать две линейные регрессии для каждого набора данных.