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

Question

Here is the dataset for repeated measures:

library(lme4)
library(nlme)

d=read.delim("http://dnett.github.io/S510/RepeatedMeasures.txt")
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)

Answer 1

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.