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

10

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:

entrez la description de l'image ici

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:

entrez la description de l'image ici

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: entrez la description de l'image ici les modèles sont significativement différents même s'ils se chevauchent. Ai-je raison de le supposer?

PaoloH
la source
Il me semble que le commentaire de l'utilisateur @ Roland à la question à laquelle vous créez un lien, répond parfaitement à votre question. Qu'est-ce que tu ne comprends pas exactement?
DeltaIV
Eh bien, deux ou trois choses, je ne savais pas si c'était une bonne réponse car c'était dans la section des commentaires, mais si ça marche alors, je dois juste être sûr d'avoir bien compris. Fondamentalement, je devrais créer un nouvel ensemble de données où je crée une colonne avec des 1 et des 0 similaires, selon les ensembles de données dont ils sont originaires Ensuite, je crée deux modèles, l'un avec toutes les données, l'autre avec un seul des ensembles de données pris en compte. Ensuite, j'applique le test anova. Est-ce que c'est ça ? De plus, je n'ai jamais utilisé le test anova, j'ai vu qu'ils parlaient de la valeur p appropriée, qu'est-ce que ce serait exactement?
PaoloH
1
[0,100]

Réponses:

15
#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.

Roland
la source
Je suis désolé, cela doit être évident, mais je ne connais pas les résultats du test Anova qu'est-ce qui nous dit que fit1 est meilleur que fit0? Est-ce le Pr (> F) qui est extrêmement bas?
PaoloH
1
La valeur de p vous indique si les modèles sont significativement différents (une valeur de p inférieure signifie «plus différent» en tenant compte de la variation, généralement p <0,05 est considéré comme significatif). Le petit RSS indique le meilleur modèle d'ajustement.
Roland
@PaoloH Btw., Vous devez éviter les ratios en tant que variables dépendantes. Les hypothèses des modèles des moindres carrés ordinaires ne sont pas valables avec une telle variable dépendante.
Roland
8

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)

entrez la description de l'image ici

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.

DeltaIV
la source
3
+1 pour le traçage. Une partie cruciale des analyses statistiques.
Roland
Juste pour être sûr, sur votre méthode: le fait que "les intervalles de confiance des deux courbes ne se chevauchent pas peut être pris comme une indication du fait que la différence entre les deux modèles est statistiquement significative". Mais je ne peux pas dire que la différence n'est pas significative si elles se chevauchent, n'est-ce pas?
PaoloH
Pour être plus précis, j'ai ajouté un exemple dans le post.
PaoloH
@PaoloH, puisque vous avez ajouté une nouvelle intrigue dans votre question, j'y ajouterai un commentaire.
DeltaIV