In plain English: I have a multiple regression or ANOVA model but the response variable for each individual is a curvilinear function of time.

  • How can I tell which of the right-hand-side variables are responsible for significant differences in the shapes or vertical offsets of the curves?
  • Is this a time-series problem, a repeated-measures problem, or something else entirely?
  • What are the best-practices for analyzing such data (preferably in R, but I'm open to using other software)?

In more precise terms: Let's say I have a model $y_{ijk} = \beta_0 + \beta_1 x_i + \beta_2 x_j + \beta_3 x_i x_j + \epsilon_k$ but $y_{ijk}$ is actually a series of data-points collected from the same individual $k$ at many time-points $t$, which were recorded as a numeric variable. Plotting the data shows that for each individual $y_{ijkt}$ is a quadratic or cyclical function of time whose vertical offset, shape, or frequency (in the cyclical case) might significantly depend on the covariates. The covariates do not change over time-- i.e., an individual has a constant body weight or treatment group for the duration of the data collection period.

So far I have tried the following R approaches:

  1. Manova


    ...where YT is a matrix whose columns are the time points, 10 of them in this example, but far more in the real data.

    Problem: this treats time as a factor, but the time-points don't exactly match for each individual. Furthermore, there are many of them relative to the sample size so the model gets saturated. It seems like the shape of the response variable over time is ignored.

  2. Mixed-model (as in Pinheiro and Bates, Mixed Effect Models in S and S-Plus)

    lme(fixed=Y~ A*B*TIME + sin(2*pi*TIME) + cos(2*pi*TIME), data=mydata, 
        random=~(TIME + sin(2*pi*TIME) + cos(2*pi*TIME))|ID), method='ML')

    ...where ID is a factor that groups data by individual. In this example the response is cyclical over time, but there could instead be quadratic terms or other functions of time.

    Problem: I'm not certain whether each time term is necessary (especially for quadratic terms) and which ones are affected by which covariates.

    • Is stepAIC() a good method for selecting them?
    • If it does remove a time-dependent term, will it also remove it from the random argument?
    • What if I also use an autocorrelation function (such as corEXP()) that takes a formula in the correlation argument-- should I make that formula for corEXP() the same as the one in random or just ~1|ID?
    • The nlme package is rarely mentioned in the context of time series outside Pinheiro and Bates... is it not considered well suited to this problem?
  3. Fitting a quadratic or trigonometric model to each individual, and then using each coefficient as a response variable for multiple regression or ANOVA.

    Problem: Multiple comparison correction necessary. Can't think of any other problems which makes me suspicious that I'm overlooking something.

  4. As previously suggested on this site (What is the term for a time series regression having more than one predictor?), there are ARIMAX and transfer function / dynamic regression models.

    Problem: ARMA-based models assume discrete times, don't they? As for dynamic regression, I heard about it for the first time today, but before I delve into yet another new method that might not pan out after all, I thought it would be prudent to ask people who have done this before for advice.

  • 5
    $\begingroup$ @f1r3br4and It may help if you add some details about exactly what is your data? I.e., how many cases measured at how many time points? Are the cases in different conditions? or something different? $\endgroup$ Commented Jun 1, 2011 at 4:29
  • $\begingroup$ Have you lost your account? If so, please register new one, possibly with your Gmail address so I could easily find what to merge. $\endgroup$
    – user88
    Commented Jun 1, 2011 at 13:08
  • 1
    $\begingroup$ @f-tussel and @jeromy-anglim: The data are various non-invasive measurements collected from lab animals over the course of their natural lifespans, often coupled with static information such as genotype, sex, or the type of diet they are fed. By 'many' points I mean from 20 to several hundred. It's not just one data set I'm trying to analyze, but to learn how to analyze a type of data that is becoming very common in my group. Are mixed models with within-group factors only recommended for smaller sample sizes? $\endgroup$
    – f1r3br4nd
    Commented Jun 1, 2011 at 13:59

2 Answers 2


As Jeromy Anglim said, it would help to know the number of time points you have for each individual; as you said "many" I would venture that functional analysis might be a viable alternative. You might want to check the R package fda and look at the book by Ramsay and Silverman.

  • 1
    $\begingroup$ Functional analysis sounds promising in the long term, but it looks like there will be quite a learning curve before I'm confident that I'm not getting meaningless or biased results. So while I'm getting up to speed on this, are the more familiar nlme-based approaches (items 2 and 3 in the OP) at least valid to use on data for the time being? $\endgroup$
    – f1r3br4nd
    Commented Jun 2, 2011 at 16:45

Since originally posing this question, I have come to the conclusion that mixed-effect models with subjects as the random blocking factor are the practical solution to this problem, i.e. option #2 in my original post. If the random argument to lme is set to ~1|ID (where ID identifies observations coming from the same test subject) then a random intercept model is fitted. If it is set to ~TIME|ID then a random slope and intercept model is fitted. Any right-sided formula containing variables that vary within the same individual can be placed between the ~ and the |ID, but overly complicated formulas will result in a saturated model and/or various numerical errors. Therefore, one can use a likelihood ratio test (anova(myModel, update(myModel,random=~TIME|ID))) to compare a random intercept model to a random slope and intercept model or other candidate random effect models. If the difference in fit is not significant then stick with the simpler model. It was overkill for me to go into random trig functions in my original post.

The other issue I raised was one of model selection. It seems like people don't like model selection of any kind, but nobody has any practical alternatives. If you blindly believe the researcher who collected the data about what explanatory variables are and are not relevant, you will often be blindly accepting their untested assumptions. If you take into account every possible bit of information, you will often end up with a saturated model. If you arbitrarily choose a particular model and variables because they're easy, you will again be accepting untested assumptions, this time your own.

So, in summary, for repeated measures it's lme models followed by trimming via MASS:::stepAIC or MuMIn:::dredge and/or nlme:::anova.lme until and unless someone has a better idea.

I'll leave this self-answer up for a while before accepting it to see if anybody has any rebuttals. Thanks for your time, and if you're reading this because you have the same sort of question I have, good luck and welcome to semi-uncharted territory.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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