Calcular "a mano" log-verosimilitud de la regresión no lineal de mínimos cuadrados generalizados (nlme)

Question

Estoy tratando de calcular el logaritmo de la probabilidad para un generalizada no lineal de mínimos cuadrados de la regresión de la función de $f(x)=\frac{\beta_1}{(1+\frac x\beta_2)^{\beta_3}}$ optimizado por el gnls función en el paquete de R nlme, el uso de la varianza de la matriz de covarianza generados por las distancias en un árbol filogenético suponiendo que el movimiento Browniano (corBrownian(phy=tree) de la ape paquete). El siguiente reproducible R código se ajusta a la gnls modelo de uso de la x,y datos y un árbol al azar con 9 taxones:

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) 

Quiero calcular el logaritmo de la probabilidad "a mano" (en R, pero sin el uso de la logLik función), basado en la estimación de los parámetros obtenidos a partir de gnls para que coincida con la salida de logLik(fit). NOTA: no estoy tratando de estimar los parámetros; sólo quiero calcular la log-verosimilitud de los parámetros estimados por el gnls función (aunque si alguien tiene una reproducible ejemplo de cómo estimar los parámetros sin gnls, yo estaría muy interesado en verlo!).

No estoy realmente seguro de cómo ir sobre hacer esto en R. El álgebra lineal notación descrita en Modelos de Efectos Mixtos en S y S-Plus (Pinheiro y Bates) está muy por encima de mi cabeza y ninguno de mis intentos han igualado logLik(fit). Aquí están los detalles descritos por Pinheiro y Bates:

La log-verosimilitud para la generalización de la no lineal de mínimos cuadrados del modelo de $y_i=f_i(\phi_i,v_i)+\epsilon_i$ donde $\phi_i=A_i\beta$ se calcula de la siguiente manera:

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

donde $N$ es el número de observaciones, y $f_i^*(\beta)=f_i^*(\phi_i,v_i)$.

$\Lambda_i$ es positiva definida, $y_i^*=\Lambda_i^{-T/2}y_i$ $f_i^*(\phi_i,v_i)=\Lambda_i^{-T/2}f_i(\phi_i,v_i)$

Fijo $\beta$$\lambda$, el estimador ML de $\sigma^2$ es

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

y el perfil de la log-verosimilitud es

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

que se utiliza con una de Gauss-Seidel algoritmo para encontrar el ML estimaciones de $\beta$$\lambda$. Una menor estimación sesgada de $\sigma^2$ se utiliza:

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

donde $p$ representa la longitud de $\beta$.

He compilado una lista de preguntas específicas que estoy mirando:

1. ¿Qué es $\Lambda_i$? Es la distancia de la matriz producida por big_lambda <- vcv.phylo(tree) en ape, o qué se necesita para ser transformado de alguna manera o con parámetros por $\lambda$ o algo completamente distinto?
2. Sería $\sigma^2$ fit$sigma^2, o la ecuación de la menor estimación sesgada (la última ecuación en este post)?
3. Es necesario el uso de $\lambda$, para calcular la log-verosimilitud, o es que sólo un paso intermedio para la estimación de parámetros? También, cómo es $\lambda$ utiliza? Es un único valor o un vector, y es multiplicado por todos los de $\Lambda_i$ o simplemente fuera de la diagonal de los elementos, etc.?
4. ¿Qué es $||y-f(\beta)||$? Tendría que ser norm(y-f(fit$coefficients,x),"F") en el paquete Matrix? Si es así, estoy confundido acerca de cómo calcular la suma de $\sum\limits_{i=1}^M||y_i^*-f_i^*(\beta)||^2$ porque norm() devuelve un único valor, no un vector.
5. Cómo se calculan las $\log|\Lambda_i|$? Es log(diag(abs(big_lambda))) donde big_lambda es $\Lambda_i$ o es logm(abs(big_lambda)) del paquete expm? Si es logm(), ¿cómo hace uno para tomar la suma de una matriz (o es implícita de que es sólo los elementos de la diagonal)?
6. Sólo para confirmar, es $\Lambda_i^{-T/2}$ calcula de la siguiente manera: t(solve(sqrtm(big_lambda)))?
7. ¿Cómo se $y_i^*$ $f_i^*(\beta)$ calculado? Es cualquiera de los siguientes:

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

y

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

o sería

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

y

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

Si todas estas preguntas son respondidas, en teoría, creo que el logaritmo de la probabilidad debe ser calculable para que coincida con la salida de logLik(fit). Cualquier ayuda en cualquiera de estas preguntas, sería muy apreciado. Si algo necesita aclaración, por favor hágamelo saber. Gracias!

ACTUALIZACIÓN: he estado experimentando con diversas posibilidades para el cálculo de la log-verosimilitud, y aquí es el mejor que han llegado tan lejos. logLik_calc es constantemente alrededor de 1 a 3 de descuento en el valor devuelto por logLik(fit). Yo estoy cerca de la solución real, o esto es puramente por casualidad. Los pensamientos?

  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

Answer 1

Vamos a empezar con el caso más sencillo, donde no hay estructura de correlación de los residuos:

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

El registro de la probabilidad puede ser fácilmente calculada a mano con:

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

Puesto que los residuos son independientes, podemos utilizar dnorm(..., log=TRUE) para obtener el registro individual probabilidad de términos (y luego sumarlos). Alternativamente, podemos usar:

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

Tenga en cuenta que fit$sigma no es el "menos estimación sesgada de $\sigma^2$", por lo que necesita para hacer la corrección manual de la primera.

Ahora, para el caso más complicado donde los residuos se encuentran correlacionados:

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

Aquí, necesitamos utilizar la distribución normal multivariante. Estoy seguro de que hay una función para este lugar, pero vamos a hacer esto a mano:

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)