Why does SAS nlmixed and R nlme give different model fit results?

Question

library(datasets)
library(nlme)
n1 <- nlme(circumference ~ phi1 / (1 + exp(-(age - phi2)/phi3)),
           data = Orange,
           fixed = list(phi1 ~ 1,
                        phi2 ~ 1,
                        phi3 ~ 1),
           random = list(Tree = pdDiag(phi1 ~ 1)),
           start = list(fixed = c(phi1 = 192.6873, phi2 = 728.7547, phi3 = 353.5323)))

I fit a nonlinear mixed effects model using nlme in R, and here's my output.

> summary(n1)
Nonlinear mixed-effects model fit by maximum likelihood
  Model: circumference ~ phi1/(1 + exp(-(age - phi2)/phi3)) 
 Data: Orange 
       AIC      BIC    logLik
  273.1691 280.9459 -131.5846

Random effects:
 Formula: phi1 ~ 1 | Tree
            phi1 Residual
StdDev: 31.48255 7.846255

Fixed effects: list(phi1 ~ 1, phi2 ~ 1, phi3 ~ 1) 
        Value Std.Error DF  t-value p-value
phi1 191.0499  16.15411 28 11.82671       0
phi2 722.5590  35.15195 28 20.55530       0
phi3 344.1681  27.14801 28 12.67747       0
 Correlation: 
     phi1  phi2 
phi2 0.375      
phi3 0.354 0.755

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.9146426 -0.5352753  0.1436291  0.7308603  1.6614518 

Number of Observations: 35
Number of Groups: 5 

I fit the same model in SAS and get the following results.

    data Orange;
    input row Tree  age circumference;
    datalines;
1     1  118            30
2     1  484            58
3     1  664            87
4     1 1004           115
5     1 1231           120
6     1 1372           142
7     1 1582           145
8     2  118            33
9     2  484            69
10    2  664           111
11    2 1004           156
12    2 1231           172
13    2 1372           203
14    2 1582           203
15    3  118            30
16    3  484            51
17    3  664            75
18    3 1004           108
19    3 1231           115
20    3 1372           139
21    3 1582           140
22    4  118            32
23    4  484            62
24    4  664           112
25    4 1004           167
26    4 1231           179
27    4 1372           209
28    4 1582           214
29    5  118            30
30    5  484            49
31    5  664            81
32    5 1004           125
33    5 1231           142
34    5 1372           174
35    5 1582           177
;


proc nlmixed data=Orange;
parms phi1=192.6873 phi2=728.7547 phi3=353.5323 vb1=991.151, s2e=61.56372;
mod = (phi1 + u1)/(1 + exp(-(age - phi2)/phi3));
model circumference ~ normal(mod, s2e);
random u1 ~ normal([0],[vb1]) subject=Tree;
run;

Can someone help me understand why I'm getting slightly different estimates? I know that the nlme uses the Lindstrom & Bates (1990) implementation. According to SAS documentation, SAS's integral approximation is based on Pinhiero & Bates (1995). I've tried changing the optimization method to Nelder-Mead to match that of nlme, but the results are still dissimilar.

I've had other cases where the standard error and parameter estimate in R vs. SAS are vastly different (I don't have a reproducible example of this, but any insight would be appreciated). I'm guessing this has to do with how nlme and nlmixed estimate the standard errors in the presence of random effects?

Answer 1

FWIW, I could reproduce the sas output using a manual optimization

########## data ################

circ <- Orange$circumference
age <- Orange$age
group <- as.numeric(Orange$Tree)
#phi1 = n1[4]$coefficients$random$Tree + 192
phi1 = 192
phi2 = 728
phi3 = 353

######### likelihood function

Likelihood <- function(x,p_age,p_circ) {
  phi1 <- x[1]
  phi2 <- x[2]
  phi3 <- x[3]
  
  fitted <- phi1/(1 + exp(-(p_age - phi2)/phi3))
  fact <- 1/(1 + exp(-(p_age - phi2)/phi3))
  resid <- p_circ-fitted
  
  sigma1 <- x[4]  #  phi1 term
  sigma2 <- x[5]  #  error term
  
  covm <- matrix(rep(0,35*35),35)  # co-variance matrix for residual terms 

  #the residuals of the group variables will be correlated in 5 7x7 blocks      
  for (k in 0:4) {
    for (l in 1:7) {
      for (m in 1:7) {
        i = l+7*k
        j = m+7*k
        if (i==j) {
          covm[i,j] <- fact[i]*fact[j]*sigma1^2+sigma2^2
        }
        else {
          covm[i,j] <- fact[i]*fact[j]*sigma1^2
        }
      }
    }
  }
  
  logd <- (-0.5 * t(resid) %*% solve(covm) %*% resid) - log(sqrt((2*pi)^35*det(covm)))
  logd
}


##### optimize

out <- nlm(function(p) -Likelihood(p,age,circ),
           c(phi1,phi2,phi3,20,8),
           print.level=1,
           iterlim=100,gradtol=10^-26,steptol=10^-20,ndigit=30) 

output

iteration = 0
Step:
[1] 0 0 0 0 0
Parameter:
[1] 192.0 728.0 353.0  30.0   5.5
Function Value
[1] 136.5306
Gradient:
[1] -0.003006727 -0.019069001  0.034154033 -0.021112696
[5] -5.669238697

iteration = 52
Parameter:
[1] 192.053145 727.906304 348.073030  31.646302   7.843012
Function Value
[1] 131.5719
Gradient:
[1] 0.000000e+00 5.240643e-09 0.000000e+00 0.000000e+00
[5] 0.000000e+00

Successive iterates within tolerance.
Current iterate is probably solution.

• So the nlmixed output is close to this optimum and it is not a different convergence thing.

• The nlme output is also close to the (different) optimum. (You can check this by changing the optimization parameters in the function call)

• I do not know exactly how nlme calculates the likelihood. It seems to be the same as in the manual optimization above. When we fill in the values from nlme output, then we get the same log-likelihood value -131.5846.

  > Likelihood(c(191.0499,722.5590,344.1621,31.48255,7.846255),age,circ)
           [,1]
  [1,] -131.5846
  

  The manual refers to an article by Lindstrom, M.J. and Bates, D.M. in Biometrics volume 46 on pp 673-687. In that article a method is described that approximates the marginal distribution of the observations (circumference in your case). The marginal distribution is integrating over all possible values of the unobserved values of the random effects. This integration is necessary because unlike with linear models, with non-linear models the random effects can not be simply linearly added. The result is that the final expression for the residuals can be different. In the code above we took the circumference minus the fitted value, but in the article they also subtract a term that arises from the linear approximation and relates to the random effect. The method is also estimating those random effects as part of the algorithm.

Answer 2

I had dealt with the same issue and agree with Martjin that you need to tweak the convergence criteria in R to make it match SAS. More specifically, you can try this combination of argument specification (in lCtr object) that I found to work pretty well in my case.

lCtr <- lmeControl(maxIter = 200, msMaxIter=200, opt='nlminb', tolerance = 1e-6, optimMethod = "L-BFGS-B")

n1 <- nlme(circumference ~ phi1 / (1 + exp(-(age - phi2)/phi3)),
           data = Orange,
           fixed = list(phi1 ~ 1,
                        phi2 ~ 1,
                        phi3 ~ 1),
           random = list(Tree = pdDiag(phi1 ~ 1)),
           start = list(fixed = c(phi1 = 192.6873, phi2 = 728.7547, phi3 = 353.5323)),
           control = lCtr)

Fair warning: this should get you the same fixed estimates between SAS and R. However, you probably wouldn't obtain the same SE of the fixed effects (which I'm still researching answers for..).