1
$\begingroup$

I have monthly values of a continuous variable from many subjects, the mean of which on plotting show curvilinear pattern with lower values in summer months. How can I analyze and report the significance of this data with p-values? Thanks for your help.

A plot of example data is given below:

enter image description here

About 5 subjects per day for 3 years are checked (each subject checked only once and is not repeated) and the data is arranged with following column headings:

subjectID, month, variable_reading
1, 1, 15
2, 1, 10
3, 2, 8
...

The question is "Does the variable have a significant seasonal variation?"

$\endgroup$
21
  • $\begingroup$ Post your data and perhaps we can help. If your data is .confidential simply scale it. $\endgroup$ – IrishStat Nov 26 '14 at 2:05
  • 1
    $\begingroup$ "Data" isn't significant or not. What questions do you seek to answer with the data? $\endgroup$ – Glen_b Nov 26 '14 at 2:17
  • 1
    $\begingroup$ Although a cosmetic detail, try to persuade your software to show integers 1 to 12 as labels, not 2.5, 7.5, 12.5. In the absence of other information, this is a matter of seasonality. But the variability within months should also be unravelled and plotted against date, assuming that you have here the composite of several years' data. The nature of the data may guide appropriate modelling: if it's climatic or climatically driven, then a sinusoid may help; if it's social or economic, then incidence of holidays, vacations, etc. may be important. Half of time series analysis could be relevant.... $\endgroup$ – Nick Cox Nov 26 '14 at 2:45
  • 1
    $\begingroup$ As @Glen_b hinted, there is no significance assessment without a precise question to be answered. Several quite different models could be fitted here. By the way, a reduction to months could be throwing away useful detail if the data arrive more finely. $\endgroup$ – Nick Cox Nov 26 '14 at 2:51
  • 1
    $\begingroup$ A seasonal ARIMA model night be appropriate , a deterministic seasonal model with dummies might be appropriate , other models may also be of interest excluding models that use predictors such as time-squared, time-cubed etc. $\endgroup$ – IrishStat Nov 26 '14 at 2:58
3
$\begingroup$

I simulated data a bit like yours, three years of daily data, 5 samples per day.

set.seed(42)
ts  <- seq(as.POSIXct('2012/1/1 00:00'), as.POSIXct('2014/12/31 23:59'), 'days')
ts  <- rep(ts, each=5) ## 5 samples every day
n   <- length(ts)

and generated a response by adding noise to a - b * cos(t):

y <- 7.5 + -2 * cos(seq(n)*2*pi/(365*5)) + rnorm(n,0,1)

The data have roughly the range and mean of yours, with a strong seasonal signal:

enter image description here

Form a data frame based on integer values and generate a month index:

d <- data.frame(sid=seq(n), month=as.POSIXlt(ts)$mon+1, y=round(y))
head(d)

  sid month y
1   1     1 7
2   2     1 5
3   3     1 6
4   4     1 6
5   5     1 6
6   6     1 5

Plot the response by month:

boxplot(y~month, data=d, xlab='month', ylab='y', main='simulated data')

enter image description here

This is roughly similar to your own, with a seasonal signal both stronger and clearer.

Fit a very simple linear model with separate means for each month:

fit <- lm(y ~ factor(month), data=d)
summary(fit)

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)      5.57419    0.04957 112.453  < 2e-16 ***
factor(month)2   0.47287    0.07173   6.592 4.74e-11 ***
factor(month)3   1.37849    0.07010  19.664  < 2e-16 ***
factor(month)4   2.41935    0.07010  34.512  < 2e-16 ***
factor(month)5   3.35054    0.07010  47.796  < 2e-16 ***
factor(month)6   3.86358    0.07068  54.661  < 2e-16 ***
factor(month)7   3.83871    0.07010  54.759  < 2e-16 ***
factor(month)8   3.23226    0.07010  46.108  < 2e-16 ***
factor(month)9   2.34358    0.07068  33.156  < 2e-16 ***
factor(month)10  1.35692    0.07068  19.197  < 2e-16 ***
factor(month)11  0.58803    0.07068   8.319  < 2e-16 ***
factor(month)12 -0.01075    0.07010  -0.153    0.878
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.069 on 5468 degrees of freedom
Multiple R-squared:  0.633, Adjusted R-squared:  0.6323
F-statistic: 857.4 on 11 and 5468 DF,  p-value: < 2.2e-16

The output here assumes January as a reference month and bundles the January mean with the intercept term. The significance of the individual month factors are then interpreted as "is the mean for this month significantly different to January?". For these data, all months are significantly different to January, except December. December and January are (statistically) indistinguishable from one another. One interpretation of this is that there is a seasonal difference between (southern hemisphere) summer months and the rest of the year.

This is not a great model by any stretch of the imagination. It doesn't even do a spectacular job at estimating the true simulated mean of 7.5. There are many shortcomings and criticsms that could be made. Importantly, it doesn't account for correlations in response between adjacent months, for example. It does, however, satisfy your initial criteria of answering the question about significant time-varying structure and reporting p-values. It's crude but "good enough" to detect a seasonal difference (in these data, even after a lot of detail has been lost to month binning). Any refinement to the correlation structure will only enhance it's sensitivity.

The point is not to defend this particular model but to say: any model you fit will have shortcomings. Finding an appropriate model for your needs depends on you having clear ideas about the goals of the analysis and the context in which it will be evaluated.

As @Nick Cox, @IrishStat (and others) have pointed out, there are many things that could be done here, almost all of them more sophisticated than this approach. But it is difficult to know what to suggest without a clear idea of what you aim to achieve. The question "Does the variable have a significant seasonal variation?" is not focused enough to guide the choice of model framework. In some contexts sARIMA may be warranted, in others the above may be adequate.

I hope that helps.

$\endgroup$
3
  • $\begingroup$ Thanks for an excellent, well illustrated answer. It is of great help to me. $\endgroup$ – rnso Nov 26 '14 at 13:47
  • $\begingroup$ @goangit I made the distinction that often the seasonal component could be deterministic as you suggest. AUTOBOX actually tests for both kinds and furthermore will detect significant change points in the 11 seasonal dummies. From my post : "HOWEVER one should also consider models that have 11 deterministic dummies as a possible alternative." $\endgroup$ – IrishStat Nov 26 '14 at 15:52
  • $\begingroup$ What I mean is that no matter what software you are using you should consider both possible kinds of seasonality and definitely check on changes in parameters over time and of course changes in the error variance over time AND deliver remedies. $\endgroup$ – IrishStat Nov 26 '14 at 16:07
1
$\begingroup$

A seasonal ARIMA model is identified via the auto-correlation function. The kth order auto-correlation is essentially a regression coefficient . The k regression coefficients are each based upon a simple two variate model ( y versus x ) reflecting a lag effect of order k in x. A classic discussion (showing the unnecessary log transformation) of the very seasonal airline series may help http://www.autobox.com/pdfs/vegas_ibf_09a.pdf starting at slide 14.

If the 12th auto-correlation and the 24th auto-correlation etc. appear to haVe a decaying pattern this is often indicative of a seasonal auto-regressive process HOWEVER one should also consider models that have 11 deterministic dummies as a possible alternative. HTH.

$\endgroup$
2
  • $\begingroup$ Just to spell it out: @IrishStat's examples assume that you are analysing monthly data. $\endgroup$ – Nick Cox Nov 26 '14 at 13:02
  • $\begingroup$ @NickCox THanks Nick I should have said 1s , 2s ,etc rather than 12 and 24 where s= frequency of the data $\endgroup$ – IrishStat Nov 26 '14 at 14:57

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.