10
$\begingroup$
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;

enter image description here

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?

$\endgroup$
4
  • $\begingroup$ It is interesting to see that the sas model uses somehow 4 degrees of freedom for the estimate of the standard error / deviation. Why not 27 or 28? How many observations are there in the dataset uses for the sas model? $\endgroup$ Mar 28 '18 at 12:51
  • $\begingroup$ @MartijnWeterings That is indeed intriguing...The Orange data set contains 35 observations. $\endgroup$
    – Adrian
    Mar 28 '18 at 16:00
  • $\begingroup$ There are some quirks in determining the DF so it might be due to that. Anyway, there might be more than the DF stuff (which, I believe, does not influence the model fitting)... I have been trying to manually fit a log-likelihood function and I can not get exactly the same as either nlme or nlmixed. I believe that the differences are in the used loglikelihood function and the used method to optimize it. $\endgroup$ Mar 28 '18 at 16:20
  • $\begingroup$ In my mind they are very close. You initiated the parms. Have you compared the trace outputs? R (or SAS) may have different convergence criteria, so the laxer one calls it quits earlier while the other skips ahead a few more iterations. $\endgroup$
    – AdamO
    Mar 28 '18 at 17:15
5
$\begingroup$

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(-(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 (although the value is nearly the same -131.6) but I suspect that it is different as the above fitting 3 parameters (the fixed effects) and 2 nuisance parameters. Using a likelihood function that uses additional parameters for the random effect, I could get a result that resembles it but not exactly. I guess that I was dealing differently with the nuisance parameters (and likely I made an error).
$\endgroup$
3
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.