A Method of Trigonometric Modelling of Seasonal Variation Demonstrated with Multiple Sclerosis Relapse Data

Summary

Combining plot analysis with trigonometric regression is a robust method for exploring complex, cyclical phenomena such as relapse onset timing in multiple sclerosis (MS). This method enabled unbiased characterisation of seasonal trends in relapse onset permitting novel inferences around the influence of seasonal variation, ultraviolet radiation (UVR) and latitude.

Abstract

This report describes a novel Stata-based application of trigonometric regression modelling to 55 years of multiple sclerosis relapse data from 46 clinical centers across 20 countries located in both hemispheres. Central to the success of this method was the strategic use of plot analysis to guide and corroborate the statistical regression modelling. Initial plot analysis was necessary for establishing realistic hypotheses regarding the presence and structural form of seasonal and latitudinal influences on relapse probability and then testing the performance of the resultant models. Trigonometric regression was then necessary to quantify these relationships, adjust for important confounders and provide a measure of certainty as to how plausible these associations were. Synchronization of graphing techniques with regression modelling permitted a systematic refinement of models until best-fit convergence was achieved, enabling novel inferences to be made regarding the independent influence of both season and latitude in predicting relapse onset timing in MS. These methods have the potential for application across other complex disease and epidemiological phenomena suspected or known to vary systematically with season and/or geographic location.

Introduction

The most common form of multiple sclerosis (MS) is Relapsing Remitting Multiple Sclerosis (RRMS). RRMS is characterized by episodic deteriorations in neurological function, followed by partial or complete recovery. Globally, the incidence and prevalence of MS increase with increasing distance away from the equator in both hemispheres.^1-3 Whether the frequency of relapse events that occur specifically in RRMS also vary with latitude, and whether there is any underlying seasonal variation in any such association, remains unclear. To date studies exploring seasonality in relapse timing have been limited to single clinical centers, limiting any inferences regarding seasonal trends in relapse timing to solitary geographical locations and thus unable to explore broader latitudinal influences. ^4-14 These studies have been further limited by small sample sizes and sparse relapse data. A 2000 meta-analysis of ten studies from clinical centres in Europe, the United States and Canada, where each study included a minimum of thirty cases reporting the season-of-onset of relapses, described a clear seasonal trend in the timing of relapse onset, with relapses peaking in spring and with a winter trough⁴. Similar cyclical annual trends in onset have been observed in subsequent, albeit smaller, studies in both Japan¹⁵ and Spain¹⁶. However, a comparable United States study failed to corroborate this pattern¹⁷. To date, these studies and observations have been limited to the northern hemisphere. The MSBase study group recently analyzed a large global dataset of MS relapses across both northern and southern hemispheres to explore seasonal trends in the timing of relapse onset in addition to latitudinal influences on the relationship between peak relapse probability and seasonal ultraviolet radiation (UVR) trough¹⁸. Central to these methods was the application of trigonometric regression to visualize and evaluate trends in the timing of relapse onset and UVR distributions.

The overall goal of this study was to test the hypothesis that temporal variation in the timing of relapse onset in MS varied predictably with season in both the northern and southern hemispheres and this seasonality was influenced by latitude. The rationale for the use of trigonometric modelling to investigate these questions was its flexibility for characterizing two- or three-dimensional phenomena that are known or suspected to describe discrete, predictable and consistent shapes or patterns, such as the annual cycle of peaks and troughs commonly observed in biological or epidemiological phenomena possessing seasonality.^19-22 A disadvantage of conventional time-series analyses, including Fourier analysis, is the presumption that time series are often characterized by stochastic processes.^21,23,24 By contrast, incorporating trigonometric functions into a regression type model has the advantage of both facilitating exploration of regular and systematic structures in periodic data whilst exploiting the regression model structure to explore other correlates or adjust for confounders of seasonality.

Trigonometric regression has previously been used widely in the medical epidemiological literature to explore temporality in topics as diverse infectious disease outbreak detection, the role of circadian rhythms in everything from autonomic nervous system dysfunction to preterm placental abruption through to seasonal correlates of congenital malformations and the timing of presentations of accident and emergency.^25-32 Such modelling typically demands larger sample sizes than more conventional time-series analyses and as such this is the first time it has been applied to a global dataset of MS relapse onset. Trigonometric regression as described here is suitable tool for investigators exploring any phenomena which is known to or suspected of cycling systematically over time. Not only can such modelling help characterise and visualize these patterns, it further permits the user to explore potential drivers and correlates of these trends.

Regarding the specific example of MS relapse onset presented here, the use of scatter and residual plots to visualize and assess how closely a hypothesized trigonometric model form fits the data constitutes the critical step in determining: 1) whether the observed data provide sufficient evidence to support a hypothesis of seasonality or other temporal trends in the timing of relapse onset; and 2) whether the frequency and arrangement of sine and cosine functions which define a particular trigonometric model is adequate to permit use of this model for subsequent inference and prediction. Regression modelling also permits control for important confounders of any observed seasonal or latitudinal effect such as patient-level propensities for relapse, particularly factors which in themselves are time-varying such as the duration of pre-relapse exposure to disease-modifying drug (DMD) treatment. Isolating independent geographic and temporal predictors and correlates of relapse onset timing in MS has the potential to guide biological investigation into the mechanisms of relapse events which in turn may inform the development of future treatment interventions aimed at preventing or delaying disease exacerbation.

The MSBase Registry

MS patients contributing relapse data to this analysis were sourced from the international MSBase registry. Established in 2004, the registry longitudinally collates demographic, disease activity, clinical examination and investigation characteristics and metrics from consenting patients attending MS clinic using an internet-based, physician-owned and operated system.³³ Member centers follow a common protocol that defines the minimum dataset required to be uploaded at agreed regular intervals to ensure outcome data such as relapse events are consistently and prospectively compiled. The date of relapse onset is included as a mandatory minimum dataset variable. In addition relevant clinical data associated with these relapse events is commonly collected including corticosteroid treatment and functional system affected. The use of the common iMed data entry system further ensures a unified approach across centers to data collection and reporting. This project holds Human Research Ethics Committee approval or exemption at each contributing center. Informed consent according to local laws from all patients included in the analysis is mandatory.

Inclusion criteria

A total of 9811 patients contributing 32,762 relapse events were included in the analysis. Clinical MS centres with a minimum of 20 registered patients consented, uploaded and tracked in the registry as of the 1^st December 2013 (date of data compilation) were eligible for inclusion in the analysis. To ensure all relapse events included in the analysis were prospectively observed, only relapse onsets dated subsequent to the first recorded patient disability assessment (using the Kurtzke Expanded Disability Status Score (EDSS)) were included in the analysis. All patients contributing relapse data to the analysis satisfied formal diagnostic criteria for MS.^34,35

Outcome measures

This study considered two primary outcomes: 1) whether there was temporal variation in the probability of relapse onset at the level of the geographic location, the hemisphere and/or globally; and 2) whether there was a relationship between latitude and the lag, in months, between the timing of seasonal UVR trough and the subsequent peak relapse probability date. The MSBase Study group hypothesized that as absolute vitamin D levels are lower in regions further away from the equator and location-specific seasonal population level vitamin D nadirs are likely reached earlier following the winter solstice in such distal locations, then the effect of low vitamin D levels on increased MS relapse probability would similarly describe such temporal and latitudinal patterns.

Relapse definition and dates

A relapse was defined as occurrence of new symptoms or exacerbation of existing symptoms persisting for at least 24 hours, in the absence of concurrent illness or fever, and occurring at least 30 days after a previous attack. This definition has previously been applied in an MSBase relapse phenotype analysis. ³⁶ The follow-up period for each eligible patient across which relapse events could be observed was defined as the period spanning the date of first EDSS assessment through to the date of the most recent EDSS assessment recorded in the registry prior to the data of data extract and compilation. In instances where the exact day of relapse onset was unavailable or unable to be determined for a particular month, clinics used either the 1^st or 15^th day of the month as default dates. Of the 32,762 relapses analysed in this report, 7913 (24.2%) and 4594 (14.0%) were recorded on the 1^st and 15^th day of the month respectively, significantly higher than the proportions recorded on any other day of the month which ranged from 0.8% through 5.6%. To correct for this, relapses recorded on either the 1^st of 15^th day of the month were randomized to a day within a 15 day interval either side of both these default dates. The internal validity of this approach was confirmed via sensitivity analyses which demonstrated that the modelled estimate of peak relapse date under default date randomization was not significantly different from a model using either the original reported dates or excluding default dates entirely.

Protocol

NOTE: Each step described corresponds to a section of Stata code with the same number in the code file provided. Stata command names have been italicised in the following protocol.

1. Prepare and Plot the Observed Relapse Onset Data

1. Open a do file by clicking on the "New Do-file Editor" button and use the generate command to calculate the number of relapse onsets dated to each of the twelve calendar months for each of the three geographic levels to be modelled: location, hemisphere and global. Action command by clicking the "Execute (do)" do-file action button in the do-file.
2. Use the swilk orsktest command to test the underlying distribution of relapse counts for normality using a Shapiro-Wilk or modified Jarque-Bera test for location aggregate relapse data or disaggregate patient level relapse data respectively. Select code and click "Execute (do)".
   NOTE: In the presence of significant skew, apply a natural log transformation and subsequently test the log-transformed relapse count variable for approximate normality by reapplying the Shapiro-Wilk or modified Jarque-Bera test as appropriate.^38,39
3. Use the generate command to create a new variable "north_month" for southern hemisphere calendar months offset by +6 to allow plotting of both northern and southern hemisphere relapses by season along the same horizontal axis. Select code and click "Execute (do)".
   1. Graph a scatterplot of observed monthly relapse onsets with relapse frequency on the y-axis and calendar month on the x-axis for each hemisphere using the twoway scatter command. Repeat for each location. Observe pattern of peaks and troughs in relapse onset over the calendar year by viewing each plot in the graph viewer the automatically opens to screen.
4. Use the radar command to draw radar plots of the distribution of relapse frequency by calendar month with each radar axis capturing a single month ordered in a clockwise manner. Select code and click "Execute (do)".
   1. Repeat for all sites. Observe pattern of peaks and troughs in relapse onset over the calendar year by viewing each plot in the graph viewer the automatically opens to screen.
5. Run the seast command to apply an Edward's test of seasonality across the observed relapse data.^40-42 Repeat for all geographic levels.

2. Model Building and Selection

1. Use the generate command to specify the annual cycle sine and cosine trigonometric functions to be used in the regression. Select code and click "Execute (do)".
2. Use the regress command specify the form of the base model with relapse count as the dependent outcome variable and the sine and cosine terms calculated in step 2.1 as the primary explanatory variables.
   1. Add location-specific UVR³⁷ to the base model as an additional adjusting covariate and use the analytic weight aweight option to weight the model for the number of patients contributed by each location. Select code and click "Execute (do)".
      NOTE: Record the model coefficient of determination (R²) and the residual error in the results window that automatically opens to screen. Ultraviolet radiation: Daily average erythemally-weighted ambient UVR for each month from 1979 to 2004 inclusive was sourced from the National Aeronautics and Space Administration Earth Probe Total Ozone Mapping Spectrometer for all individual locations included in the analysis. ³⁷
3. Store the model predicted monthly log(relapse) using the predict command. Convert log relapses back to integer relapse counts by exponentiating the log(relapse) term using the generate command. Select code and click "Execute (do)". Repeat for all sites.
4. Overlay the exponentiated predicted monthly relapse estimates from 2.3 over the observed monthly relapse data using the twoway scatter command. Select code and click "Execute (do)".
   1. Repeat for all sites. View each plot in the graph viewer.
5. Use the regress command to expand the model specified in 2.2 by adding an additional harmonic sine/cosine pair. Select code and click "Execute (do)".
   NOTE: Record the residual error and the coefficient of determination. Save and transform model estimates as per 2.3 and plot model estimates over observed data as per 2.4. Repeat for all sites.
6. Use the regress command to further expand the model specified in 2.2 by adding two additional harmonic sine/cosine pairs. Select code and click "Execute (do)".
   NOTE: Record the residuals and the coefficient of determination. Compare this model directly with the base model using a likelihood ratio test. Use the estat ic post-estimation command to generate Akaike and Bayesian Information Criteria. Save and transform model estimates as per 2.3 and plot model estimates over observed data as per 2.4. Repeat for all geographic levels.

3. Estimating Peak Relapse Probability

1. Use the non-linear combination of estimators function (nlcom) to calculate the point estimate and 95% confidence interval for the phase-shift, using the best-fitting model identified from steps 2.1 through 2.6. Select code and click "Execute (do)".
   1. Convert these point estimates and associated confidence intervals to numbers representing calendar dates of peak relapse frequency (T_max) and trough relapse frequency (T_min) , where 1=1^st January and 365=31^st December and T_max = phase-shift + (365/4) and T_min = phase-shift + ((365/4)*3). Repeat for all geographic levels. Match T_max and T_min to a calendar date via the Excel look-up file.
2. Use the generate command to calculate peak-to-trough difference (T_max minus T_min) for each location, standardized for every 100 patients per site. Use a Wilcoxon rank-sum test to compare standardized peak-to-trough difference by latitude range. Select code and click "Execute (do)".

4. Modelling Ultraviolet Radiation Data

1. Run use command to load the UVR data. Calculate median monthly UVR for each location using the egen command. Select code and click "Execute (do)".
2. Graph a scatterplot of monthly UVR (y-axis) by calendar month (x-axis) for each location using the twoway scatter function. View each plot in the graph viewer the automatically opens to screen.
3. Repeat step 1.2 for the UVR data and use the regress command to specify a base model of location-level annual UVR trend where monthly UVR is specified as the dependent outcome variables and the sine and cosine trigonometric functions specified in step 2.1 are incorporated into the model as the explanatory variables.
4. Repeat steps 2.4 through 2.6 for the UVR model and limited to the location-specific models only. This involves rerunning the twoway scatter command to overlay predicted estimates on observed data and using the regress command to run the expanded harmonic model alternatives.
5. Using the best-fitting model of location-specific monthly UVR identified in steps 4.2 through 4.4 use the generate command to calculate the phase-shift point estimate and associated 95% confidence interval for UVR by again applying the double-angle formulae specified in step 3.1. Calculate T_min (date of trough UVR) for each location using the formula described in step 3.1. Select code and click "Execute (do)".

5. Modelling UVR-trough-to-relapse-peak Lag

1. Append the model-estimated dates of seasonal UVR trough from step 4.5 and relapse peak dates from step 3.1 for each location using the merge command. Use the generate command to calculate the time lapsed in months between UVR trough date and subsequent relapse peak date. Select code and click "Execute (do)".
2. Use the sktest command to test the UVR-trough-to-relapse-peak lag variable for significant departures from normality using a Shapiro-Wilk test Select code and click "Execute (do)".
3. Append location-level latitude data to dataset using the merge command. Convert relative latitude to absolute latitude using the abs(x) function. Select code and click "Execute (do)".
4. Using the regress command, test the linearity of the relationship between lag and absolute latitude by running both linear and quadratic regressions and comparing residuals. Select code and click "Execute (do)".
5. Using regress, specify a linear mean regression model with UVR-trough-to-relapse-peak lag as the dependent outcome variable and absolute latitude in units of 10 degrees as the predictor variable. Weight the model for the number of patients contributed by each location using the aweights regress option. Select code and click "Execute (do)".
6. Use the twoway scatter command to plot absolute latitude on the y-axis against UVR-trough-to-relapse-lag in months on the x-axis. Overlay a line of best fit using lfit graph option. Visualise the relative patient weights of each location using the aweight analytical weights option. Select code and click "Execute (do)".

6. Sensitivity Analyses of Patient-level Relapse Propensities

1. Use the mepossion command to specify a mixed-effects Poisson regression where monthly relapse count is the dependent outcome variable, the sine and cosine trigonometric functions specified in step 2.1 are again incorporated into the model as the fixed variables, baseline EDSS, age at MS onset and prior exposure to MS specific disease-modifying treatment are included as potential confounders and unique patient identifier is specified as a random effect. Select code and click "Execute (do)".
2. Repeat steps 2.4 through 2.6 to identify the best-fitting Poisson model. This involves rerunning the twoway scatter command to overlay predicted estimates on observed data and using the regress command to run the expanded harmonic model alternatives.
3. Use the non-linear combination of estimators function (nlcom) to calculate the point estimate and 95% confidence interval for the phase-shift and calculate the date of peak relapse frequency. Compare results with the primary analysis.
4. Use the generate command to recalculate UVR-trough-to-relapse-peak lag in months for each location as described in step 5.1, using the patient-level Poisson model estimates of peak relapse date derived in step 6.3. Select code and click "Execute (do)".
5. Use the regress command to remodel absolute latitude as a predictor of lag as described in step 5.5 and compare results with the primary analysis. Select code and click "Execute (do)".

Results

The application of trigonometric regression to 32,762 relapse events sourced from 46 clinical centers across 20 countries was the basis for providing a defensible statistical argument for the observation that the timing of relapse onset in MS is cyclic and seasonal across both hemispheres and that the duration between seasonal UVR trough and subsequent relapse peak correlates with latitude. Critical to this was the reliance upon plot analysis to guide the necessarily iterative process of model development, evaluation and...

Discussion

The protocol described herein details a systematic regression based technique, guided by visual plot analysis, of global MS relapse onset data. It takes as a starting point a relatively simple descriptive analysis of relapse data from 20 countries across both hemispheres, allowing the user to explore theories regarding the temporality of relapse onset timing in MS and testing these theories formally through the use of trigonometric models. Through a step-wise process of first plotting global relapse onset data and then s...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Stata SE Version 13	StataCorp, College Station, Texas	Version 13	Statistical analysis software used for analysis
Microsoft Excel 2010	Microsoft	2010	Spreadsheet program for calendar date look-up