4

I often am in the situation of having data sets consisting of an independent variable, a dependent variable, and a factor with multiple levels - for instance, calibration curves for an instrument measured on different days. I would like to know whether the data are best described by a single fit line or by a different fit line for each factor level; i.e., whether the data are best described by a single calibration curve based on all of the data, or whether there are significant differences in the individual curves for each day.

From what I understand, ANCOVA will tell me separately whether the factor interacts with the slope and then whether it interacts with the intercept. What I want to know is whether the factor has a significant effect on the slope and the intercept of the line.

An example in R:

require(reshape)  # for melt()
#Set up some data
set.seed(0)
x <- seq(from=0, to=2, length.out=50)
y1 <- x + rnorm(length(x)) + 0.3
y2 <- 1.2*x + rnorm(length(x)) 

d <- data.frame(x=x, day1=y1, day2=y2)#, day3=y3, day4=y4)

dm <- melt(d, id.vars="x")

ggplot(dm, aes(x=x, y=value, colour=variable)) + geom_point() + geom_smooth(method="lm", se=T)

enter image description here

m <- lm(value ~ x*variable, data=dm)
summary(m)

enter image description here

Here, the effect of the factor and the interaction of the factor with the independent variable are each marginally significant - at alpha = 0.05 we would fail to reject the hypothesis that the factor has an effect on slope or intercept. However, taken together, perhaps they matter. Is there a good way to assess this?

Stephan Kolassa
  • 95,027
  • 13
  • 197
  • 357
Drew Steen
  • 674
  • 6
  • 18

2 Answers2

3

What you need is model simplification. You can use stepwise deletion method, which removes term by term, starting from the interactions of highest order. Or you can compare AIC of all possible models and select the model with lowest AIC. This is the easiest approach in this case.

m1 <- lm(value ~ x*variable, data=dm)    # your original model
m2 <- lm(value ~ 0 + variable + x:variable, data=dm) # the same model, parametrized
        # in more sane way - shows the coefficients you used for computation

m3 <- lm(value ~ 0 + variable + x, data=dm) # simplification - global slope
m4 <- lm(value ~ x, data=dm)   # further simplification - global intercept
m5 <- lm(value ~ variable, data=dm) # another variant - no slope, only categories

Let's compare AIC:

> AIC(m1)
[1] 265.857
> AIC(m2)
[1] 265.857
> AIC(m3)
[1] 267.6916
> AIC(m4)
[1] 266.0254
> AIC(m5)
[1] 313.5748

This shows that the first model is the best one - so you actually cannot use global slope and global intercept, you should use per-category slope and intercept.

Tomas
  • 5,735
  • 11
  • 52
  • 93
  • +1 for choosing models using AIC (an extremely good background monograph on AIC is this: http://www.springer.com/statistics/statistical+theory+and+methods/book/978-0-387-95364-9). However, beware of stepwise procedures. They are fine if your end goal is prediction for new data, but they cannot be used if the end goal is inference (since they filter for p values, the p values for the final model will be too small, since the filtering invalidates the assumptions on calculating those p values). – Stephan Kolassa Jan 18 '13 at 08:28
  • @StephanKolassa, thanks for your comment, but I don't get the note about stepwise procedures. What do you mean with filtering? *"... p values for the final model will be too small"* - you mean too big? Because by simplification the p-values (of parameters and whole regression) will usualy grow, as the model will fit less.. – Tomas Jan 18 '13 at 09:43
  • Stepwise procedures usually remove "insignificant" predictors or add "significant" predictors by checking p-values; this is what I meant by "filtering". So if you determine your model by removing everything where the p value is large, you will artificially bias the p values of the remaining terms downward. See http://www.stata.com/support/faqs/statistics/stepwise-regression-problems/ as well as http://onlinelibrary.wiley.com/doi/10.1111/j.1365-2656.2006.01141.x/pdf – Stephan Kolassa Jan 18 '13 at 09:50
  • @StephanKolassa, "biasing p-values" is a new concept to me. Biasing compared to what? There is no "true p-value" as far as I know, or is it? We usualy want to see most parsimonious model and look what's significant, that's all. – Tomas Jan 18 '13 at 09:56
  • Stepwise methods will bias the final coefficient estimates compared to the estimates from a non-"optimized" model (see links above); I think one can call the effect on p-values "biasing". And yes, we want a parsimonious model, but if we use stepwise methods, we need to account for that ("data snooping") in assessing significance in the final model, and that is non-trivial. – Stephan Kolassa Jan 18 '13 at 09:59
  • @StephanKolassa, but what do you call non-"optimized" model? Which one? I keep comparing tand referring to some "canonical" model, which is "correct", but there is none, AFAIK... – Tomas Jan 18 '13 at 11:01
  • The non-"optimized" model would be the initial model, or whatever model makes sense based on theory. I'll happily agree that there is no "true" model. We can go into deep philosophical discussions here, but all this has been covered elsewhere, so let's just head over here: http://stats.stackexchange.com/questions/13686/what-are-modern-easily-used-alternatives-to-stepwise-regression or here: http://stats.stackexchange.com/questions/5360/stepwise-logistic-regression-and-sampling/38333#38333 – Stephan Kolassa Jan 18 '13 at 11:06
  • @StephanKolassa, thanks. In summary, is it believed that the p-values of the initial model are more "correct" than those of the most parsimonious model? – Tomas Jan 18 '13 at 11:11
  • 1
    The distributional assumptions underlying the calculation of p-values include "no model selection based on the data". Therefore, only the original model really fulfills those assumptions. (However, it may not fulfill normality requirements, which in turn we only can assess by looking at the data... which is part of why NHST is problematic.) – Stephan Kolassa Jan 18 '13 at 11:32
1

It may depend on what "matters" means to you. However, we can say that a model containing the factor and its interaction with the covariate explains significantly more variation than a model containing the covariate alone:

anova(update(m,.~variable),m)

Analysis of Variance Table

Model 1: value ~ variable
Model 2: value ~ x * variable
  Res.Df     RSS Df Sum of Sq     F    Pr(>F)    
1     98 126.854                                 
2     96  75.631  2    51.224 32.51 1.655e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

So, taken together, the factor and the covariate indeed seem to make a difference.

Stephan Kolassa
  • 95,027
  • 13
  • 197
  • 357
  • So should my general approach be to compare the model with the factor to the model without the factor, to see whether the former explains significantly more variance than the latter? (Also: is there an name for this process of comparing two models?) – Drew Steen Jan 17 '13 at 16:17
  • 1
    It may make sense to just throw the variable out and compare the model with the factor only to the model without anything. It makes little sense to, e.g., compare the model with the interaction to a model with the interaction *but without the main effect from the factor*. As I wrote, it depends on what "matters" means to you. You may want to look at AIC or cross-validation to compare models (which I would simply call "comparing models" ;-). – Stephan Kolassa Jan 17 '13 at 16:23
  • @DrewSteen Maybe I misunderstand your question but I would try Stephan's approache with `Model 1: value ~ x` instead of `Model 1: value ~ variable` – Stéphane Laurent Jan 17 '13 at 21:32