UPDATED:Should I use lme() or gls() for a repeated measure desing?

Question

I have a problem which I have been overthinking for a long time.

Here my data structure:

• Y= quantitative variable
• T=time, 5 different cathegorical times (equally spaced)
• Group= I divide individuals in two groups based on an specific characteristic of my interest. (It's like a treatment)
• I=individual, 47 in total, but I have a lot of measures by individual (random). At first I try with a mixed effect model but because the non-normality and heterogeneity of residuals I decided try with gls to model the structure of residuals. I made the following model with nlme package:
  m1_gls<-gls(Y ~ Time*group, correlation = corAR1(form = ~ 1 | individual), bd, weights = varIdent(form = ~ 1 | Time))

I chose this model based my data and the observation of its distribution in the plot. The summary of this model is the following:

Generalized least squares fit by REML

Model: Y ~ TIME * group
Data: bd
AIC BIC logLik
1073695 1073877 -536828.3

Correlation Structure: AR(1)
Formula: ~1 | individual
Parameter estimate(s):
Phi
0.2018184
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | TIME
Parameter estimates:
0 1 2 3 4 5
1.0000000 0.9351717 0.7442982 0.6263877 0.5363402 0.4993560

Coefficients:
Value Std.Error t-value p-value (Intercept) 84.96208 0.5257654 161.59693 0.0000
TIME1 1.53680 0.7319807 2.09951 0.0358
TIME2 -18.47460 0.7040417 -26.24077 0.0000
TIME3 -33.28677 0.6986190 -47.64652 0.0000
TIME4 -40.68250 0.6884976 -59.08881 0.0000
TIME5 -47.62906 0.7114190 -66.94937 0.0000
groupL 2.38166 0.6887940 3.45773 0.0005
TIME1:groupL -6.40065 0.9600161 -6.66724 0.0000
TIME2:groupL -10.19750 0.9331267 -10.92831 0.0000
TIME3:groupL -4.98490 0.9298458 -5.36099 0.0000
TIME4:groupL -8.12306 0.9298256 -8.73611 0.0000
TIME5:groupL -3.27478 0.9769064 -3.35220 0.0008

Correlation:
(Intr) TIME1 TIME2 TIME3 TIME4 TIME5 groupL T1:L T2:L T3:L TIME1 -0.718
TIME2 -0.747 0.536
TIME3 -0.753 0.540 0.563
TIME4 -0.764 0.548 0.570 0.576
TIME5 -0.738 0.530 0.551 0.556 0.565
groupL -0.763 0.548 0.570 0.574 0.583 0.564
TIME1:groupL 0.547 -0.762 -0.409 -0.412 -0.418 -0.404 -0.717
TIME2:groupL 0.563 -0.405 -0.754 -0.425 -0.430 -0.416 -0.738 0.530
TIME3:groupL 0.565 -0.406 -0.423 -0.751 -0.432 -0.418 -0.741 0.531 0.548
TIME4:groupL 0.565 -0.406 -0.422 -0.426 -0.740 -0.418 -0.741 0.531 0.547 0.550 TIME5:groupL 0.538 -0.386 -0.402 -0.405 -0.412 -0.728 -0.704 0.505 0.520 0.522 T4:L TIME1
TIME2
TIME3
TIME4
TIME5
groupL
TIME1:groupL
TIME2:groupL
TIME3:groupL
TIME4:groupL
TIME5:groupL 0.523

Standardized residuals:
Min Q1 Med Q3 Max -1.7182284 -0.8158252 -0.0727800 0.6805442 5.6095940

Residual standard error: 47.59601
Degrees of freedom: 106806 total; 106794 residual

The plot of residuals is:

I have validated this model with cross validation. The R2 wasn't the better 0.17, while RMSE was 39.5. One of the objetive of my model is get only one value by individual that represent him. My question is, how can I extract an unique value or get the solution by individual because gls() make a covariance matrix by individual and don't allow one include random effect. For that reason I have tried with lme() including individuals like random effects and extract each coefficient (ranef()) by individual. something like this: m1_lme<-lme(Y ~ Time*group, random=~1|individual, correlation = corAR1(form = ~ 1 | individual), bd, weights = varIdent(form = ~ 1 | Time)) However, I found low variability explained by individuals (almost 5%) so I'm not sure if I should include the individuals like random effect in my model. Also, the residuals have a pattern.

Answer 1

If the time variable T is categorical discrete, you could define it as an ordered() factor in R. Then automatically, contr.poly() will be used to expand the ordered time variables using orthogonal polynomials. From the output, you could then see how many polynomial terms you need to keep.

If the time variable is continuous, you can still use orthogonal polynomials for it using, e.g., poly(T, 3) for cubic polynomials. Alternatively, you could use splines via the bs() and ns() function provided by the splines package.