Estadística

No sé si visteis el vídeo que colgué el otro día. Trataba el problema de determinar si dos poblaciones

beer  <- c(27, 20, 21, 26, 27, 31, 24,
        21, 20, 19, 23, 24,
        18, 19, 24, 29, 18, 20, 17,
        31, 20, 25, 28, 21, 27)
water <- c(21, 22, 15, 12, 21, 16, 19,
        15, 22, 24, 19, 23, 13,
        22, 20, 24, 18, 20)

tienen o no la misma media. Más concretamente, si la población beer tiene una media superior a la de water como en efecto sucede:

mean(beer)
#[1] 23.2
mean(water)
#[1] 19.22222

¿Pero es esta diferencia significativa?

Muchos plantearían un t-test:

t.test(beer, water, alternative = "greater")
# Welch Two Sample t-test
#
# data:  beer and water
# t = 3.3086, df = 39.271, p-value = 0.001007
# alternative hypothesis: true difference in means is greater than 0
# 95 percent confidence interval:
#   1.952483      Inf
# sample estimates:
#   mean of x mean of y
# 23.20000  19.22222

Pero en el vídeo se propone una alternativa basada en remuestreos:

Hace unos años, Juanjo Gibaja y yo organizamos un “curso de estadística moderna con R”. Queríamos mostrar en él que otra estadística es posible, que con la ayuda de los ordenadores (¡y de R!) los problemas clásicos de la estadística pueden afrontarse de otra manera. Y que esta manera es más natural y accesible.

Hoy uno de nuestros antiguos alumnos nos ha agradecido que le señalásemos el camino de esos superpoderes:

El otro día hablaba con una colega sobre una charla a la que habíamos asistido. Yo le decía que sí, que estaba bien, pero que todo lo que habían contado era mentira. Debí haber sido más preciso y decir que no era verdad, que es distinto. Pero las canapescas circunstancias no eran propicias para el distingo. Mi interlocutora me escuchaba, pienso, entre sorprendida e incrédula. Todavía está en la edad en la que hay que creérselo todo —sí, esa edad y esa obligación existe— y tiempo tendrá de dejarse envenenar por el nihilismo. Es lo suficientemente lista como para eso.

Todos hemos comido macarrones con cosas de la nevera. Estás en casa, tienes hambre y, si no hay otra cosa, son estupendos. Distinto es ir a un bodorrio de alto copete y decirle al camarero:

—Oiga, esto del solomillo y tal… ¿No tendrán Vds. un platazo de macarrones con cosas de la nevera?

Viene esto a que cierta gente trabaja con grandes datos. Y quieren construir modelos. Y por algún motivo que no comprendo del todo, optan por la regresión logística. Hay mil motivos por los que estaría desaconsejado ajustar regresiones logísticas con todos los datos. Aun así, hay gente —sí, la hay— que lo hace.

Hoy voy a hablar de esa especie de oxímoron que es el el bootstrap bayesiano. Comenzaré planteando un pequeño problema bien conocido: tenemos números $latex x_1, \dots, x_n$ y hemos calculado su media. Pero nos preguntamos cómo podría variar dicha media (de realizarse otras muestras).

La respuesta de Efron (1979) es esta:

replicate(n, mean(sample(x, length(x), replace = TRUE)))

Es decir, crear muestras de $latex x_i$ con reemplazamiento y hacer la media de cada una de ellas para obtener su presunta distribución (o una muestra de la presunta distribución de esa media).

Me gusta criticar. Bien lo saben quienes me siguen. Pero hoy toca aplaudir un artículo tan raro como valiente. Que no hace sino criticar por mí. Se titula On the Near Impossibility of Measuring the Returns to Advertising. Sus autores, quiero subrayarlo aquí, trabajan en Google y Microsoft.

Los métodos data driven gozan del mayor de los predicamentos. Véase una pequeña muestra extraída de una reciente conversación en Twitter:

Hoy he encontrado esto en Twitter:

Míralo bien. Vuelve a mirarlo. Efectivamente, los ricos votaron en contra de la independencia; los pobres, a favor. ¿Verdad?

Muchos, yo incluido, estamos inclinados a pensarlo así. Los resultados de una pequeña muestra que he hecho en la oficina han sido contundentes: todos, a pesar de sus doctorados, han estado de acuerdo unánimemente con el juicio anterior.

Así que ha sucedido lo siguiente:

Voy a escribir sobre un artículo como no debe hacerse: sin haberlo leído. Los bayesianos dirían que esta opinión que aquí voy a vertir es mi prior para cuando encuentre el tiempo y bajo la cual matizaré lo que en el se diga. Lo advierto, en todo caso, para que quien me lea no renuncie al sanísimo escepticismo.

Voy a hablar de Inferring causal impact using Bayesian structural time-series models y del paquete de R que lo acompaña, CausalImpact, cuyos autores trabajan en Google.

Tuve que saltarme una diapositiva en el DataBeers de Madrid del pasado jueves.

(A propósito, aquí están las 1+20 diapositivas.)

La decimonona, de la que trata la entrada, viene a hablar de lo siguiente. Tenemos una base de datos con sujetos (ids) que hacen cosas en determinados momentos. No es inhabitual calcular la frecuencia de esos sujetos así:

select id, count(*) as freq
from mytabla
where fecha between current_date - 7 and current_date
group by id
;

Esa variable se utiliza frecuentemente ya sea como descriptor de los sujetos o como alimento de otros modelos.

Una de las mayores contrariedades de estar sentado cerca de alguien que es más matemático que un servidor (de Vds., no de silicio) es que oye siempre preguntar por qué. Una letanía de preguntas me condujo a leer papelotes que ahora resumo.

Primero, unos datos:

set.seed(1234)

n <- 100

x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)

y <- 0.3 + 0.2 * x1 + 0.5 * (x2 > 0) + 0.2 * rnorm(n)

Luego, un modelo:

modelo <- lm(y ~ x1)
summary(modelo)

# Call:
#   lm(formula = y ~ x1)
#
# Residuals:
#   Min      1Q  Median      3Q     Max
# -0.9403 -0.2621  0.0420  0.2299  0.6877
#
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  0.55632    0.03364  16.538  < 2e-16 ***
#   x1           0.21876    0.03325   6.579 2.34e-09 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.3323 on 98 degrees of freedom
# Multiple R-squared:  0.3063,  Adjusted R-squared:  0.2992
# F-statistic: 43.28 on 1 and 98 DF,  p-value: 2.341e-09

Pocos que no entiendan cómo se han generado los datos advertirían lo malo de su especificación: hemos omitido una variable explicativa cuyo efecto ha ido a incrementar el error de manera que los tests habituales de bondad de ajuste no advierten.

Dos motivos me han tenido missing estas últimas semanas. Uno es una estancia en la Universidad de Santa Catalina del Burgo de Osma. Oportunamente ubicada en las estribaciones de la muy generosa en caldos de calidad Ribera del Duero, ha sido reconvertida a la sazón en un hotel propicio para la evasión y la agrafía.

El segundo es que en horas intempestivas he estado purgando de missings unas matrices enormes y de la, se conoce, mayor trascendencia. Es un asunto delicado, jamás bien resuelto, para el que el paquete [VIM](http://cran.r-project.org/web/packages/VIM/index.html) puede proporcionar ayuda. Sobre todo en los aspectos gráficos.

Muchos cuidamos de nuestro ordenador casi como una mascota: le ponemos un nombre (a menudo escribo desde tiramisu), le hacemos algo de mantenimiento, etc. Hay quienes, incluso, decoran sus máquinas con pegatinas.

Pero llega un momento en que hay que comenzar a tratar a las máquinas no tanto como mascotas sino como rebaños. Desde una pantalla aneja a esta en la que escribo estoy manejando un clúster de más de 200 GB y 50 núcleos distribuido en varias máquinas que ni sé dónde están. Además, solo espero que crezca. Ya no cuido de una mascota; cuido de un rebaño.

Partamos el tiempo en, p.e., días y contemos una serie de eventos que suceden en ellos. Es posible que esos recuentos se distribuyan según un proceso de Poisson de parámetro $latex \lambda$, que es un valor que regula la intensidad.

Si los días son homogéneos, i.e., no hay variaciones de intensidad diaria, estimar $latex \lambda$ (por máxima verosimilitud), es tan fácil como calcular la media de los sucesos por día. Pero puede suceder que la intensidad varíe en el tiempo (p.e., se reduzca los fines de semana). O que fluctúe de cualquier manera. O que haya periodos de gran intensidad y otros de calma. Es decir, que el proceso no sea homogéneo y que $latex \lambda$ varíe en el tiempo.