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:gnls
nlme
corBrownian(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 logLik
funkcji) na podstawie oszacowanych parametrów uzyskanych z, gnls
więc pasuje do wyniku logLik(fit)
. UWAGA: Nie próbuję oszacować parametrów; Chcę po prostu obliczyć logarytmiczne prawdopodobieństwo parametrów oszacowanych przez gnls
funkcję (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 β
gdzie jest liczbą obserwacji, a .
ma wartość dodatnią, i
Dla stałych i estymator ML to
a profilowane prawdopodobieństwo dziennika wynosi
który jest używany z algorytmem Gaussa-Seidela do znajdowania oszacowań ML dla i . Zastosowano mniej stronnicze oszacowanie :
gdzie oznacza długość .
Przygotowałem listę konkretnych pytań, przed którymi stoję:
- 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?
big_lambda <- vcv.phylo(tree)
ape
- Byłoby be lub równanie całkiem miarodajne oszacowania (ostatniego równania tej wiadomości)?
fit$sigma^2
- 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.?
- 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 ( β ) | | 2)
norm(y-f(fit$coefficients,x),"F")
Matrix
norm()
- 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)?Λ i
log(diag(abs(big_lambda)))
big_lambda
logm(abs(big_lambda))
expm
logm()
- Wystarczy, aby potwierdzić, czy obliczane tak: ?
t(solve(sqrtm(big_lambda)))
- Jak obliczane są i ? Czy to jedno z poniższych: f ∗ i ( β )
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_calc
jest 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
Odpowiedzi:
Zacznijmy od prostszego przypadku, w którym nie ma struktury korelacji dla reszt:
Prawdopodobieństwo dziennika można następnie łatwo obliczyć ręcznie za pomocą:
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ć:Zauważ, żeσ2
fit$sigma
nie jest to „mniej stronnicze oszacowanie ” - dlatego najpierw musimy ręcznie dokonać korekty.Teraz bardziej skomplikowany przypadek, w którym reszty są skorelowane:
W tym przypadku musimy użyć wielowymiarowego rozkładu normalnego. Jestem pewien, że jest gdzieś funkcja, ale zróbmy to ręcznie:
źródło
vcv
funkcją) - ale musisz uzyskać macierz korelacji, a następnie użyć aby przekształcić ją w macierz var-cov.