Máxima Verosimilitud (Likelihood)
Objetivos:
- Relacionar el likelihood de datos con parámetros de una distribución
- Familiarizarse con curvas de likelihood y log-likelihood
- Comprobar que la forma de la curva de likelihood depende del tamaño muestreal
- Aprender a usar herramientas para encontrar óptimos
Vimos que las funciones de densidad nos dan la probabilidad de observar un valor determinado de una variable aleatoria (para variables continuas nos dan la probabilidad alrededor de un valor). Por ejemplo, para una distribución de Poisson, la probabilidad para un valor \(y\) es \(\lambda^{y} \exp(- \lambda) / y!\) donde \(!\) es el factorial de \(y\). Entonces, para una variable aleatoria que asumimos tiene una distribución de Poisson con \(\lambda=2\) la probabilidad de observar un cero (\(y=0\)) es de \(2^{0} \exp(-2)/0! = \exp(-2)\). En R usamos la función dpois para hacer este tipo de cuentas. Para este ejemplo escribimos dpois(0, lambda=2).
Normalmente tenemos un conjunto de datos, por ejemplo números de renovales por metro cuadrado de muestreo en un bosque: \(y = 1, 0, 2, 1, 0, 0, 1, 1, 2, 2\). Si asumimos que estos datos se pueden describir como el resultado de un proceso aleatorio de tipo Poisson, y que además las observaciones son independientes, podemos calcular la probabilidad del conjunto de datos como el producto de las probabilidades de las observaciones:
\[ p(y)= \prod_i^n{ \frac{\lambda^{y_i} \exp(- \lambda)} {y_i!} } \]
Podemos ver que el producto (\(\prod\)) se calcula usando \(i\) como índice para cada observación del set de datos \(y\), desde \(i=1\) hasta \(i=n\) que es el tamaño de \(y\). Esta probabilidad de observar los datos es lo que llamamos likelihood o verosimilitud. Para hacer la cuenta de arriba tenemos que definir además el parámetro \(\lambda\). La probabilidad de los datos depende entonces no sólo de los valores presentes en el set de datos sino también de la distribución utilizada y el valor del o los parámetros de la distribución. De manera más general, la expresión es un modelo de datos, es decir una descripción probabilística de como llegar a los datos. Más adelante veremos como combinando procesos determinísticos y estocásticos podemos armar modelos de datos realistas y satisfactorios.
Veamos ahora como calcular la probabilidad de estos datos en R asumiendo \(\lambda = 1\):
y <- c(1, 0, 2, 1, 0, 0, 1, 1, 2, 2)
lambda <- 1
p.y <- dpois(y, lambda = lambda)
p.y ## [1] 0.3678794 0.3678794 0.1839397 0.3678794 0.3678794 0.3678794 0.3678794
## [8] 0.3678794 0.1839397 0.1839397prod(p.y)## [1] 5.674991e-06Como las probabilidades son números que van entre \(0\) y \(1\), el producto de varias probabilidades resulta en números muy pequeños. Vemos en este caso que la probabilidad del conjunto de datos es muy baja (5.67e-06, donde e-06 es \(10^{-6}\) = \(0.000001\)) pero esto no debería sorprendernos porque esa es la probabilidad de encontrar exactamente esos datos. Entonces, para evitar problemas numéricos preferimos trabajar con la suma de los logaritmos de las probabilidades, es decir con el logaritmo del likelihood. En R, podemos usar la opción log = TRUE dentro de las funciones de densidad.
y <- c(1, 0, 2, 1, 0, 0, 1, 1, 2, 2)
lambda <- 1
ll.y <- dpois(y, lambda = lambda, log = TRUE)
sum(ll.y)## [1] -12.07944Una vez que logramos hacer estas cuentas, podemos preguntarnos cómo cambia la probabilidad de observar el set de datos si cambiamos el valor de \(\lambda\).
Probemos con \(\lambda=2\):
y <- c(1, 0, 2, 1, 0, 0, 1, 1, 2, 2)
lambda <- 2
sum(dpois(y, lambda = lambda, log = TRUE))## [1] -15.14797Con un \(\lambda\) de \(2\) la log likelihood de los datos empeora, pero podemos seguir probando distintos valores de \(\lambda\) hasta encontrar el que hace más posible al set de datos. Esta es la idea detrás de los análisis de máxima verosimilitud (a falta de mejor nombre) o de maximum likelihood (si estamos dispuestos a usar palabras en inglés). Normalmente tenemos más de un parámetro en nuestro modelo de datos, y usamos aglún algoritmo de búsqueda para encontrar la combinación de parámetros que hace más probables a nuestros datos (por tradición, los algoritmos buscan minimizar el menos logaritmo de la verosimilitud de nuestros datos). Esta idea de buscar combinaciones de parámetros que maximizan la verosimilitud de los datos se conoce como maximum likelihood estimation y muchas veces aparece abreviado como MLE.
Ejemplo: Remoción de Frutos
Mariano Rodriguez-Cabal estudió la remoción de frutos de quintral por parte del monito del monte en bosque continuo y en fragmentos en 3 sitios pareados con la idea de ver si la fragmentación afectabla la interacción entre este marsupial y el muérdago. El trabajo está dispoible [aquí]((https://github.com/jmmorales/cursoMyD/raw/master/papers/Rodriguez-Cabal_07.pdf)
Los datos son de número de frutos removidos en la temporada reproductiva de un total disponible. Podemos usar una distribución Binomial para capturar la variabilidad (estocasticidad) en el proceso de remoción de frutos. Es decir, asumimos que la distribución Binomial tiene sentido para representar el proceso que genera los datos. Dicho de otro modo, asumimos que los datos son una variable aleatoria con distribución Binomial.
url <- "https://github.com/jmmorales/cursoMyD/raw/master/Data/quintral.txt"
quintral <- read.table(url, header = TRUE)Para ver como están organizados los datos podemos hacer:
str(quintral) ## 'data.frame': 70 obs. of 4 variables:
## $ bosque : chr "c" "c" "c" "c" ...
## $ sitio : chr "Campanario" "Campanario" "Campanario" "Campanario" ...
## $ Removidos: int 30 40 11 33 35 25 48 13 25 31 ...
## $ Frutos : int 46 51 32 49 46 37 53 31 33 44 ...Vemos que se trata de un data frame que es la manera preferida de R para organizar datos. Los data frames tienen las observaciones en filas y las variables en columnas. Las variables individuales se pueden acceder escribiendo el nombre del data frame seguido del signo $ y el nombre de la variable. Por ejemplo, si queremos ver los valores de frutos disponibles en las tres primeras observaciones hacemos:
quintral$Frutos[1:3]## [1] 46 51 32Ahora calculemos el likelihood para la primera observación del set de datos asumiendo que la tasa de remoción es \(0.5\) (si no nos acordamos cómo es la distribución Binomial o cómo usarla en R podemos pedir ayuda con ?dbinom). En lo que sigue, vamos a llamar \(\theta\) a la tasa (probabilidad) de remoción de frutos.
theta <- 0.5
dbinom(quintral$Removidos[1], size=quintral$Frutos[1], prob=theta)## [1] 0.01408998Normalmente trabajamos con log-likelihoods:
dbinom(quintral$Removidos[1], size=quintral$Frutos[1], prob=theta, log=TRUE)## [1] -4.262292¿Cómo cambia la probabilidad de observar ese dato cuando cambiamos el parámetro de probabilidad de éxito (remoción) en la Binomial? Como en este caso el likelihood depende de un solo parámetro, podemos responder a esta pregunta gráficamente.
op <- par(mfrow=c(1,2),lwd=2, bty="l", las=1)
theta <- seq(0, 1, length = 1000)
plot(theta, dbinom(quintral$Removidos[1], size=quintral$Frutos[1], prob=theta, log=FALSE),
type="l", ylab="L", xlab = expression(theta))
plot(theta,-dbinom(quintral$Removidos[1], size=quintral$Frutos[1], prob=theta, log=TRUE),
type="l", ylab="-logL", xlab = expression(theta))
par(op)También podemos ver cómo cambia el Likelihood cuando cambiamos el número de observaciones que consideramos.
op<-par(mfrow = c(2,3), lwd=2, bty="l", cex.lab=1.5, las = 1, mar = c(5,4,1,1))
theta <- seq(0, 1, length = 1000)
plot(theta, dbinom(quintral$Removidos[1], size=quintral$Frutos[1], prob=theta,log=FALSE),
type="l",xlab="", ylab="L", main="n = 1")
nLL2 <- array(0, length(theta))
for (i in 1:length(theta)) {
nLL2[i] <- -sum(dbinom(quintral$Removidos[1:2], size=quintral$Frutos[1:2],
prob=theta[i],log=TRUE))
}
plot(theta,exp(-nLL2), type = "l", xlab = "", ylab="", main="n=2")
nLL10 <- array(0, length(theta))
for (i in 1:length(theta)) {
nLL10[i] <- -sum(dbinom(quintral$Removidos[1:10], size=quintral$Frutos[1:10],
prob=theta[i],log=TRUE))
}
plot(theta,exp(-nLL10), type = "l",xlab = "", ylab = "", main="n=10")
plot(theta, -dbinom(quintral$Removidos[1], size=quintral$Frutos[1], prob=theta,log=TRUE),
type="l",xlab= "", ylab="-logL", main="")
plot(theta, nLL2, type = "l", xlab=expression(theta), ylab="", main="", ylim=c(0,100))
plot(theta, nLL10, type = "l", xlab="", ylab="", main="", ylim=c(0,100))
par(op)Pregunta (1): ¿Cómo se interpretan estos cambios?
Buscando el máximo (o más bien el mínimo…)
Queremos encontrar el valor de tasa de remoción por fruto que maximiza el likelihood (o minimiza el negative log-likelihood) del conjunto de datos. Una opción es hacerlo a lo bruto:
theta <- seq(0, 1, length = 50)
nLL <- array(NA, length(theta))
for (i in 1:length(theta)) {
nLL[i] <- -sum(dbinom(quintral$Removidos, size=quintral$Frutos, prob=theta[i], log=TRUE))
}
theta[which(nLL == min(nLL))]## [1] 0.5918367En este caso, existe una solución analítica para la tasa de remoción y podemos compararla con la que acabamos de calcular:
sum(quintral$Removidos)/sum(quintral$Frutos)## [1] 0.5968072Pregunta (2): ¿Qué pasa si cambiamos la resolución del vector de tasas de remoción theta <- seq(0, 1, length = 50) ?
Búsquedas más eficientes usando optim y mle2 de Ben Bolker
Claramente no podemos recurrir a búsquedas gráficas o de fuerza bruta cuando tenemos modelos con varios parámetros. Para esos casos, podemos hacer uso de los algoritmos de búsqueda que están disponibles en R, por ejemplo a través de la función optim. Primero tenemos que definir una función para “optimizar”, (para optim optimizar viene a ser encontrar el mínimo). Además, tenemos que definir por dónde empezar a buscar y con qué método.
binNLL <- function(theta,k,N){
-sum(dbinom(k, prob=theta, size=N, log=TRUE))
}
O1 <- optim(fn=binNLL, par=list(theta = 0.5), N=quintral$Frutos, k=quintral$Removidos, method="BFGS")Cuando ejecutamos esta función, en este caso aparecen algunos “warnings” de NaNs producidos en dbinom pero nada de qué preocuparse (esto ocurre cuando optim prueba algún valor de tasa de remoción que no tiene sentido, por ejemplo -0.1). El objeto de salida de optim, que en este caso se llama O1 (por falta de creatividad) contiene:
- en
$par, el valor óptimo del parámetro “p” - en
$value, el valor del log likelihood para $par - en
$counts, hay algunas cosas crípticas acerca de cómo se desarrolló la búsqueda
Antes de poder decir algo acerca de los valores de los parámetros, hay que asegurarse de que el algoritmo haya convergido ( $convergence = 0 )
La función optim en R es una función general para optimización, pero para MLE nos conviene usar mle2 desarrollada por Ben Bolker y parte del paquete bbmle. En mle2 la función de negative log-likelihood se llama “minuslogl” en vez de “fn” (en optim).
library(bbmle)
m1 <- mle2(minuslogl=binNLL, start = list(theta = 0.5), data = list(N=quintral$Frutos, k=quintral$Removidos), method="BFGS")En m1 tenemos la función usada, los valores de los coeficientes y el log-likelihood estimado. Podemos ver más detalles con summary(m1)
summary(m1)## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = binNLL, start = list(theta = 0.5), method = "BFGS",
## data = list(N = quintral$Frutos, k = quintral$Removidos))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## theta 0.5967997 0.0099246 60.133 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 717.3864Ahora tenemos errores standares y valores de ‘p’ para los parámetros (coeficientes), poniendo a prueba la hipótesis nula de que el coeficiente es igual a cero. También aparece la devianza (-2 log-likelihood). Más adelante veremos como usar éstas y otras salidas de mle2 para obtener intervalos de confianza y hacer comparaciones entre modelos.
Comparación (cruda) con los datos
Una vez que ajustamos el modelo es natural preguntarse qué tan bueno es el modelo para representar el proceso ecológico que nos interesa. Una forma rústica pero efectiva para hacer esta evaluación es graficar la distibución de valores observados con los predichos por el modelo:
k <- quintral$Removidos
hist(k, breaks = 40, freq=FALSE, ylim = c(0,0.15), main = "")
points(0:50, dbinom(0:50, size=round(mean(quintral$Frutos)), prob=m1@coef),
pch=16, col="darkgray")
Una limitación de la comparación que acabamos de hacer es que el número de frutos disponibles es variable. Para hacer una mejor comparación podemos usar un bootstrap (remuestreo)
B <- 1000
res <- matrix(0, length(quintral$Frutos), B)
for (i in 1:B){
n <- sample(quintral$Frutos, length(quintral$Frutos), replace=TRUE)
res[,i] <- rbinom(length(quintral$Frutos), size=n, prob=m1@coef)
}
hist(k, breaks = 40, freq=F, ylim = c(0,0.15), main = "")
lines(density(res), lwd=2)
Una alternativa más estándar es graficar los residuos (\(y - \hat{y}\)):
plot(m1@coef * quintral$Frutos - quintral$Removidos, xlab = "", ylab = "Residuos")
abline(h=0, lwd=2, lty=2)
Pregunta (3): ¿Qué se puede decir sobre el ajuste del modelo a los datos?
Ejercicio
- Estimar las tasas de remoción de frutos por sitio y tipo de bosque. Pista: pueden usar la función
whichpara separar los datos correspondientes a cada caso.
Ejemplo: Deposición de polen en flores de pomelo
Natacha Chacoff estudió la deposición de polen en flores de pomelo en función del número de visitas de abejas. El trabajo terminado se puede ver aquí
apis <- read.table("https://github.com/jmmorales/cursoMyD/raw/master/Data/polen_Apis.csv", header = TRUE, sep = ",")
plot(apis$visitas, apis$granos, xlab="Número de Visitas",
ylab="Granos de Polen Depositados", pch=19, cex.lab = 1.2)
title(main= expression(paste("Eficiencia de ", italic(Apis~mellifera))))
Podemos tratar de descibir esta relación entre visitas y deposición de polen como una función lineal asumiendo una distribución de Poisson para los granos de polen. Una forma de escribir este modelo es:
\[ \begin{aligned} & y_i \sim \textrm{Poisson}(\lambda_i) \\ & \lambda_i = \beta_0 + \beta_1 \times v_i\\ & \textrm{para } i \textrm{ de } 1 \textrm{ a } n\\ \end{aligned} \]
Donde \(y\) son las observaciones (número de granos de polen depositados) y \(v\) la variable predictora o la “covariable” (número de visitas de abejas) y los subíndices \(i\) indican cada una de las \(n\) observaciones. El script que sigue define y ajusta este modelo usando mle2.
m2 <- function (pars) {
b0 <- pars[1]
b1 <- pars[2]
lambda <- b0 + b1 * visitas
-sum(dpois(granos, lambda, log=TRUE))
}
parnames(m2) <- c("b0", "b1")
fit.m2 <- mle2(m2, start = c(b0 = 10, b1 = 3),
data = list(visitas = apis$visitas, granos = apis$granos))
summary(fit.m2)## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = m2, start = c(b0 = 10, b1 = 3), data = list(visitas = apis$visitas,
## granos = apis$granos))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## b0 14.87137 0.36203 41.078 < 2.2e-16 ***
## b1 6.85330 0.26371 25.988 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 5416.507plot(apis$visitas, apis$granos,
xlab="Número de Visitas",
ylab="Granos de Polen Depositados")
title(main= expression(paste("Eficiencia de ", italic(Apis~mellifera))))
x = 0:9 # vector de valores de referencia
lines(x, fit.m2@coef[1] + x * fit.m2@coef[2], lwd=3)
Pregunta (4): ¿Qué inconvenientes podríamos enocontrar tratando de ajustar este modelo asumiendo que las observaciones tienen una distribución de Poisson?
Pensando un poco, quizás sea esperable que en algún momento la cantidad de polen por estigma se sature ya que no podemos amontonar infinitos granos de polen. Cambiar el modelo anterior a un modelo de saturación no es muy complicado:
m3 <- function (pars) {
y0 <- pars[1]
a <- pars[2]
b <- pars[3]
lambda <- y0 + a*(1-exp(-b*visitas))
-sum(dpois(granos, lambda, log=TRUE))
}
parnames(m3) <- c("y0", "a", "b")
fit.m3 <- mle2(m3, start = c(y0 = 10, a = 41, b = 0.16),
data = list(visitas = apis$visitas, granos = apis$granos))
summary(fit.m3)## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = m3, start = c(y0 = 10, a = 41, b = 0.16), data = list(visitas = apis$visitas,
## granos = apis$granos))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## y0 14.091487 0.363584 38.7572 < 2.2e-16 ***
## a 40.857524 4.488943 9.1018 < 2.2e-16 ***
## b 0.291066 0.052784 5.5143 3.502e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 5363.868plot(apis$visitas, apis$granos,
xlab="Número de Visitas",
ylab="Granos de Polen Depositados",
pch=19)
title(main= expression(paste("Eficiencia de ", italic(Apis~mellifera))))
lines(x, fit.m3@coef[1] + fit.m3@coef[2]*(1-exp(-fit.m3@coef[3]* x)), lwd=3)
Pregunta (5): ¿Cómo ajustan estos modelos a los datos?
Pregunta (6): ¿Qué otros modelos podríamos probar?
sessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] bbmle_1.0.24 knitr_1.36
##
## loaded via a namespace (and not attached):
## [1] magrittr_2.0.1 MASS_7.3-54 lattice_0.20-45
## [4] R6_2.5.1 rlang_0.4.12 fastmap_1.1.0
## [7] stringr_1.4.0 highr_0.9 tools_4.1.2
## [10] grid_4.1.2 xfun_0.28 jquerylib_0.1.4
## [13] htmltools_0.5.2 yaml_2.2.1 digest_0.6.28
## [16] numDeriv_2016.8-1.1 Matrix_1.3-4 bdsmatrix_1.3-4
## [19] sass_0.4.0 evaluate_0.14 rmarkdown_2.11
## [22] stringi_1.7.5 compiler_4.1.2 bslib_0.3.1
## [25] jsonlite_1.7.2 mvtnorm_1.1-3Juan Manuel Morales. 6 de Septiembre del 2015. Última actualización: 2021-12-02