I am doing statistical analysis of a natural experiment that consists of multiple years of measurements. I have two independent variables that are physically related to Y.

  1. I am interested in whether there are differences between years that are not explained by the model.
  2. I want to test whether a third predictor can reduce the differences between years.

I added Year as a categorial predictor in my linear regression (I am using natural cubic splines). My model looks like this:

Linear Regression Model

 ols(formula = depend ~ rcs(X1, 3) + rcs(X2, 4) + Yearf, data = data,
     x = T, y = T)

                 Model Likelihood    Discrimination
                       Ratio Test           Indexes
 Obs     869    LR chi2    805.76    R2       0.604
 sigma0.6327    d.f.           10    R2 adj   0.600
 d.f.    858    Pr(> chi2) 0.0000    g        0.848


      Min       1Q   Median       3Q      Max
 -2.75550 -0.40733  0.01893  0.42597  1.70108

            Coef    S.E.   t      Pr(>|t|)
 Intercept   1.5327 0.2067   7.41 <0.0001
 X1          1.0437 0.0525  19.89 <0.0001
 X1'        -0.8147 0.0686 -11.88 <0.0001
 X2          1.2507 0.1670   7.49 <0.0001
 X2'        -2.4775 0.6915  -3.58 0.0004
 X2''        3.2983 1.2123   2.72 0.0066
 Yearf=2016  0.2475 0.0814   3.04 0.0024
 Yearf=2017  0.1620 0.0802   2.02 0.0437
 Yearf=2018  0.0440 0.0862   0.51 0.6096
 Yearf=2019 -0.5260 0.0829  -6.34 <0.0001
 Yearf=2020  0.1457 0.0813   1.79 0.0734

Effect plots for the predictors look like this:

Enter image description here

After adding X3, the effect plot looks like this:

Enter image description here

Regarding my questions formulated at the beginning of this post, I would interpret the results as follows, but I am not sure if the whole approach is valid:

  1. If X1,X2 (and X3 in the second plot) are set to their average values, the mean response of Y for each year would be the value within the plot for Yearf.

  2. Visually it appears that the differences between the years reduces when X3 is added to the model while the standard error of the mean response increases. Is it valid to do this visually? As far as I know, an anova only tells me if the difference between year x and my base year (here 2015) is significant, but I don’t want to compare it for a specific year. I want to get an idea of the whole picture. Does this make sense to you?


1 Answer 1


Well done. I think your interpretation is OK. Running anova() will give you a more general test for a time trend than how you described it---a test of whether the time trend is flat, i.e., where there is a difference between any pair of years. It would also be useful to examine the difference in $R^2$ due to yearf with and without adjustment for x3.

The analysis would be slightly better were you to measure time with more resolution. Then you could model year + fraction of a year instead of integer year.

  • 1
    $\begingroup$ Thank you!! Using validate I get the following corrected R2: X1,X2 = 0.5377 X1,X2,Yearf = 0.5932 X1,X2,X3 = 0.6039 X1,X2,X3,Yearf = 0.6201 I'd conclude from this that X3 cannot explain all but some of the differences between the years. I have data with resolution of 30min-1 so I could go into a lot of more time-dependent details, but I also want to keep my readers interested. ;-) $\endgroup$
    – Felix Phl
    Jan 25 at 14:30

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.