I'm fitting a linear model where the response is a function both of time and of static covariates (i.e. ones that are independent of time). The ultimate goal is to identify significant effects of the static covariates.

Is this the best general strategy for selecting variables (in R, using the nlme package)? Anything I can do better?

  1. Break the data up by groups and plot it against time. For continuous covariates, bin it and plot the data in each bin against time. Use the group-specific trends to make an initial guess at what time terms to include-- time, time^n, sin(2*pi*time)+cos(2*pi*time), log(time), exp(time), etc.
  2. Add one term at a time, comparing each model to its predecessor, never adding a higher order in the absence of lower order terms. Sin and cos are never added separately. Is it acceptable to pass over a term that significantly improves the fit of the model if there is no physical interpretation of that term?.
  3. With the full dataset, use forward selection to add static variables to the model and then relevant interaction terms with each other and with the time terms. I've seen some strong criticism of stepwise regression, but doesn't forward selection ignore significant higher order terms if the lower order terms they depend on are not significant? And I've noticed that it's hard to pick a starting model for backward elimination that isn't saturated, or singular, or fails to converge. How do you decide between variable selection algorithms?
  4. Add random effects to the model. Is this as simple as doing the variable selection using lm() and then putting the final formula into lme() and specifying the random effects? Or should I include random effects from the very start?. Compare the fits of models using a random intercept only, a random interaction with the linear time term, and random interaction with each successive time term.
  5. Plot a semivariogram to see if an autoregressive error term is needed. What should a semivariogram look like if the answer is 'no'? A horizontal line? How straight, how horizontal? Does including autoregression in the model again require checking potential variables and interactions to make sure they're still relevant?
  6. Plot the residuals to see if the variance changes as a function of fitted value, time, or any of the other terms. If it does, weigh the variances appropriately (for lme(), use the weights argument to specify a varFunc()) and compare to the unweighted model to see if this improves the fit. Is this the right sequence in which to do this step, or should it be done before autocorrelation?.
  7. Do summary() of the fitted model to identify significant coefficients for numeric covariates. Do Anova() of the fitted model to identify significant effects for qualitative covariates.

2 Answers 2


Fully data-driven model selection will result in standard errors and P-values that are too small, confidence intervals that are too narrow, and overstated effects of remaining terms in the model.

For time effects I usually model using restricted cubic splines. A detailed case study in the context of generalized least squares for correlated serial data may be found at http://biostat.mc.vanderbilt.edu/RmS - see the two attachments at the bottom named course2.pdf and rms.pdf. This uses the R rms package. The case study contains information about the choice of basis functions for the time component.

  • 1
    $\begingroup$ @Frank H: In the presence of static covariates you are suggesting that we condition their coefficients on a deterministic impact of "memory" or "time".There are two diametrically different ways of conditioning or adjusting for "memory" or "time". The second approach would be to use a suitable lag structure for the Y variate as causative/support variables. I have found that trying both ways i.e. conditioning upon dummy time trends/level shifts vs lagged memory,one can determine an optimal strategy to evaluate the significance of the static covariates. $\endgroup$
    – IrishStat
    Commented Jun 3, 2011 at 17:55
  • $\begingroup$ I may be not understanding where you're coming from, but a serial data model handles the correlation structure through a specific correlation function (usually) of the absolute difference between two times. For example it is common to use the continuous time AR1 correlation structure. You don't actually need to form lags to carry out the analysis using generalized least squares or mixed effects models. I specify the mean time profile through fixed effects of time (e.g., quadratic, cubic spline, etc.). $\endgroup$ Commented Jun 4, 2011 at 4:19
  • $\begingroup$ @Frank H: I apologize if I wasn't clear but if your forge a model as y(t)= constant+phi1*y(t-1)+phi2*y(t-2)+....phik*y(t-k) PLUS the covariates where the number and the nature of the "phis" have been analytically determined to be suffcient you are then controlling for "time dependency" . THis model might also actually include "level shifts and or time trends and/or pulses/seasonal pulses as necessary. By comparing the resultant adequacy of this approach vs your assumed approach one might be able to conclude about the best way "to condition for time" so that the covariates can be evaluated. $\endgroup$
    – IrishStat
    Commented Jun 4, 2011 at 18:42
  • 1
    $\begingroup$ That is one time series approach to the problem. For many (but not all) problems it is simpler and just as well-fitting to use a "forward look" approach where between-time correlations are described in the covariance/correlation matrix for a vector of y values across time. What I've described is the basis for mixed effects models, generalized least squares, and GEE. $\endgroup$ Commented Jun 5, 2011 at 14:03
  • 1
    $\begingroup$ Yes I think you are right if there are time-varying covariates. But the statement of the problem explicitly used the phrase "static covariates". In that setting there are two ways to model the dynamic response variable y - the way you've suggested and by specifying the mean profile function and correlation functions. $\endgroup$ Commented Jun 12, 2011 at 13:31

Fitting models that include time, time-squared, time-cubed , sines , cosines et all are not very useful in my opinion as they assume deterministic structure that often is inappropriate. Using historical values that are lagged values of the output and possibly covariate series is the approach to take. When constructing these models one needs to verify a few things ( the Gaussian Assumptions ) before declaring victory : 1. that the mean of the errors is "near zero" for all time intervals otherwise empirically identifiable structure like Pulses, Level Shifts, local time trends and/or seasonal pulses might be needed (N.B. This is not guaranteed by simply including a constant in the model ) ; 2. that there is no provable auto-correlative structure in the residuals for ALL LAGS ; 3. that the parameters of your model are stable/constant over time ; 4. that the variance of the errors is constant over time thusly no structural breaks in variance requiring weighted estimation and no dependence of the error variablility on the level of the series and no stochastic variance adaptation in effect. The suggested approach is called a Transfer Function and is a super-set of Regression and ARIMA modelling in conjunction with Intervention Detection


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.