10
$\begingroup$

So first of all I did some research on this forum, and I know extremely similar questions have been asked but they usually haven't been answered properly or sometimes the answer are simply not detailed enough for me to understand. So this time my question is : I have two sets of data, on each, I do a polynomial regression like so :

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

The polynomial regressions plots are:

enter image description here

The coefficients are :

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

And now I want to know, if there is a way to use an R function to do a test that would tell me whether or not there is a statistical significance in the difference between the two polynomials regression knowing that the relevant interval of days is [1,100].

From what I understood I can not apply directly the anova test because the values come from two different sets of data nor the AIC, which is used to compare model/true data.

I tried to follow the instructions given by @Roland in the related question but I probably misunderstood something when looking at my results :

Here is what I did :

I combined both my datasets into one.

f is the variable factor that @Roland talked about. I put 1s for the first set and 0s for the other one.

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

My data looks like this now :

enter image description here

The red one is fit3ANOVA which is still working but I have a problem with the blue one fit3ANOVACN the model has weird results. I don't know if the fit model is correct, I do not understand what @Roland meant exactly.

Considering @DeltaIV solution I suppose that in that case : enter image description here The models are significantly different even though they overlap. Am I right to assume so ?

$\endgroup$
  • $\begingroup$ It seems to me that user @Roland's comment to the question you are linking to, answers your question perfectly. What is exactly that you don't understand? $\endgroup$ – DeltaIV Aug 22 '16 at 10:29
  • $\begingroup$ Well a couple of things, I was not sure wether or not this was a proper answer since it was in the comment section, but if it's working then, I just need to be sure I understood. Basically, I should create a new dataset where I create a column with like 1s and 0s depending on which datasets they originally came from ? Then after that I create two models one with every data the other one with only one of the datasets taken into account. Then I apply the anova test. Is that it ? Also I've never used the anova test, I saw they talked about proper p-value what would that be exactly ? $\endgroup$ – PaoloH Aug 22 '16 at 12:13
  • 1
    $\begingroup$ In your case the two regressions are on the same interval . This is the best case, to interpret confidence bands for linear regression. In this case, the two regressions are not statistically different if they are completely inside each other's confidence band over the whole interval ($[0,100]$ in your case)- not if they just overlap in a small subinterval. $\endgroup$ – DeltaIV Aug 23 '16 at 19:24
14
$\begingroup$
#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

As you see, fit1 is significantly better than fit0, i.e. the effect of the grouping variable is significant. Since the grouping variable represents the respective datasets, the polynomial fits to the two datasets can be considered significantly different.

$\endgroup$
  • $\begingroup$ I'm sorry this must be obvious but i'm not familiar with the Anova test results what is telling us that fit1 is better than fit0 ? Is it the Pr(>F) which is extremely low ? $\endgroup$ – PaoloH Aug 22 '16 at 13:55
  • 1
    $\begingroup$ The p-value tells you if the models are significantly different (lower p-value means "more different" taking into account the variation, usually p < 0.05 is considered significant). The smaller RSS indicates the better fitting model. $\endgroup$ – Roland Aug 22 '16 at 13:57
  • $\begingroup$ @PaoloH Btw., you should avoid ratios as dependent variables. Assumptions of ordinary least squares models don't hold with such a dependent variable. $\endgroup$ – Roland Aug 22 '16 at 14:49
8
$\begingroup$

@Ronald 's answer is the best and it's widely applicable to many similar problems (for example, is there a statistically significant difference between men and women in the relationship between weight and age?). However, I'll add another solution which, while not as quantitative (it doesn't provide a p-value), gives a nice graphical display of the difference.

EDIT: according to this question, it looks like predict.lm, the function used by ggplot2 to compute the confidence intervals, doesn't compute simultaneous confidence bands around the regression curve, but only pointwise confidence bands. These last bands are not the right ones to assess if two fitted linear models are statistically different, or said in another way, whether they could be compatible with the same true model or not. Thus, they are not the right curves to answer your question. Since apparently there's no R builtin to get simultaneous confidence bands (strange!), I wrote my own function. Here it is:

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)

enter image description here

The inner bands are those computed by default by geom_smooth: these are pointwise 95% confidence bands around the regression curves. The outer, semitransparent bands (thanks for the graphics tip, @Roland ) are instead the simultaneous 95% confidence bands. As you can see, they're larger than the pointwise bands, as expected. The fact that the simultaneous confidence bands from the two curves don't overlap can be taken as an indication of the fact that the difference between the two models is statistically significant.

Of course, for a hypothesis test with a valid p-value, @Roland approach must be followed, but this graphical approach can be viewed as exploratory data analysis. Also, the plot can give us some additional ideas. It's clear that the models for the two data set are statistically different. But it also looks like two degree 1 models would fit the data nearly as well as the two quadratic models. We can easily test this hypothesis:

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

The difference between the degree 1 model and the degree 2 model is not significant, thus we may as well use two linear regressions for each data set.

$\endgroup$
  • 2
    $\begingroup$ +1 for plotting. A crucial part of statistical analyses. $\endgroup$ – Roland Aug 22 '16 at 19:49
  • $\begingroup$ Just to make sure, on your method : the fact that the "confidence intervals from the two curves don't overlap can be taken as an indication of the fact that the difference between the two models is statistically significant." But I can't say that the difference isn't significant if they do overlap right ? $\endgroup$ – PaoloH Aug 23 '16 at 8:07
  • $\begingroup$ To be more specific I added an example in the post. $\endgroup$ – PaoloH Aug 23 '16 at 8:38
  • $\begingroup$ @PaoloH, since you added a new plot in your question, I'll add a comment there. $\endgroup$ – DeltaIV Aug 23 '16 at 19:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.