Vergleichen Sie die statistische Signifikanz des Unterschieds zwischen zwei Polynomregressionen in R.

Also habe ich zuerst in diesem Forum recherchiert und ich weiß, dass extrem ähnliche Fragen gestellt wurden, aber sie wurden normalerweise nicht richtig beantwortet oder manchmal sind die Antworten einfach nicht detailliert genug, um von mir verstanden zu werden. Diesmal lautet meine Frage also: Ich habe zwei Datensätze, für jeden mache ich eine Polynomregression wie folgt:

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

Die Polynomregressionsdiagramme sind:

Die Koeffizienten sind:

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

Und jetzt möchte ich wissen, ob es eine Möglichkeit gibt, eine R-Funktion zu verwenden, um einen Test durchzuführen, der mir sagt, ob der Unterschied zwischen der Regression der beiden Polynome eine statistische Signifikanz aufweist oder nicht, wenn ich weiß, dass das relevante Intervall von Tagen [ist. 1.100].

Soweit ich verstanden habe, kann ich den Anova-Test nicht direkt anwenden, da die Werte aus zwei verschiedenen Datensätzen stammen, noch aus dem AIC, der zum Vergleichen von Modell- / True-Daten verwendet wird.

Ich habe versucht, den Anweisungen von @Roland in der entsprechenden Frage zu folgen, aber ich habe wahrscheinlich etwas falsch verstanden, als ich meine Ergebnisse betrachtet habe:

Folgendes habe ich getan:

Ich habe meine beiden Datensätze zu einem zusammengefasst.

fist der variable Faktor, über den @Roland gesprochen hat. Ich habe 1s für den ersten Satz und 0s für den anderen gesetzt.

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

Meine Daten sehen jetzt so aus:

Das rote ist fit3ANOVAdas, was noch funktioniert, aber ich habe ein Problem mit dem blauen, fit3ANOVACNdas Modell hat seltsame Ergebnisse. Ich weiß nicht, ob das Anpassungsmodell korrekt ist, ich verstehe nicht, was @Roland genau bedeutete.

In Anbetracht der @ DeltaIV-Lösung nehme ich in diesem Fall Folgendes an: Die Modelle unterscheiden sich erheblich, obwohl sie sich überlappen. Habe ich Recht, dies anzunehmen?

#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

Wie Sie sehen, fit1ist signifikant besser als fit0, dh der Effekt der Gruppierungsvariablen ist signifikant. Da die Gruppierungsvariable die jeweiligen Datensätze darstellt, kann die Polynomanpassung an die beiden Datensätze als signifikant unterschiedlich angesehen werden.

Die Antwort von @Ronald ist die beste und auf viele ähnliche Probleme anwendbar (gibt es beispielsweise einen statistisch signifikanten Unterschied zwischen Männern und Frauen in der Beziehung zwischen Gewicht und Alter?). Ich werde jedoch eine andere Lösung hinzufügen, die zwar nicht so quantitativ ist (keinen p- Wert liefert), aber eine schöne grafische Darstellung des Unterschieds liefert.

BEARBEITEN : Nach dieser Frage sieht es so aus predict.lm, als würde die Funktion ggplot2zur Berechnung der Konfidenzintervalle keine simultanen Konfidenzbänder um die Regressionskurve berechnen , sondern nur punktweise Konfidenzbänder. Diese letzten Bänder sind nicht die richtigen, um zu beurteilen, ob zwei angepasste lineare Modelle statistisch unterschiedlich sind oder auf andere Weise gesagt werden, ob sie mit demselben wahren Modell kompatibel sein könnten oder nicht. Daher sind sie nicht die richtigen Kurven, um Ihre Frage zu beantworten. Da anscheinend kein R eingebaut ist, um simultane Vertrauensbänder zu erhalten (seltsam!), Schrieb ich meine eigene Funktion. Hier ist es:

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)

Die inneren Bänder sind diejenigen, die standardmäßig berechnet werden durch geom_smooth: Dies sind punktweise 95% -Konfidenzbänder um die Regressionskurven. Die äußeren, halbtransparenten Bänder (danke für den Grafiktipp @Roland) sind stattdessen die gleichzeitigen 95% -Konfidenzbänder. Wie Sie sehen können, sind sie erwartungsgemäß größer als die punktweisen Bänder. Die Tatsache, dass sich die simultanen Konfidenzbänder aus den beiden Kurven nicht überlappen, kann als Hinweis darauf angesehen werden, dass der Unterschied zwischen den beiden Modellen statistisch signifikant ist.

Natürlich muss für einen Hypothesentest mit einem gültigen p- Wert der @ Roland-Ansatz befolgt werden, aber dieser grafische Ansatz kann als explorative Datenanalyse angesehen werden. Die Handlung kann uns auch einige zusätzliche Ideen geben. Es ist klar, dass die Modelle für die beiden Datensätze statistisch unterschiedlich sind. Es sieht aber auch so aus, als würden zwei Grad-1-Modelle fast genauso gut zu den Daten passen wie die beiden quadratischen Modelle. Wir können diese Hypothese leicht testen:

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

Der Unterschied zwischen dem Grad 1-Modell und dem Grad 2-Modell ist nicht signifikant, daher können wir auch zwei lineare Regressionen für jeden Datensatz verwenden.