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

12

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)=β1(1+xβ2)β3gnlsnlmecorBrownian(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 βyi=fi(ϕi,vi)+ϵiϕi=Aiβ

l(β,σ2,δ|y)=12{Nlog(2πσ2)+i=1M[||yifi(β)||2σ2+log|Λi|]}

gdzie jest liczbą obserwacji, a .Nfi(β)=fi(ϕi,vi)

Λi ma wartość dodatnią, iyi=ΛiT/2yifi(ϕi,vi)=ΛiT/2fi(ϕi,vi)

Dla stałych i estymator ML toβλσ2

σ^(β,λ)=i=1M||yifi(β)||2/N

a profilowane prawdopodobieństwo dziennika wynosi

l(β,λ|y)=12{N[log(2π/N)+1]+log(i=1M||yifi(β)||2)+i=1Mlog|Λi|}

który jest używany z algorytmem Gaussa-Seidela do znajdowania oszacowań ML dla i . Zastosowano mniej stronnicze oszacowanie :βλσ2

σ2=i=1M||Λ^iT/2[yifi(β^)]||2/(Np)

gdzie oznacza długość .pβ

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?Λibig_lambda <- vcv.phylo(tree)apeλ
  2. Byłoby be lub równanie całkiem miarodajne oszacowania (ostatniego równania tej wiadomości)?σ2fit$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.?λλΛ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 ( β ) | | 2)||yf(β)||norm(y-f(fit$coefficients,x),"F")Matrixi=1M||yifi(β)||2norm()
  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)?Λ ilog|Λi|log(diag(abs(big_lambda)))big_lambdaΛilogm(abs(big_lambda))expmlogm()
  6. Wystarczy, aby potwierdzić, czy obliczane tak: ?ΛiT/2t(solve(sqrtm(big_lambda)))
  7. Jak obliczane są i ? Czy to jedno z poniższych: f i ( β )yifi(β)

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
Eric
źródło
w twojej definicji funkcji brakuje po prawej stronie. xf(x)x
Glen_b

Odpowiedzi:

10

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.σ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)
Wolfgang
źródło
Prawdopodobieństwo logarytmiczne dla nieskorelowanych reszt działało idealnie, jednak nie mogę ustalić wielowymiarowego rozkładu normalnego. W takim przypadku co to jest S? Próbowałem S <- vcv.phylo (drzewo) i dostałem około -700 dla prawdopodobieństwa dziennika, podczas gdy logLik (dopasowanie) wynosił około -33.
Eric
Przepraszam - popsułem kod, kiedy wkleiłem kod. Teraz jest kompletny. S jest macierzą wariancji-kowariancji reszt. Byłeś na dobrej drodze (z vcvfunkcją) - ale musisz uzyskać macierz korelacji, a następnie użyć aby przekształcić ją w macierz var-cov. σ^2
Wolfgang