FWIW，我可以使用手动优化重现 sas 输出

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

输出

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.

• 所以 nlmixed 输出接近这个最优值，它不是一个不同的收敛事物。

• nlme 输出也接近（不同的）最优值。（您可以通过更改函数调用中的优化参数来检查这一点）

• 我不知道 nlme 是如何计算可能性的。好像和上面手动优化的一样。当我们从 nlme 输出中填充值时，我们得到相同的对数似然值 -131.5846。

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

  手册_参考 Lindstrom, MJ 和 Bates, DM 在 Biometrics 第 46 卷第 673-687 页上的一篇文章。在那篇文章中描述了一种近似观察值的边际分布的方法（在您的情况下为周长）。边际分布是对未观察到的随机效应值的所有可能值进行积分。这种整合是必要的，因为与线性模型不同，非线性模型不能简单地线性添加随机效应。结果是残差的最终表达式可能不同。在上面的代码中，我们将周长减去拟合值，但在本文中，它们还减去了一个由线性近似产生的与随机效应相关的项。该方法还将这些随机效应估计为算法的一部分。

我处理过同样的问题并同意 Martjin 的观点，即您需要调整 R 中的收敛标准以使其与 SAS 匹配。更具体地说，您可以尝试这种我发现在我的情况下工作得很好的参数规范组合（在 lCtr 对象中）。

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)

公平警告：这应该让您在 SAS 和 R 之间获得相同的固定估计值。但是，您可能不会获得相同的固定效应 SE（我仍在研究答案......）。