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?
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.
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?.
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?
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.
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?
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?.
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.
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.
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