Hitung log-kemungkinan "dengan tangan" untuk regresi kuadrat terkecil nonlinier umum (nlme)

Saya mencoba menghitung log-kemungkinan untuk regresi kuadrat terkecil nonlinear umum untuk fungsi dioptimalkan olehfungsi dalam paket R, menggunakan matriks varians kovarians yang dihasilkan oleh jarak pada pohon filogenetik dengan asumsi gerakan Brown (daripaket). Kode R yang dapat direproduksi berikut cocok dengan model gnls menggunakan x, data y dan pohon acak dengan 9 taksa:$f(x)=\frac{\beta_1}{(1+\frac x\beta_2)^{\beta_3}}$gnlsnlmecorBrownian(phy=tree)ape

require(ape)
require(nlme)
require(expm)
tree <- rtree(9)
x <- c(0,14.51,32.9,44.41,86.18,136.28,178.21,262.3,521.94)
y <- c(100,93.69,82.09,62.24,32.71,48.4,35.98,15.73,9.71)
data <- data.frame(x,y,row.names=tree$tip.label)
model <- y~beta1/((1+(x/beta2))^beta3)
f=function(beta,x) beta[1]/((1+(x/beta[2]))^beta[3])
start <- c(beta1=103.651004,beta2=119.55067,beta3=1.370105)
correlation <- corBrownian(phy=tree)
fit <- gnls(model=model,data=data,start=start,correlation=correlation)
logLik(fit) 

Saya ingin menghitung log-likelihood "dengan tangan" (dalam R, tetapi tanpa menggunakan logLikfungsi) berdasarkan pada estimasi parameter yang diperoleh gnlssehingga cocok dengan output dari logLik(fit). CATATAN: Saya tidak mencoba memperkirakan parameter; Saya hanya ingin menghitung log-kemungkinan parameter yang diperkirakan oleh gnlsfungsi (walaupun jika seseorang memiliki contoh yang dapat direproduksi tentang cara memperkirakan parameter tanpa gnls, saya akan sangat tertarik melihatnya!).

Saya tidak begitu yakin bagaimana cara melakukan ini di R. Notasi aljabar linier yang dijelaskan dalam Model Efek Campuran di S dan S-Plus (Pinheiro dan Bates) sangat banyak di kepala saya dan tidak ada upaya saya yang cocok logLik(fit). Berikut detail yang dijelaskan oleh Pinheiro dan Bates:

Log-kemungkinan untuk model kuadrat terkecil nonlinier umum mana ϕ i = A i β dihitung sebagai berikut:$y_i=f_i(\phi_i,v_i)+\epsilon_i$$\phi_i=A_i\beta$

$l(\beta,\sigma^2,\delta|y)=-\frac 12 \Bigl\{ N\log(2\pi\sigma^2)+\sum\limits_{i=1}^M{\Bigl[\frac{||y_i^*-f_i^*(\beta)||^2}{\sigma^2}+\log|\Lambda_i|\Bigl]\Bigl\}}$

di mana adalah jumlah pengamatan, dan f ∗ i ( β ) = f ∗ i ( ϕ i , v i ) .$N$$f_i^*(\beta)=f_i^*(\phi_i,v_i)$

pasti-positif, y ∗ i = Λ - T / 2 i y i dan f ∗ i ( ϕ i , v i ) = Λ - T / 2 i f i ( ϕ i , v i$\Lambda_i$$y_i^*=\Lambda_i^{-T/2}y_i$$f_i^*(\phi_i,v_i)=\Lambda_i^{-T/2}f_i(\phi_i,v_i)$

Untuk tetap dan λ , estimator ML σ 2 adalah$\beta$$\lambda$$\sigma^2$

$\hat\sigma(\beta,\lambda)=\sum\limits_{i=1}^M||y_i^*-f_i^*(\beta)||^2 / N$

dan kemungkinan log yang diprofilkan adalah

$l(\beta,\lambda|y)=-\frac12\Bigl\{N[\log(2\pi/N)+1]+\log\Bigl(\sum\limits_{i=1}^M||y_i^*-f_i^*(\beta)||^2\Bigl)+\sum\limits_{i=1}^M\log|\Lambda_i|\Bigl\}$

$\beta$$\lambda$$\sigma^2$

$\sigma^2=\sum\limits_{i=1}^M\Bigl|\Bigl|\hat\Lambda_i^{-T/2}[y_i-f_i(\hat\beta)]\Bigl|\Bigl|^2/(N-p)$

$p$$\beta$

Saya telah menyusun daftar pertanyaan spesifik yang saya hadapi:

1. $\Lambda_i$big_lambda <- vcv.phylo(tree)ape$\lambda$
2. $\sigma^2$fit$sigma^2
3. $\lambda$$\lambda$$\Lambda_i$
4. $||y-f(\beta)||$norm(y-f(fit$coefficients,x),"F")Matrix$\sum\limits_{i=1}^M||y_i^*-f_i^*(\beta)||^2$norm()
5. $\log|\Lambda_i|$log(diag(abs(big_lambda)))big_lambda$\Lambda_i$logm(abs(big_lambda))expmlogm()
6. $\Lambda_i^{-T/2}$t(solve(sqrtm(big_lambda)))
7. $y_i^*$$f_i^*(\beta)$

y_star <- t(solve(sqrtm(big_lambda))) %*% y

dan

f_star <- t(solve(sqrtm(big_lambda))) %*% f(fit$coefficients,x)

atau apakah itu

y_star <- t(solve(sqrtm(big_lambda))) * y

dan

f_star <- t(solve(sqrtm(big_lambda))) * f(fit$coefficients,x) ?

Jika semua pertanyaan ini dijawab, secara teori, saya pikir kemungkinan log harus dapat dihitung agar sesuai dengan keluaran logLik(fit). Bantuan apa pun dari semua pertanyaan ini akan sangat dihargai. Jika ada yang butuh klarifikasi, beri tahu saya. Terima kasih!

UPDATE : Saya telah bereksperimen dengan berbagai kemungkinan untuk perhitungan log-likelihood, dan di sini adalah yang terbaik yang saya dapatkan sejauh ini. logLik_calcsecara konsisten sekitar 1 hingga 3 dari nilai yang dikembalikan oleh logLik(fit). Entah saya dekat dengan solusi yang sebenarnya, atau ini murni karena kebetulan. Adakah pikiran?

  C <- vcv.phylo(tree) # variance-covariance matrix
  tC <- t(solve(sqrtm(C))) # C^(-T/2)
  log_C <- log(diag(abs(C))) # log|C|
  N <- length(y)
  y_star <- tC%*%y 
  f_star <- tC%*%f(fit$coefficients,x)
  dif <- y_star-f_star  
  sigma_squared <-  sum(abs(y_star-f_star)^2)/N
  # using fit$sigma^2 also produces a slightly different answer than logLik(fit)
  logLik_calc <- -((N*log(2*pi*(sigma_squared)))+
       sum(((abs(dif)^2)/(sigma_squared))+log_C))/2

Mari kita mulai dengan kasus sederhana di mana tidak ada struktur korelasi untuk residu:

fit <- gnls(model=model,data=data,start=start)
logLik(fit)

Kemungkinan log kemudian dapat dengan mudah dihitung dengan tangan dengan:

N <- fit$dims$N
p <- fit$dims$p
sigma <- fit$sigma * sqrt((N-p)/N)
sum(dnorm(y, mean=fitted(fit), sd=sigma, log=TRUE))

Karena residu bersifat independen, kami hanya dapat menggunakan dnorm(..., log=TRUE)untuk mendapatkan istilah kemungkinan log individual (dan kemudian menjumlahkannya). Atau, kita bisa menggunakan:

sum(dnorm(resid(fit), mean=0, sd=sigma, log=TRUE))

fit$sigma$\sigma^2$

Sekarang untuk kasus yang lebih rumit di mana residu berkorelasi:

fit <- gnls(model=model,data=data,start=start,correlation=correlation)
logLik(fit)

Di sini, kita perlu menggunakan distribusi normal multivariat. Saya yakin ada fungsi untuk ini di suatu tempat, tapi mari kita lakukan ini dengan tangan:

N <- fit$dims$N
p <- fit$dims$p
yhat <- cbind(fitted(fit))
R <- vcv(tree, cor=TRUE)
sigma <- fit$sigma * sqrt((N-p)/N)
S <- diag(sigma, nrow=nrow(R)) %*% R %*% diag(sigma, nrow=nrow(R))
-1/2 * log(det(S)) - 1/2 * t(y - yhat) %*% solve(S) %*% (y - yhat) - N/2 * log(2*pi)