Comparer la signification statistique de la différence entre deux régressions polynomiales dans R

Donc, tout d'abord, j'ai fait des recherches sur ce forum, et je sais que des questions extrêmement similaires ont été posées, mais elles n'ont généralement pas été répondues correctement ou parfois la réponse n'est tout simplement pas assez détaillée pour que je puisse comprendre. Donc cette fois, ma question est: j'ai deux ensembles de données, sur chacun, je fais une régression polynomiale comme ceci:

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

Les tracés de régressions polynomiales sont:

Les coefficients sont:

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

Et maintenant, je veux savoir, s'il existe un moyen d'utiliser une fonction R pour effectuer un test qui me dirait s'il y a ou non une signification statistique dans la différence entre les deux régressions polynomiales sachant que l'intervalle de jours pertinent est [ 1,100].

D'après ce que j'ai compris, je ne peux pas appliquer directement le test anova car les valeurs proviennent de deux ensembles de données différents ni de l'AIC, qui est utilisé pour comparer le modèle / les données réelles.

J'ai essayé de suivre les instructions données par @Roland dans la question connexe mais j'ai probablement mal compris quelque chose en regardant mes résultats:

Voici ce que j'ai fait :

J'ai combiné mes deux ensembles de données en un seul.

fest le facteur variable dont @Roland a parlé. J'ai mis 1s pour le premier set et 0s pour l'autre.

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

Mes données ressemblent maintenant à ceci:

Le rouge est fit3ANOVAqui fonctionne toujours, mais j'ai un problème avec le bleu, fit3ANOVACNle modèle a des résultats étranges. Je ne sais pas si le modèle d'ajustement est correct, je ne comprends pas exactement ce que @Roland voulait dire.

En considérant la solution @DeltaIV, je suppose que dans ce cas: les modèles sont significativement différents même s'ils se chevauchent. Ai-je raison de le supposer?

#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

Comme vous le voyez, fit1est nettement mieux que fit0, c'est-à-dire que l'effet de la variable de regroupement est significatif. Étant donné que la variable de regroupement représente les ensembles de données respectifs, les ajustements polynomiaux aux deux ensembles de données peuvent être considérés comme sensiblement différents.

La réponse de @Ronald est la meilleure et elle est largement applicable à de nombreux problèmes similaires (par exemple, existe-t-il une différence statistiquement significative entre les hommes et les femmes dans la relation entre le poids et l'âge?). Cependant, j'ajouterai une autre solution qui, bien qu'elle ne soit pas aussi quantitative (elle ne fournit pas de valeur p ), donne un bel affichage graphique de la différence.

EDIT : selon cette question , il semble que predict.lmla fonction utilisée par ggplot2pour calculer les intervalles de confiance, ne calcule pas les bandes de confiance simultanées autour de la courbe de régression, mais uniquement les bandes de confiance ponctuelles. Ces dernières bandes ne sont pas les bonnes pour évaluer si deux modèles linéaires ajustés sont statistiquement différents, ou disent d'une autre manière, s'ils pourraient être compatibles avec le même vrai modèle ou non. Ce ne sont donc pas les bonnes courbes pour répondre à votre question. Puisqu'apparemment il n'y a pas de R intégré pour obtenir des bandes de confiance simultanées (étrange!), J'ai écrit ma propre fonction. C'est ici:

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)

Les bandes internes sont celles calculées par défaut par geom_smooth: ce sont des bandes de confiance ponctuelles à 95% autour des courbes de régression. Les bandes externes semi-transparentes (merci pour la pointe graphique, @Roland) sont plutôt les bandes de confiance simultanées à 95%. Comme vous pouvez le voir, elles sont plus grandes que les bandes ponctuelles, comme prévu. Le fait que les bandes de confiance simultanées des deux courbes ne se chevauchent pas peut être considéré comme une indication du fait que la différence entre les deux modèles est statistiquement significative.

Bien sûr, pour un test d'hypothèse avec une valeur p valide , l'approche @Roland doit être suivie, mais cette approche graphique peut être considérée comme une analyse exploratoire des données. De plus, l'intrigue peut nous donner quelques idées supplémentaires. Il est clair que les modèles des deux ensembles de données sont statistiquement différents. Mais il semble également que deux modèles de degré 1 correspondraient presque aussi bien aux données qu'aux deux modèles quadratiques. On peut facilement tester cette hypothèse:

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 différence entre le modèle de degré 1 et le modèle de degré 2 n'est pas significative, nous pouvons donc aussi utiliser deux régressions linéaires pour chaque ensemble de données.