7
$\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. enter image description here

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$
  • $\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$ – Martijn Weterings 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$ – Martijn Weterings 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
3
$\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.