Simulando y analizando datos

La mejor forma de saber si un método de análisis de datos funciona bien es simular datos con parámetros conocidos y luego aplicar el análisis en cuestión para ver si es capaz de recuperar los valores correctos. De esta forma, tenemos control sobre el origen de los datos y conocemos los valores verdaderos de los parámetros. Esta situación dista mucho de lo que ocurre cuando analizamos datos de campo ya que en esos casos, el modelo que genera datos es desconocido.

Para familiarizarnos con el análisis de datos usando máxima verosimilitud vamos a simular la parte estocástica de los datos usando funciones r y luego analizar esos datos usando las correspondientes funciones d. Por ejemplo, para modelar la remoción de frutos en una planta en un lapso de tiempo definido (e.g. un día) podemos asumir que el proceso de remoción es parecido al que resulta en una distribución Binomial, es decir que cada fruto en una planta tiene la misma probabilidad de ser removido, independientemente de lo que ocurre con los otros frutos. Entonces, para simular la remoción de frutos en una planta que tiene 100 frutos donde la probabilidad de remoción por fruto por día es θ = 0.2 hacemos:

set.seed(1234)
removidos <- rbinom(1, size = 100, prob = 0.2)
removidos

[1] 15

Si queremos ver la probabilidad de observar ese valor de número de frutos removidos dado que hay 100 frutos disponibles y con θ = 0.2 hacemos:

dbinom(removidos, size = 100, prob = 0.2)

## [1] 0.04806179

Es más ilustrativo ver cómo cambia la probabilidad de observar ese valor cuando cambiamos el valor de θ.

pvec <- seq(0, 1, by = 0.01)
op <- par(cex.lab = 1.3, font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(pvec, dbinom(removidos, size = 100, prob = pvec), type = "l", xlab = expression(theta), 
    ylab = "Probabilidad", lwd = 3)
abline(v = removidos/100, lty = 2, lwd = 2)

par(op)

Vemos que "el valor más probable" para θ es 0.15 pero que el valor verdadero de θ = 0.2 está dentro de los posibles. Este ejemplo ilustra un aspecto importante del análisis de datos en general: como trabajamos con muestras, las mejores estimaciones no son necesariamente las correctas. Una cuestión importante entonces es como cuantificar la incertidumbre alrededor de nuestras estimaciones.

Una opción para definir intervalos de confianza es buscar un corte en escala de − log likelihood, por ejemplo dos unidades más grandes que el mínimo. El gráfico de abajo muestra el mínimo de negative log likelihood y los puntos de corte. Allí se puede ver que el valor verdadero queda incluido dentro de esta especie de intervalo de confianza.

Comunmente se usa una prueba de likelihood ratio test (LRT) para definir intervalos de confianza a partir de perfiles de likelihood como los de la figura de arriba. Para hacer un LRT consideramos una función de likelihood con n parámetros y encontramos el mejor ajuste, llamamos a este valor de likelihood como ℒ_a con "a" de absoluto. Ahora, hacemos otro ajuste fijando una cantidad f de parámetros y maximizamos el ajuste cambiando los otros. Llamamos a este nuevo valor de likelihood ℒ_r con "r" de restringido. El LRT se basa en resultados teóricos que dicen que para muestras grandes, la distribución de − 2logℒ_r/ℒ_a, es decir la devianza, sigue una distribución χ_f² con grados de libertad igual al número de parámetros fijos. Entonces, para límites de confianza univariados el corte se hace en min(neg log lik + χ₁²(1 − α)/2, donde α es el nivel de confianza seleccionado (dividimos por dos porque es un test de una sola vía). El valore de χ₁²(1 − α)/2 para un intervalo de 95% es 1.9207294 que es un poco menor que el 2 que usamos antes.

pvec <- seq(0, 1, by = 0.001)
nll <- -dbinom(removidos, size = 100, prob = pvec, log = TRUE)
op <- par(cex.lab = 1.3, font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(pvec, nll, type = "l", ylim = c(0, 20), xlab = expression(theta), ylab = "- log lik")
abline(h = min(nll) + qchisq(0.95, 1)/2, lty = 2)
tmp <- which.min(nll)
points(pvec[tmp], nll[tmp], pch = 16)
izq <- nll[1:tmp]
der <- nll[tmp:length(nll)]
tmp1 <- which(abs(izq - (min(nll) + qchisq(0.95, 1)/2)) == min(abs(izq - (min(nll) + 
    qchisq(0.95, 1)/2))))
points(pvec[tmp1], nll[tmp1], pch = 16, col = 2)
abline(v = pvec[tmp1], lty = 2)
tmp2 <- which(abs(der - (min(nll) + qchisq(0.95, 1)/2)) == min(abs(der - (min(nll) + 
    qchisq(0.95, 1)/2))))
points(pvec[tmp + tmp2], nll[tmp + tmp2], pch = 16, col = 2)
abline(v = pvec[tmp + tmp2], lty = 2)
abline(v = 0.2)

Como vimos en la clase pasada, la forma en la que la curva de negative-log-likelihood se aleja del óptimo depende del número de observaciones. Veamos esto nuevamente con datos simulados:

set.seed(1234)
pvec <- seq(0, 1, by = 0.001)
removidos <- rbinom(10, size = 100, prob = 0.2)
nll <- numeric(length(pvec))
for (i in 1:length(pvec)) {
    nll[i] <- -sum(dbinom(removidos, size = 100, prob = pvec[i], log = TRUE))
}
op <- par(cex.lab = 1.3, font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(pvec, nll, type = "l", xlab = expression(theta), ylab = "neg log lik", 
    ylim = c(0, 200))

abline(h = min(nll) + qchisq(0.95, 1)/2, lty = 2)
tmp <- which.min(nll)
points(pvec[tmp], nll[tmp], pch = 16)
izq <- nll[1:tmp]
der <- nll[tmp:length(nll)]
tmp1 <- which(abs(izq - (min(nll) + qchisq(0.95, 1)/2)) == min(abs(izq - (min(nll) + 
    qchisq(0.95, 1)/2))))
points(pvec[tmp1], nll[tmp1], pch = 16, col = 2)
abline(v = pvec[tmp1], lty = 2)
tmp2 <- which(abs(der - (min(nll) + qchisq(0.95, 1)/2)) == min(abs(der - (min(nll) + 
    qchisq(0.95, 1)/2))))
points(pvec[tmp + tmp2], nll[tmp + tmp2], pch = 16, col = 2)
abline(v = pvec[tmp + tmp2], lty = 2)
abline(v = 0.2)

par(op)

Vemos que con más observaciones, el valore de MLE se parece más al verdadero (0.2) y el intervalo de confianza es menor. Sin embargo, esperamos que si repetimos este proceso muchas veces, aproximadamente el 5% de las veces el intervalo de confianza calculado de esta manera no va a incluir al valor verdadero. Podemos comprobar esto haciendo simulaciones. Para simplificar un poco las cosas vamos a ajustar el modelo con mle2 para calcular los intervalos de confianza de manera más fácil.

library(bbmle)
nrep <- 500
n <- 10
cmin <- numeric(nrep)
cmax <- numeric(nrep)

binomNLL <- function(p, k, N) {
    -sum(dbinom(k, prob = p, size = N, log = TRUE))
}

cuenta <- 0

for (i in 1:nrep) {
    removidos <- rbinom(n, size = 100, prob = 0.2)
    fit <- mle2(minuslogl = binomNLL, start = list(p = 0.5), data = list(N = 100, 
        k = removidos))
    ci <- confint(fit)
    cmin[i] <- ci[1]
    cmax[i] <- ci[2]
    if (0.2 > cmin[i] & 0.2 < cmax[i]) 
        cuenta <- cuenta + 1
}

cuenta/nrep

[1] 0.946

op <- par(cex.lab = 1.3, font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(1:nrep, rep(0.2, nrep), type = "l", xlab = "réplica", ylab = "ci")
arrows(1:nrep, cmin, 1:nrep, cmax, length = 0)
lines(1:nrep, rep(0.2, nrep), type = "l", lwd = 3, col = 2)

par(op)

Pruebas de Potencia

También podemos usar simulaciones para ver la potencia de nuestros experimentos. Supongamos por ejemplo que queremos cuantificar el efecto del número de frutos en una planta sobre la remoción. Vamos a simular el proceso que nos interesa medir y cuántas muestras vamos a tomar. Una de las cosas que tenemos que definir es el tamaño del efecto que suponemos que existe. Por ejemplo, podemos pensar que el efecto de la abundancia de frutos incrementa la remoción de frutos en un 10% en escala de desvío estándar de abundancia de frutos.

library(bbmle)
set.seed(1234)
n <- 20
beta <- 0.1
nreps <- 1000

res <- matrix(NA, nreps, 4)
fnll <- function(pars) {
    beta <- pars[1]
    pr <- plogis(beta * nfset)
    res <- -sum(dbinom(removidos, size = fset, prob = pr, log = TRUE))
}

# definimos nombres para los parámetros
parnames(fnll) <- c("beta")


for (i in 1:nreps) {
    fset <- rpois(n, lambda = 20)  # cantidad de frutos
    nfset <- (fset - mean(fset))/sd(fset)  # cantidad de frutos estandarizada
    pr <- plogis(beta * nfset)  # probabilidad de remoción en función de la abundancia de frutos
    removidos <- rbinom(n, size = fset, prob = pr)  # datos simulados
    fit <- mle2(fnll, start = c(beta = 0), data = list(removidos = removidos, 
        fset = fset))
    
    tmp <- summary(fit)
    ci <- confint(fit)
    res[i, ] <- tmp@coef
}

En la matriz res, tenemos en la primera columna al valor estimado de beta para cada réplica, en la segunda columna el error estándard de beta, en la tercera el valor de z y en la cuarta el valor de p para la hipótesis nula de que beta es igual a cero. Podemos ver por ejemplo cuántas veces encontramos valores significativos para beta. Para esto hacemos length(which(res[,4]<0.05))/nreps * 100, que nos da 17.5. Vemos entonces que una muestra de n = 20 resulta en estimaciones de efecto significativo del número de frutos sobre la tasa de remoción solamente en un porcentaje bajo de las veces. Más preocupante todavía es ver qué valores tiene beta cuando es significativo:

hist(res[which(res[, 4] < 0.05), 1], 30, main = "", xlab = "beta", freq = FALSE)

Esto es lo que nuestro amigo Andrew Gelman llama errores tipo M (de magnitud) y tipo S (de signo). Mete miedo...

sessionInfo()

## R version 3.2.2 (2015-08-14)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.10.5 (Yosemite)
## 
## 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.17 knitr_1.11  
## 
## loaded via a namespace (and not attached):
##  [1] lattice_0.20-33      digest_0.6.8         MASS_7.3-44         
##  [4] mime_0.4             grid_3.2.2           knitrBootstrap_1.0.0
##  [7] formatR_1.2          magrittr_1.5         evaluate_0.7.2      
## [10] stringi_0.5-5        rmarkdown_0.8        tools_3.2.2         
## [13] stringr_1.0.0        markdown_0.7.7       numDeriv_2014.2-1   
## [16] yaml_2.1.13          htmltools_0.2.6

Juan Manuel Morales. 6 de Septiembre del 2015. Última actualización: 2015-09-28