Máxima Verosimilitud (Likelihood)

Objetivos:

  • Aprender a usar herramientas para encontrar óptimos

Búsquedas más eficientes usando “optim” y “mle2” de Ben Bolker

Vimos en la introducción a MLE en qué consiste el likelihood de un set de datos y cómo podemos buscar combinaciones de parámetros que maximicen la probabilidad de nuestros datos. Cuando tenemos modelos con varios parámetros, claramente no podemos recurrir a búsquedas gráficas o de fuerza bruta como en los ejemplos que vimos antes. Entonces, tenemos que hacer uso de los algoritmos de búsqueda que están disponibles en R, por ejemplo con la función optim. Primero tenemos que definir una función para “optimizar”. Para la función optim, optimizar quiere decir encontrar el mínimo. Además, tenemos que definir por dónde empezar a buscar y con qué método. Veamos por ejemplo cómo haríamos con el ejemplo de la remoción de frutos de quintral:

url <- "https://sites.google.com/site/modelosydatos/quintral.txt"
quintral <- read.table(url, header = TRUE)

# escribimos una función para el - log likelihood
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 “theta”
  • 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

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

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), 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.3864

Ahora tenemos errores standares y valores de ‘p’ para los coeficientes (poniendo a prueba la hipótesis nula de que el coeficiente es igual a cero, cosa que en este caso no tiene mucho sentido). 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 = 21, bg = "gray")

Una limitación de la comparación que acabamos de hacer es que graficamos la distribución de valores de frutos removidos para el valor promedio de frutos disponibles, de manera que estamos subestimando la variablidad en los datos que se debe a que el número de frutos disponibles en las plantas muestreadas es variable. Para hacer una mejor comparación podemos usar un bootstrap (remuestreo):

nboot <- 1000
res <- matrix(0, length(quintral$Frutos), nboot)

for (i in 1:nboot) {
    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 = FALSE, ylim = c(0, 0.15), main = "")
lines(density(res), lwd = 2)

Pregunta (1): ¿Qué se puede decir sobre el ajuste del modelo a los datos?

Ejercicio

  1. Estimar las tasas de remoción de frutos por sitio y tipo de bosque. Pista: pueden usar la función which para 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://sites.google.com/site/modelosydatos/polen_Apis.csv", 
    header = T, sep = ",")

plot(apis$visitas, apis$granos, xlab = "Número de Visitas", ylab = "Granos de Polen Depositados", 
    pch = 21, bg = "gray", cex.lab = 1.3, las = 1)
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.

poislin <- function(pars) {
    b0 <- pars[1]
    b1 <- pars[2]
    lambda <- b0 + b1 * visitas
    -sum(dpois(granos, lambda, log = TRUE))
}

parnames(poislin) <- c("b0", "b1")

fit.lin <- mle2(poislin, start = c(b0 = 10, b1 = 3), data = list(visitas = apis$visitas, 
    granos = apis$granos))

summary(fit.lin)
## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = poislin, 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.507
plot(apis$visitas, apis$granos, xlab = "Número de Visitas", ylab = "Granos de Polen Depositados", 
    pch = 21, bg = "gray", cex.lab = 1.3, las = 1)
title(main = expression(paste("Eficiencia de ", italic(Apis ~ mellifera))))

x = 0:9  # vector de valores de referencia
lines(x, fit.lin@coef[1] + x * fit.lin@coef[2], lwd = 3)

Pregunta (2): ¿Qué inconvenientes podríamos enocontrar tratando de ajustar este modelo asumiendo que las observaciones tienen una distribución de Poisson? ***

Pensando un poco sobre modelo, 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:

poissat <- 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(poissat) <- c("y0", "a", "b")
fit.sat <- mle2(poissat, start = c(y0 = 10, a = 41, b = 0.16), data = list(visitas = apis$visitas, 
    granos = apis$granos))

summary(fit.sat)
## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = poissat, 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.868
plot(apis$visitas, apis$granos, xlab = "Número de Visitas", ylab = "Granos de Polen Depositados", 
    pch = 21, bg = "gray", cex.lab = 1.3, las = 1)
title(main = expression(paste("Eficiencia de ", italic(Apis ~ mellifera))))

lines(x, fit.sat@coef[1] + fit.sat@coef[2] * (1 - exp(-fit.sat@coef[3] * x)), 
    lwd = 3)

Pregunta (3): ¿Cómo ajustan estos modelos a los datos?

Pregunta (4): ¿Qué otros modelos podríamos probar?

sessionInfo()
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.3
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/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.20 knitr_1.20  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.0        bookdown_0.8      lattice_0.20-38  
##  [4] digest_0.6.18     MASS_7.3-51.1     grid_3.5.2       
##  [7] formatR_1.5       magrittr_1.5      evaluate_0.12    
## [10] stringi_1.2.4     rmarkdown_1.10.15 tools_3.5.2      
## [13] stringr_1.3.1     numDeriv_2016.8-1 xfun_0.4         
## [16] yaml_2.2.0        compiler_3.5.2    htmltools_0.3.6

Juan Manuel Morales. 6 de Septiembre del 2015. Última actualización: 2019-03-02