R
為什麼 SAS nlmixed 和 R nlme 給出不同的模型擬合結果？
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)))我在 R 中擬合了一個非線性混合效應模型
nlme，這是我的輸出。> 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我在 SAS 中擬合相同的模型並得到以下結果。
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;
有人可以幫我理解為什麼我得到的估計略有不同嗎？我知道
nlme使用 Lindstrom & Bates (1990) 實現。根據 SAS 文檔，SAS 的積分近似值基於 Pinhiero & Bates (1995)。我嘗試將優化方法更改為 Nelder-Mead 以匹配nlme，但結果仍然不同。我遇到過其他情況，其中 R 與 SAS 中的標準誤差和參數估計有很大不同（我沒有可重複的例子，但任何見解都會受到讚賞）。我猜這與在存在隨機效應的情況下如何
nlme估計nlmixed標準誤差有關？
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(-(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 是如何計算可能性的（儘管該值幾乎相同 -131.6），但我懷疑它與上述擬合 3 個參數（固定效應）和 2 個有害參數不同。使用對隨機效應使用附加參數的似然函數，我可以得到與它相似但不完全相似的結果。我想我對討厭的參數的處理方式不同（很可能我犯了一個錯誤）。