Berechnen Sie die log-Wahrscheinlichkeit „von Hand“ für die verallgemeinerte nichtlineare Regression der kleinsten Quadrate (nlme)

12

Ich versuche, die log-Wahrscheinlichkeit für eine verallgemeinerte nichtlineare Regression der kleinsten Quadrate für die Funktion f ( x ) = β 1 zu berechnenoptimiert durch dieFunktion im R-Paketunter Verwendung der Varianz-Kovarianz-Matrix, die durch Abstände auf einem phylogenetischen Baum unter Annahme einer Brownschen Bewegung (aus demPaket) erzeugt wird. Der folgende reproduzierbare R-Code passt unter Verwendung von x-, y-Daten und einem zufälligen Baum mit 9 Taxa zum gnls-Modell: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) 

Ich möchte die log-Wahrscheinlichkeit "von Hand" berechnen (in R, aber ohne Verwendung des logLik Funktion) basierend auf den geschätzten Parametern erhalten wurden, gnlssodass sie mit der Ausgabe von übereinstimmen logLik(fit). HINWEIS: Ich versuche nicht, Parameter zu schätzen. Ich möchte nur die Log-Wahrscheinlichkeit der von der gnlsFunktion geschätzten Parameter berechnen (obwohl gnlsich sehr interessiert wäre, wenn jemand ein reproduzierbares Beispiel für die Schätzung von Parametern ohne hat!).

Ich bin mir nicht sicher, wie ich das in R machen soll. Die lineare Algebra-Notation, die in Mixed-Effects-Modellen in S und S-Plus (Pinheiro und Bates) beschrieben ist, geht mir sehr weit über den Kopf und keiner meiner Versuche hat gepasst logLik(fit). Hier sind die von Pinheiro und Bates beschriebenen Details:

Die log-Wahrscheinlichkeit für das verallgemeinerte nichtlineare Modell der kleinsten Quadrate wo φ i = A i βyi=fi(ϕi,vi)+ϵiϕi=Aiβ wird wie folgt berechnet:

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

wobei die Anzahl der Beobachtungen ist und f i ( β ) = fN.fi(β)=fi(ϕi,vich)

Λichyich=Λich-T/2yichfich(ϕich,vich)=Λich-T/2fich(ϕich,vich)

βλσ2

σ^(β,λ)=ich=1M||yich-fich(β)||2/N

und die profilierte Log-Wahrscheinlichkeit ist

l(β,λ|y)=-12{N[Log(2π/N)+1]+Log(ich=1M||yich-fich(β)||2)+ich=1MLog|Λich|}

βλσ2

σ2=ich=1M||Λ^ich-T/2[yich-fich(β^)]||2/(N-p)

pβ

Ich habe eine Liste mit spezifischen Fragen zusammengestellt, mit denen ich konfrontiert bin:

  1. Λichbig_lambda <- vcv.phylo(tree)apeλ
  2. σ2 wäre fit$sigma^2, oder die Gleichung für die weniger voreingenommen Schätzung (die letzte Gleichung in diesem Beitrag)?
  3. Ist es notwendig zu verwenden λλΛich
  4. ||y-f(β)||norm(y-f(fit$coefficients,x),"F")Matrixich=1M||yich-fich(β)||2, weil norm()ein einzelner Wert zurückgegeben wird, kein Vektor.
  5. Wie rechnet man Log|Λich|? Ist es log(diag(abs(big_lambda)))wo big_lambdaistΛich, oder ist es logm(abs(big_lambda))aus dem Paket expm? Wenn ja logm(), wie nimmt man die Summe einer Matrix (oder impliziert man, dass es sich nur um die diagonalen Elemente handelt)?
  6. Nur um zu bestätigen, ist Λich-T/2berechnet wie folgt aus : t(solve(sqrtm(big_lambda)))?
  7. Wie bist yich und fich(β)berechnet? Ist es einer der folgenden:

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

und

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

oder wäre es

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

und

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

Wenn alle diese Fragen beantwortet sind, sollte theoretisch die Log-Wahrscheinlichkeit kalkulierbar sein, um mit der Ausgabe von übereinzustimmen logLik(fit). Jede Hilfe zu diesen Fragen wäre sehr dankbar. Wenn etwas geklärt werden muss, lassen Sie es mich bitte wissen. Vielen Dank!

UPDATE : Ich habe mit verschiedenen Möglichkeiten für die Berechnung der Log-Wahrscheinlichkeit experimentiert, und hier ist das Beste, das ich mir bisher ausgedacht habe. logLik_calcist konsistent etwa 1 bis 3 vom zurückgegebenen Wert von logLik(fit). Entweder bin ich der eigentlichen Lösung nahe, oder das ist rein zufällig. Irgendwelche Gedanken?

  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
quelle
Ihre Definition der Funktion f(x) fehlt ein xauf der rechten Seite.
Glen_b -Reinstate Monica

Antworten:

10

Beginnen wir mit dem einfacheren Fall, in dem es keine Korrelationsstruktur für die Residuen gibt:

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

Die Log-Wahrscheinlichkeit kann dann einfach von Hand berechnet werden mit:

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

Da die Residuen unabhängig sind, können wir einfach dnorm(..., log=TRUE)die einzelnen Log-Likelihood-Terme abrufen (und sie dann zusammenfassen). Alternativ könnten wir verwenden:

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

Beachten Sie, dass dies fit$sigmanicht die "weniger voreingenommene Schätzung von" istσ2"- also müssen wir die Korrektur zuerst manuell vornehmen.

Nun zum komplizierteren Fall, in dem die Residuen korreliert sind:

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

Hier müssen wir die multivariate Normalverteilung verwenden. Ich bin mir sicher, dass es dafür irgendwo eine Funktion gibt, aber machen wir das einfach von Hand:

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
quelle
Die Log-Wahrscheinlichkeit für die unkorrelierten Residuen funktionierte einwandfrei, jedoch kann ich die multivariate Normalverteilung nicht herausfinden. Was ist in diesem Fall S? Ich habe S <- vcv.phylo (Baum) ausprobiert und ungefähr -700 für die log-Wahrscheinlichkeit erhalten, während logLik (fit) ungefähr -33 war.
Eric
Entschuldigung - Ich habe es durcheinander gebracht, als ich den Code kopiert habe. Jetzt ist es fertig. S ist die Varianz-Kovarianz-Matrix der Residuen. Sie waren auf dem richtigen Weg (mit der vcvFunktion) - aber Sie müssen die Korrelationsmatrix erhalten und dann verwendenσ^2um dies in die var-cov Matrix umzuwandeln.
Wolfgang,