Bandingkan signifikansi statistik dari perbedaan antara dua regresi polinomial dalam R

Jadi pertama-tama saya melakukan riset di forum ini, dan saya tahu pertanyaan yang sangat mirip telah ditanyakan tetapi biasanya mereka belum dijawab dengan benar atau kadang-kadang jawabannya tidak cukup detail untuk saya pahami. Jadi kali ini pertanyaan saya adalah: Saya punya dua set data, masing-masing, saya melakukan regresi polinomial seperti:

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

Plot regresi polinomial adalah:

Koefisiennya adalah:

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

Dan sekarang saya ingin tahu, apakah ada cara untuk menggunakan fungsi R untuk melakukan tes yang akan memberi tahu saya apakah ada signifikansi statistik dalam perbedaan antara dua regresi polinomial mengetahui bahwa interval hari yang relevan adalah [ 1.100].

Dari apa yang saya pahami saya tidak bisa langsung menerapkan tes anova karena nilai-nilai berasal dari dua set data yang berbeda atau AIC, yang digunakan untuk membandingkan model / data yang sebenarnya.

Saya mencoba mengikuti instruksi yang diberikan oleh @Roland dalam pertanyaan terkait tetapi saya mungkin salah mengerti sesuatu ketika melihat hasil saya:

Inilah yang saya lakukan:

Saya menggabungkan kedua set data saya menjadi satu.

fadalah faktor variabel yang dibicarakan @Roland. Saya menempatkan 1 untuk set pertama dan 0 untuk yang lainnya.

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

Data saya terlihat seperti ini sekarang:

Yang merah fit3ANOVAmasih berfungsi tapi saya punya masalah dengan yang biru fit3ANOVACNmodelnya punya hasil yang aneh. Saya tidak tahu apakah model yang cocok sudah benar, saya tidak mengerti apa arti sebenarnya dari @Rand.

Mempertimbangkan solusi @DeltaIV, saya kira dalam kasus itu: Modelnya sangat berbeda walaupun tumpang tindih. Apakah saya benar berasumsi demikian?

#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

Seperti yang Anda lihat, fit1secara signifikan lebih baik daripada fit0, yaitu pengaruh variabel pengelompokan adalah signifikan. Karena variabel pengelompokan mewakili masing-masing dataset, kecocokan polinomial dengan dua dataset dapat dianggap berbeda secara signifikan.

Jawaban @Ronald adalah yang terbaik dan secara luas dapat diterapkan untuk banyak masalah serupa (misalnya, apakah ada perbedaan yang signifikan secara statistik antara pria dan wanita dalam hubungan antara berat badan dan usia?). Namun, saya akan menambahkan solusi lain yang, meskipun bukan sebagai kuantitatif (tidak memberikan nilai- p ), memberikan tampilan grafis yang bagus tentang perbedaannya.

EDIT : menurut pertanyaan ini , sepertinya predict.lm, fungsi yang digunakan oleh ggplot2untuk menghitung interval kepercayaan, tidak menghitung band - band kepercayaan simultan di sekitar kurva regresi, tetapi hanya band-band percaya diri yang searah. Pita terakhir ini bukan yang tepat untuk menilai apakah dua model linier yang sesuai secara statistik berbeda, atau dikatakan dengan cara lain, apakah mereka dapat kompatibel dengan model yang sama atau tidak. Jadi, itu bukan kurva yang tepat untuk menjawab pertanyaan Anda. Karena ternyata tidak ada R builtin untuk mendapatkan band kepercayaan simultan (aneh!), Saya menulis fungsi saya sendiri. Ini dia:

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)

Pita dalam adalah pita yang dihitung secara default dengan geom_smooth: ini adalah pita kepercayaan 95% pointwise di sekitar kurva regresi. Pita terluar, semitransparan (terima kasih atas tip grafisnya, @Roland) sebagai gantinya band kepercayaan 95% simultan . Seperti yang Anda lihat, mereka lebih besar dari band-band pointwise, seperti yang diharapkan. Fakta bahwa pita kepercayaan simultan dari dua kurva tidak tumpang tindih dapat diambil sebagai indikasi fakta bahwa perbedaan antara kedua model secara statistik signifikan.

Tentu saja, untuk uji hipotesis dengan nilai- p yang valid , pendekatan @Rand harus diikuti, tetapi pendekatan grafis ini dapat dilihat sebagai analisis data eksplorasi. Juga, alur ceritanya dapat memberi kita beberapa ide tambahan. Jelas bahwa model untuk dua set data berbeda secara statistik. Tetapi juga terlihat bahwa dua model derajat 1 akan cocok dengan data hampir sama halnya dengan dua model kuadratik. Kami dapat dengan mudah menguji hipotesis ini:

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

Perbedaan antara model derajat 1 dan model derajat 2 tidak signifikan, sehingga kami juga dapat menggunakan dua regresi linier untuk setiap set data.