Porównaj istotność statystyczną różnicy między dwiema regresjami wielomianowymi w R.

Najpierw więc przeprowadziłem badania na tym forum i wiem, że zadawano bardzo podobne pytania, ale zwykle nie udzielono na nie prawidłowej odpowiedzi lub czasami odpowiedzi nie są wystarczająco szczegółowe, aby je zrozumieć. Więc tym razem moje pytanie brzmi: mam dwa zestawy danych, na każdym wykonuję regresję wielomianową tak:

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

Wykresy regresji wielomianowej to:

Współczynniki są następujące:

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

A teraz chcę wiedzieć, czy istnieje sposób użycia funkcji R do wykonania testu, który powiedziałby mi, czy istnieje różnica statystyczna w różnicy między regresją dwóch wielomianów, wiedząc, że odpowiedni przedział dni wynosi [ 1100].

Z tego, co zrozumiałem, nie mogę bezpośrednio zastosować testu anova, ponieważ wartości pochodzą z dwóch różnych zestawów danych ani z AIC, który służy do porównywania danych modelowych / rzeczywistych.

Próbowałem postępować zgodnie z instrukcjami podanymi przez @Roland w powiązanym pytaniu, ale prawdopodobnie źle zrozumiałem, patrząc na moje wyniki:

Oto co zrobiłem:

Połączyłem oba zestawy danych w jeden.

fjest zmiennym czynnikiem, o którym mówił @Roland. Wstawiam 1 dla pierwszego zestawu i 0 dla drugiego.

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

Moje dane wyglądają teraz tak:

Ten czerwony fit3ANOVAwciąż działa, ale mam problem z niebieskim, fit3ANOVACNktóry ma dziwne wyniki. Nie wiem, czy model dopasowania jest poprawny, nie rozumiem dokładnie, co @Roland miał na myśli.

Biorąc pod uwagę rozwiązanie @DeltaIV, przypuszczam, że w takim przypadku: Modele różnią się znacznie, mimo że się pokrywają. Czy mam rację zakładać, że tak?

#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

Jak widać, fit1jest znacznie lepszy niż fit0, tj. Wpływ zmiennej grupującej jest znaczący. Ponieważ zmienna grupująca reprezentuje odpowiednie zestawy danych, dopasowanie wielomianowe do dwóch zestawów danych można uznać za znacznie różne.

Odpowiedź @Ronald jest najlepsza i ma szerokie zastosowanie do wielu podobnych problemów (na przykład, czy istnieje statystycznie istotna różnica między mężczyznami i kobietami w relacji między wagą a wiekiem?). Dodam jednak inne rozwiązanie, które choć nie jest tak ilościowe (nie zapewnia wartości p ), daje ładny graficzny obraz różnicy.

EDYCJA : zgodnie z tym pytaniem wygląda na to predict.lm, że funkcja używana ggplot2do obliczania przedziałów ufności nie oblicza równoczesnych pasm ufności wokół krzywej regresji, a jedynie pasma pewności punktowej. Te ostatnie pasma nie są odpowiednie do oceny, czy dwa dopasowane modele liniowe są statystycznie różne lub, mówiąc inaczej, czy mogą być kompatybilne z tym samym prawdziwym modelem, czy nie. Dlatego nie są to odpowiednie krzywe, aby odpowiedzieć na twoje pytanie. Ponieważ najwyraźniej nie ma wbudowanego R, aby uzyskać jednoczesne przedziały ufności (dziwne!), Napisałem własną funkcję. Oto on:

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)

Wewnętrzne pasma to te obliczane domyślnie przez geom_smooth: są to punktowe 95% przedziały ufności wokół krzywych regresji. Zewnętrzne, półprzezroczyste pasma (dzięki za końcówkę graficzną, @Roland) są równoczesnymi 95% pasmami pewności. Jak widać, są one większe niż pasma punktowe, zgodnie z oczekiwaniami. Fakt, że równoczesne pasma ufności z dwóch krzywych się nie pokrywają, można uznać za wskazówkę, że różnica między dwoma modelami jest statystycznie znacząca.

Oczywiście w przypadku testu hipotez z prawidłową wartością p należy zastosować podejście @Roland, ale to podejście graficzne można postrzegać jako analizę danych eksploracyjnych. Fabuła może dać nam dodatkowe pomysły. Oczywiste jest, że modele dwóch zestawów danych są statystycznie różne. Wygląda jednak na to, że dwa modele 1 stopnia pasowałyby do danych prawie tak dobrze, jak dwa modele kwadratowe. Możemy łatwo przetestować tę hipotezę:

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

Różnica między modelem stopnia 1 a modelem stopnia 2 nie jest znacząca, dlatego równie dobrze możemy zastosować dwie regresje liniowe dla każdego zestawu danych.