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: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)
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, gnls
sodass 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 gnls
Funktion geschätzten Parameter berechnen (obwohl gnls
ich 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 β wird wie folgt berechnet:
wobei die Anzahl der Beobachtungen ist und f ∗ i ( β ) = f.
und die profilierte Log-Wahrscheinlichkeit ist
Ich habe eine Liste mit spezifischen Fragen zusammengestellt, mit denen ich konfrontiert bin:
big_lambda <- vcv.phylo(tree)
ape
- wäre
fit$sigma^2
, oder die Gleichung für die weniger voreingenommen Schätzung (die letzte Gleichung in diesem Beitrag)? - Ist es notwendig zu verwenden
norm(y-f(fit$coefficients,x),"F")
Matrix
, weilnorm()
ein einzelner Wert zurückgegeben wird, kein Vektor.- Wie rechnet man ? Ist es
log(diag(abs(big_lambda)))
wobig_lambda
ist, oder ist eslogm(abs(big_lambda))
aus dem Paketexpm
? Wenn jalogm()
, wie nimmt man die Summe einer Matrix (oder impliziert man, dass es sich nur um die diagonalen Elemente handelt)? - Nur um zu bestätigen, ist berechnet wie folgt aus :
t(solve(sqrtm(big_lambda)))
? - Wie bist und 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_calc
ist 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
Antworten:
Beginnen wir mit dem einfacheren Fall, in dem es keine Korrelationsstruktur für die Residuen gibt:
Die Log-Wahrscheinlichkeit kann dann einfach von Hand berechnet werden mit:
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:Beachten Sie, dass diesσ2 "- also müssen wir die Korrektur zuerst manuell vornehmen.
fit$sigma
nicht die "weniger voreingenommene Schätzung von" istNun zum komplizierteren Fall, in dem die Residuen korreliert sind:
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:
quelle
vcv
Funktion) - aber Sie müssen die Korrelationsmatrix erhalten und dann verwenden