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

Question

Así que en primer lugar hice algunas investigaciones en este foro, y sé que extremadamente similar Se han formulado preguntas, pero normalmente no se han respondido adecuadamente o, a veces, la respuesta simplemente no es lo suficientemente detallada como para que yo la entienda. Así que esta vez mi pregunta es : Tengo dos conjuntos de datos, en cada uno, hago una regresión polinómica así :

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

Los gráficos de las 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 alguna forma de usar una función de R para hacer un test que me diga si hay o no significación estadística en la diferencia entre los dos polinomios de regresión sabiendo que el intervalo relevante de días es [1,100].

Por lo que he entendido no puedo aplicar directamente el test anova porque los valores proceden de dos conjuntos de datos diferentes ni el AIC, que se utiliza para comparar modelo/datos verdaderos.

Intenté seguir las instrucciones dadas por @Roland en la pregunta relacionada pero probablemente entendí algo mal al ver mis resultados :

Esto es lo que hice:

He combinado mis dos conjuntos de datos en uno solo.

f es el factor variable del que hablaba @Roland. Puse 1s para el primer conjunto 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 tienen ahora este aspecto :

El rojo es fit3ANOVA que sigue funcionando pero tengo un problema con el azul fit3ANOVACN el modelo tiene resultados extraños. No sé si el modelo de ajuste es correcto, no entiendo lo que @Roland quería decir exactamente.

Teniendo en cuenta la solución de @DeltaIV supongo que en ese caso : Los modelos son significativamente diferentes aunque se solapen. ¿Estoy en lo cierto?

Answer 1

#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 ves, fit1 es 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.

Answer 2

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

EDITAR según esta pregunta Parece que predict.lm la función utilizada por ggplot2 para calcular los intervalos de confianza, no calcula bandas de confianza simultáneas alrededor de la curva de regresión, sino sólo bandas de confianza puntuales. Estas últimas bandas no son las adecuadas para evaluar si dos modelos lineales ajustados son estadísticamente diferentes, o dicho de otro modo, si podrían ser compatibles con el mismo modelo verdadero o no. Por lo tanto, no son las curvas adecuadas para responder a su pregunta. Como aparentemente no existe una función de R para obtener bandas de confianza simultáneas (¡extraño!), escribí mi propia función. Aquí la tienes:

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 interiores son las calculadas por defecto por geom_smooth : estos son puntualmente Bandas de confianza del 95% alrededor de las curvas de regresión. Las bandas exteriores semitransparentes (gracias por el consejo sobre los gráficos, @Roland ) son en cambio el simultánea Bandas de confianza del 95%. Como puede ver, son mayores que las bandas puntuales, como era de esperar. El hecho de que las bandas de confianza simultáneas de las dos curvas no se solapen puede tomarse como una indicación de que la diferencia entre los dos modelos es estadísticamente significativa.

Por supuesto, para una prueba de hipótesis con un p -valor, debe seguirse el enfoque @Roland, pero este enfoque gráfico puede considerarse un análisis exploratorio de datos. Además, el gráfico puede darnos algunas ideas adicionales. Está claro que los modelos de 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 comprobar 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 de grado 2 no es significativa, por lo que podemos utilizar dos regresiones lineales para cada conjunto de datos.