Oblicz prawdopodobieństwo logarytmiczne „ręcznie” dla uogólnionej regresji nieliniowej metodą najmniejszych kwadratów (NLM)

Próbuję obliczyć prawdopodobieństwo logarytmiczne dla uogólnionej regresji nieliniowej metodą najmniejszych kwadratów dla funkcji zoptymalizowanej przez funkcja w pakiecie R , przy użyciu macierzy kowariancji wariancji generowanej przez odległości na drzewie filogenetycznym przy założeniu ruchu Browna ( z pakietu). Poniższy odtwarzalny kod R pasuje do modelu GNSS przy użyciu danych x, y i losowego drzewa z 9 taksonami:$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) 

Chciałbym obliczyć logarytm prawdopodobieństwa „ręcznie” (w R, ale bez użycia logLikfunkcji) na podstawie oszacowanych parametrów uzyskanych z, gnlswięc pasuje do wyniku logLik(fit). UWAGA: Nie próbuję oszacować parametrów; Chcę po prostu obliczyć logarytmiczne prawdopodobieństwo parametrów oszacowanych przez gnlsfunkcję (chociaż jeśli ktoś ma powtarzalny przykład sposobu oszacowania parametrów bez gnls, byłbym bardzo zainteresowany jego zobaczeniem!).

Nie jestem do końca pewien, jak to zrobić w R. Liniowa notacja algebry opisana w Modelach efektów mieszanych w S i S-Plus (Pinheiro i Bates) jest bardzo ważna i żadna z moich prób się nie zgadza logLik(fit). Oto szczegóły opisane przez Pinheiro i Batesa:

prawdopodobieństwa dla uogólnionego nieliniowego modelu najmniejszych kwadratów gdzie oblicza się w następujący sposób:ϕ i = A i$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\}}$

gdzie jest liczbą obserwacji, a .$N$$f_i^*(\beta)=f_i^*(\phi_i,v_i)$

$\Lambda_i$ ma wartość dodatnią, 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)$

Dla stałych i estymator ML to$\beta$$\lambda$$\sigma^2$

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

a profilowane prawdopodobieństwo dziennika wynosi

$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\}$

który jest używany z algorytmem Gaussa-Seidela do znajdowania oszacowań ML dla i . Zastosowano mniej stronnicze oszacowanie :$\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)$

gdzie oznacza długość .$p$$\beta$

Przygotowałem listę konkretnych pytań, przed którymi stoję:

1. Co to jest ? Jest to matryca odległość wyprodukowany przez w , czy też trzeba jakoś przekształcone lub parametryzowane przez , czy coś zupełnie innego?$\Lambda_i$big_lambda <- vcv.phylo(tree)ape$\lambda$
2. Byłoby be lub równanie całkiem miarodajne oszacowania (ostatniego równania tej wiadomości)?$\sigma^2$fit$sigma^2
3. Czy konieczne jest użycie do obliczenia prawdopodobieństwa logarytmicznego, czy to tylko pośredni krok do oszacowania parametru? W jaki sposób używany jest ? Czy jest to pojedyncza wartość czy wektor i czy jest ona pomnożona przez wszystkie czy tylko elementy o przekątnej itp.?$\lambda$$\lambda$$\Lambda_i$
4. Co to jest? Czy to będzie w pakiecie ? Jeśli tak, nie jestem pewien, jak obliczyć sumę , ponieważ zwraca pojedynczą wartość, a nie wektor.M ∑ i = 1 | | y ∗ i - f ∗ i ( β ) | | $||y-f(\beta)||$norm(y-f(fit$coefficients,x),"F")Matrix$\sum\limits_{i=1}^M||y_i^*-f_i^*(\beta)||^2$norm()
5. Jak obliczyć? Czy to gdzie jest , czy pochodzi z paczki ? Jeśli tak, to jak wziąć sumę macierzy (czy sugeruje się, że są to tylko elementy ukośne)?Λ $\log|\Lambda_i|$log(diag(abs(big_lambda)))big_lambda$\Lambda_i$logm(abs(big_lambda))expmlogm()
6. Wystarczy, aby potwierdzić, czy obliczane tak: ?$\Lambda_i^{-T/2}$t(solve(sqrtm(big_lambda)))
7. Jak obliczane są i ? Czy to jedno z poniższych: f ∗ i ( β $y_i^*$$f_i^*(\beta)$

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

i

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

czy by to było

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

i

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

Jeśli na wszystkie te pytania udzielono odpowiedzi, teoretycznie myślę, że log-prawdopodobieństwo powinno być obliczalne, aby dopasować wynik logLik(fit). Każda pomoc na którekolwiek z tych pytań byłaby bardzo mile widziana. Jeśli coś wymaga wyjaśnienia, daj mi znać. Dzięki!

AKTUALIZACJA : Eksperymentowałem z różnymi możliwościami obliczania prawdopodobieństwa logarytmicznego i oto najlepsze, jakie do tej pory wymyśliłem. logLik_calcjest konsekwentnie o 1 do 3 mniejszy od wartości zwracanej przez logLik(fit). Albo jestem blisko rzeczywistego rozwiązania, albo to przez przypadek. jakieś pomysły?

  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

Zacznijmy od prostszego przypadku, w którym nie ma struktury korelacji dla reszt:

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

Prawdopodobieństwo dziennika można następnie łatwo obliczyć ręcznie za pomocą:

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

Ponieważ reszty są niezależne, możemy po prostu użyć, dnorm(..., log=TRUE)aby uzyskać poszczególne warunki prawdopodobieństwa dziennika (a następnie je podsumować). Alternatywnie możemy użyć:

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

Zauważ, że fit$sigmanie jest to „mniej stronnicze oszacowanie ” - dlatego najpierw musimy ręcznie dokonać korekty.$\sigma^2$

Teraz bardziej skomplikowany przypadek, w którym reszty są skorelowane:

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

W tym przypadku musimy użyć wielowymiarowego rozkładu normalnego. Jestem pewien, że jest gdzieś funkcja, ale zróbmy to ręcznie:

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)