Compare la significación estadística de la diferencia entre dos regresiones polinómicas en R

Entonces, antes que nada, investigué un poco en este foro, y sé que se han hecho preguntas extremadamente similares, pero por lo general no se han respondido correctamente o, a veces, las respuestas simplemente no son lo suficientemente detalladas como para que yo las entienda. Entonces esta vez mi pregunta es: tengo dos conjuntos de datos, en cada uno, hago una regresión polinómica de la siguiente manera:

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

Las gráficas de regresiones polinómicas son:

Los coeficientes son:

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

Y ahora quiero saber, si hay una manera de usar una función R para hacer una prueba que me diga si hay una significación estadística en la diferencia entre las dos regresiones de polinomios sabiendo que el intervalo de días relevante es [ 1.100].

Por lo que entendí, no puedo aplicar directamente la prueba de anova porque los valores provienen de dos conjuntos diferentes de datos ni del AIC, que se utiliza para comparar datos modelo / verdaderos.

Traté de seguir las instrucciones dadas por @Roland en la pregunta relacionada, pero probablemente entendí mal algo al mirar mis resultados:

Aquí esta lo que hice :

Combiné mis dos conjuntos de datos en uno.

fes el factor variable del que habló @Roland. Puse 1s para el primer set y 0s para el otro.

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)) )

Mis datos se ven así ahora:

El rojo es el fit3ANOVAque aún funciona, pero tengo un problema con el azul, fit3ANOVACNel modelo tiene resultados extraños. No sé si el modelo de ajuste es correcto, no entiendo lo que @Roland quiso decir exactamente.

Considerando la solución @DeltaIV, supongo que en ese caso: los modelos son significativamente diferentes a pesar de que se superponen. ¿Tengo razón en asumir eso?

#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

Como puede ver, fit1es significativamente mejor que fit0, es decir, el efecto de la variable de agrupación es significativo. Dado que la variable de agrupación representa los conjuntos de datos respectivos, los ajustes polinómicos a los dos conjuntos de datos pueden considerarse significativamente diferentes.

La respuesta de @Ronald es la mejor y es ampliamente aplicable a muchos problemas similares (por ejemplo, ¿existe una diferencia estadísticamente significativa entre hombres y mujeres en la relación entre peso y edad?). Sin embargo, agregaré otra solución que, aunque no es tan cuantitativa (no proporciona un valor p ), ofrece una buena visualización gráfica de la diferencia.

EDITAR : según esta pregunta , parece que predict.lmla función utilizada ggplot2para calcular los intervalos de confianza no calcula las bandas de confianza simultáneas alrededor de la curva de regresión, sino solo las bandas de confianza puntuales. Estas últimas bandas no son las correctas para evaluar si dos modelos lineales ajustados son estadísticamente diferentes, o si se dice de otra manera, si podrían ser compatibles con el mismo modelo verdadero o no. Por lo tanto, no son las curvas correctas para responder a su pregunta. Como aparentemente no hay R incorporado para obtener bandas de confianza simultáneas (¡extraño!), Escribí mi propia función. Aquí está:

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)

Las bandas internas son las calculadas por defecto por geom_smooth: estas son bandas puntuales de confianza del 95% alrededor de las curvas de regresión. Las bandas externas semitransparentes (gracias por la punta de los gráficos, @Roland) son en cambio las bandas de confianza simultáneas del 95%. Como puede ver, son más grandes que las bandas puntiagudas, como se esperaba. El hecho de que las bandas de confianza simultáneas de las dos curvas no se superpongan puede tomarse como una indicación del hecho de que la diferencia entre los dos modelos es estadísticamente significativa.

Por supuesto, para una prueba de hipótesis con un valor p válido , se debe seguir el enfoque @Roland, pero este enfoque gráfico se puede ver como un análisis de datos exploratorio. Además, la trama puede darnos algunas ideas adicionales. Está claro que los modelos para los dos conjuntos de datos son estadísticamente diferentes. Pero también parece que dos modelos de grado 1 se ajustarían a los datos casi tan bien como los dos modelos cuadráticos. Podemos probar fácilmente esta hipótesis:

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

La diferencia entre el modelo de grado 1 y el modelo de grado 2 no es significativa, por lo tanto, también podemos usar dos regresiones lineales para cada conjunto de datos.