Unstructured model vs random slope model for repeated measures based on R functions lmer, lme and gls

Here is the dataset for repeated measures:

library(lme4)
library(nlme)

d$$Program = factor(d$$Program)
d$$Subj = factor(d$$Subj)
d$$Timef = factor(d$$Time)


I have built an unstructured model using gls function and I am trying to reproduce the results using a random slope (or random slope + intercept) model fitted by R functions lmer or lme. Here is the code and output:

gls function:

> d.gls <- gls(Strength ~ Program * Timef, data = d,
+     correlation = corSymm(form =~ 1 | Subj),
+     weight = varIdent(form = ~ 1 | Timef))
> getVarCov(d.gls)
Marginal variance covariance matrix
[,1]   [,2]    [,3]    [,4]    [,5]    [,6]    [,7]
[1,] 8.7801 8.7571  8.9656  8.1984  8.6781  8.2203  8.4169
[2,] 8.7571 9.4730  9.4631  8.5686  9.2012  8.7307  8.6875
[3,] 8.9656 9.4631 10.7080  9.9266 10.6660 10.0700 10.2140
[4,] 8.1984 8.5686  9.9266 10.0770 10.6000  9.8987 10.0430
[5,] 8.6781 9.2012 10.6660 10.6000 12.0950 11.3440 11.3640
[6,] 8.2203 8.7307 10.0700  9.8987 11.3440 11.7560 11.6500
[7,] 8.4169 8.6875 10.2140 10.0430 11.3640 11.6500 12.7100
Standard Deviations: 2.9631 3.0778 3.2723 3.1745 3.4778 3.4287 3.5651
> logLik(d.gls)
'log Lik.' -617.4479 (df=49)


lmer function:

   > d.lm <- lmer(Strength ~ Program * Timef +(0+Timef|Subj),d,control = lmerControl(check.nobs.vs.nRE = "ignore"))
Warning message:
In checkConv(attr(opt, "derivs"), opt$$par, ctrl = control$$checkConv,  :
Model failed to converge with max|grad| = 0.00271822 (tol = 0.002, component 1)
> #VarCorr(d.lm)
> as.matrix(Matrix::bdiag(VarCorr(d.lm)))
Timef2   Timef4    Timef6    Timef8   Timef10   Timef12   Timef14
Timef2  8.534639 8.756711  8.965454  8.198261  8.677920  8.220413  8.416961
Timef4  8.756711 9.227384  9.462754  8.568308  9.201094  8.730927  8.687653
Timef6  8.965454 9.462754 10.462654  9.926419 10.665922 10.070205 10.213872
Timef8  8.198261 8.568308  9.926419  9.832025 10.599504  9.898893 10.043531
Timef10 8.677920 9.201094 10.665922 10.599504 11.849795 11.344650 11.364012
Timef12 8.220413 8.730927 10.070205  9.898893 11.344650 11.511176 11.650515
Timef14 8.416961 8.687653 10.213872 10.043531 11.364012 11.650515 12.465186
> logLik(d.lm)
'log Lik.' -617.4479 (df=50)


lme function is super slow and fails to converge. My questions are:

1. Considering that gls function does not produce a warning whereas lmer function does, do we trust the gls result rather than lmer result?
2. Are the two models by gls and lmer are the same unstructured model ? If not, why their log likelihood is the same (-617)? If yes, why the variance-covariance estimates are different?
3. Is it possible to reproduce the gls unstructured model using lme/lmer function?

UPDATE I was able to generate the results using lme() through a different optimizer "optim":

    > d.lme <- lme(Strength ~ Program * Timef, random = ~ -1+Timef | Subj, d,control = lmeControl(maxIter = 50, msMaxIter = 50, msVerbose = TRUE,opt='optim'))
initial  value 1202.821021
iter  10 value 1202.809062
iter  20 value 1202.789065
final  value 1202.784342
converged
> #VarCorr(d.lme)
> getVarCov(d.lme)
Random effects variance covariance matrix
Timef2 Timef4  Timef6  Timef8 Timef10 Timef12 Timef14
Timef2  8.5759 8.7549  8.9628  8.1949  8.6764  8.2180  8.4140
Timef4  8.7549 9.2691  9.4601  8.5649  9.1991  8.7284  8.6841
Timef6  8.9628 9.4601 10.5030  9.9229 10.6630 10.0680 10.2100
Timef8  8.1949 8.5649  9.9229  9.8725 10.5970  9.8965 10.0400
Timef10 8.6764 9.1991 10.6630 10.5970 11.8910 11.3420 11.3600
Timef12 8.2180 8.7284 10.0680  9.8965 11.3420 11.5510 11.6480
Timef14 8.4140 8.6841 10.2100 10.0400 11.3600 11.6480 12.5060
Standard Deviations: 2.9285 3.0445 3.2408 3.1421 3.4483 3.3987 3.5363
> logLik(d.lme)
'log Lik.' -617.4481 (df=50)

• I can't reproduce your results. There is no Timef variable in the data frame (I guess it's the factor recoding of Time). You might want to add library information too (I can see them on tags.) May 17 at 18:05
• Sorry, I missed 3 lines to transform the variables to factors. Please let me know if it still does not work. Thanks! May 17 at 19:09
• These seem to be different models. There's a critical error message if you omit the control argument to lmer(): "Error: number of observations (=399) <= number of random effects (=399) for term (0 + Timef | Subj); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable." Please edit the question to explain in more detail why you think that should be the same model as with gls; please don't respond only in a comment, as comments are easy to overlook and can be deleted.
– EdM
May 17 at 19:52
• @EdM, I do not know if these are the same model. If they are not, is it possible to reproduce the results in lmer and lme (Q3) ? May 17 at 20:19

The gls and lme models have an important difference. As this Stack Overflow answer quotes from Pinheiro and Bates:

The gls function... can be veiwed as an lme function without the argument random.

Although the matrix of error correlations among observations in your gls model has no imposed internal structure, it has a very strict structure within the model: except for your allowing different variances within Timef groups via the VarIdent() you supplied to the weights argument, the very same correlation matrix is assumed for all individuals.

The lme model that eventually worked for you assumes by default no within-individual correlations. It instead models individuals as having random Timef coefficients that have a strictly normal distribution within each value of Timef.

Although in this case the fixed-effect coefficients seem to be the same in the gls and lme models (as in the linked Stack Overflow answer), I don't think that you can necessarily expect the "unstructured" error-correlation matrix in gls to be exactly "reproduced" by the variance-covariance matrix of the sets of Timef coefficients that are normally distributed among individuals, although they do come close here. (That might be the case asymptotically in the limit of large numbers of observations.)

Maybe more important is the difference in the ways to think about the models. This answer emphasizes the marginal (gls) versus conditional (lme) interpretations, and has a link to why some results reported for lmer models aren't the same as for corresponding lme models.