Esta entrada es casi una referencia para mí. Cada vez tiro más de lme4 en mis modelos y en uno en concreto que tengo entre manos toca simular escenarios. Para lo cual, simulate.merMod.

Véamoslo en funcionamiento. Primero, datos (ANOVA-style) y el modelo que piden a gritos:

library(plyr)
library(lme4)

a <- c(0,0,0, -1, -1, 1, 1, -2, 2)
factors <- letters[1:length(a)]

datos <- ldply(1:100, function(i){
    data.frame(x = factors, y = a + rnorm(length(a)))
})
modelo <- lmer(y ~ (1 | x), data = datos)

El resumen del modelo está niquelado:

summary(modelo)

# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ (1 | x)
# Data: datos
#
# REML criterion at convergence: 2560.3
#
# Scaled residuals:
#     Min      1Q  Median      3Q     Max
# -3.6798 -0.6442 -0.0288  0.6446  3.3582
#
# Random effects:
#     Groups   Name        Variance Std.Dev.
# x        (Intercept) 1.5197   1.2328
# Residual             0.9582   0.9789
# Number of obs: 900, groups:  x, 9
#
# Fixed effects:
#     Estimate Std. Error t value
# (Intercept) -0.009334   0.412212  -0.023

En particular,

El primero es el chino. Es el que se aplicó a otras crisis víricas (SARS, etc.), a la viruela y a la polio. Consiste en aplicar medidas drásticas hasta que el virus desaparezca. De hecho, hay provincias en china que llegaron a tener un número importante de casos,

pero donde ya no quedan casos activos:

El otro es el escenario RU: el virus va a seguir entre nosotros y todos, en algún momento u otro vamos a pasar por él (o, más propiamente, a la inversa). En cuyo caso:

Desafortunadamente, el concepto de interacción, muy habitual en modelización estadística, no ha penetrado la literatura del llamado ML. Esencialmente, el concepto de interacción recoge el hecho de que un fenómeno puede tener un efecto distinto en subpoblaciones distintas que se identifican por un nivel en una variable categórica.

El modelo lineal clásico,

$$ y \sim x_1 + x_2 + \dots$$

no tiene en cuenta las interacciones (aunque extensiones suyas, sí, por supuesto).

[Este es un aviso para todos aquellos que depositan una excesiva fe en lo que nos cuenta el INE.]

La causa de muerte no es la causa de muerte. Al menos, necesariamente. Lo que el INE llama causa de muerte es una imagen distorsionada de la causa de muerte por culpa de un embudo administrativo.

Comiendo con unos epidemiólogos en el ISCIII hace un tiempo, me decían, con cierta envidia, cómo en otros países como Dinamarca, se registraban hasta ocho causas de muerte: la última, la concomitante, la… Y bromeaban diciendo que, al final, todos nos morimos de parada cardiorrespiratoria.

Este buscador de oro busca pepitas en su tramo de río. El río arrastra piedrecitas, muchas piedrecitas, y pepitas de oro, pocas pepitas de oro.

Tiene un artilugio que toma barro del río y que hace lo siguiente:

• Descarta casi todas las piedrecitas (y el resto las mete en una caja)
• Detecta casi todas las pepitas (y las mete en la misma caja)

Al final del día, ¿qué encontrará en la caja?

Existe una respuesta naive a la pregunta anterior que es la que tal vez en la que estéis pensando.

Pero, ¡ah!, existe también otra ciencia, la política, que no lo es menos que otras, en las que la evidencia dice cosas tales como:

Our results show that voters significantly reward disaster relief spending, holding the incumbent presidential party accountable for actions taken after a disaster. In contrast, voters show no response at all, on average, to preparedness spending, even though investing in preparedness produces a large social benefit.

[En esta entrada deambulo peligrosamente por los límites de un NDA; sin embargo, me siento obligado a exponerme a las posibles consecuencias debido a la gravedad de las circunstancias actuales.]

En España existe un mecanismo de monitorización de la mortalidad diaria por todas las causas. Su existencia no es explícitamente pública, pero sí que existen indicios implícitos de su existencia en informes de salud pública: véanse, p.e., referencias a MoMo y EuroMOMO aquí. [Nota: MoMo es el acrónimo de mortality monitoring].

[Retomando un tema que dejé inconcluso y que tampoco remataré hoy aquí.]

Imagina que quieres saber cuánto le dura a la gente el portátil. Para eso preguntas por ahí cuándo se compraron el último.

Lo que obtienes es un conjunto de datos donde todas las observaciones están censuradas. Y no, el análisis de la supervivencia clásico no funciona.

Buscando en la literatura he encontrado, sin embargo, Survival Analysis of Backward Recurrence Times, donde se discute el problema y al que, bueno, otro día con menos penas volveré.

El código usado en

Coronavirus: los nuevos casos diarios se estabilizan en muchos países menos en... pic.twitter.com/XOwxyccsZG

— Carlos Gil Bellosta (@gilbellosta) March 10, 2020

es

library(reshape2)
library(ggplot2)
library(plyr)

url <- "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv"
cvirus <- read.table(url, sep = ",", header = T)

cvirus$Lat <- cvirus$Long <- NULL
cvirus$Province.State <- NULL

cvirus <- melt(cvirus, id.vars = "Country.Region")

colnames(cvirus) <- c("país", "fecha", "casos")
cvirus$fecha <- as.Date(as.character(cvirus$fecha),
    format = "X%m.%d.%y")

tmp <- cvirus[cvirus$país %in% c("Italy", "Spain",
    "France", "Germany", "South Korea", "UK"),]

foo <- function(x){
    x <- x[order(x$fecha),]
    data.frame(fecha = x$fecha[-1],
        casos = diff(x$casos))
}

res <- ddply(tmp, .(país), foo)

res$país <- reorder(res$país, res$casos, function(x) -max(x))

res <- res[res$fecha > as.Date("2020-02-15"),]

ggplot(res, aes(x = fecha, y = casos)) +
    geom_point(size = 0.5) + geom_line(alpha = 0.3) +
    facet_wrap(~país, scales = "free_y") +
    ggtitle("Coronavirus: new daily cases") +
    theme_bw()

ggsave("/tmp/new_daily_cases.png", width = 12,
    height = 8, units = "cm")

Por referencia y afán de completar dos entradas que hice hace un tiempo sobre el método delta, esta y esta, dejo constar mención al paquete propagate, que contiene métodos para la propagación de la incertidumbre.

Para desavisados: si $latex x \sim N(5,1)$ e $latex y \sim N(10,1)$, ¿cómo sería la distribución de $latex x/y$? Etc.