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.
Modele różnią się znacznie, mimo że się pokrywają. Czy mam rację zakładać, że tak?