R

Un modelo jerárquico para lo de Casillas

Vuelvo a lo de Casillas inspirándome en el primer ejemplo de este artículo de Gelman et al.

El planteamiento es el siguiente: el número de paradas, $latex n_i$ que realiza el $latex i$-ésimo portero tiene una distribución binomial

$$ n_i \sim B(N_i, p_i)$$

donde $latex N_i$ es el número de disparos entre los palos y $latex p_i$ es la habilidad innata del portero. Estas habilidades innatas siguen una distribución dada, la de habilidades innatas de los porteros de primera división, que podemos suponer que sigue una distribución beta

Efectos en regresiones logísticas

Rescato y reconvierto un comentario de mi buen amigo José Luis Cañadas en una entrada mía reciente en la de hoy.

Sugiere José Luis el uso del paquete effects de R para estudiar el efecto de (que el caso concreto de interés, aunque hay otros) las variables de un modelo logístico.

Nos copia el código

library(effects)
mod.cowles <- glm(volunteer ~ sex + neuroticism*extraversion,
    data = Cowles, family = binomial)
eff.cowles <- allEffects(mod.cowles,
    xlevels = list(extraversion = seq(0, 24, 6)),
    given.values = c(sexmale = 0.5))
plot(eff.cowles, type = "response")

que genera

¿Son normales las alturas (de los individuos)?

Diríase que sí. La altura de un individuo está sujeta a multitud de factores que suman y restan. Está la genética (que es el resultado de la suma y resta del impacto de muchos genes individuales). Está la dieta, está… Diríase, insisto, que la altura es el promedio de muchos efectos pequeños y no demasiado dependientes entre ellos.

Y en efecto, (una vez descargados los microdatos de la Encuesta Nacional de Salud de 2011),

Una interpretación (rápida y sucia) de los coeficientes de la regresión logística

Los coeficientes de la regresión logística tienen una interpretación recta en términos de odds ratio. Que es un concepto sobre el que puede que alguien tenga algún tipo de intuición. Pero yo no.

¿Cómo podemos interpretar, aunque sea de manera rápida y grosera, los coeficientes? En términos de la variación de la probabilidad cuando la variable correspondiente cambia de valor (p.e., en una unidad). El problema es que la probabilidad depende del valor del resto de las variables: la relación no es lineal. No obstante, esa intuición es posible (en algunos casos: véase la nota).

R Consortium

R

Acaba de nacer el R Consortium con no sé qué objetivos. Los declarados son

trabajar con y dar soporte a la R Foundation y a las organizaciones clave que desarrollan, mantienen, distribuyen y usan R a través de la identificación, desarrollo e implementación de proyectos de infraestructura.

Vamos, nada que la R Foundation no viniese haciendo ya de oficio.

Solo que es probable que lo quieran hacer de otra manera. R se está volviendo demasiado importante, pensarán, como para que su desarrollo siga en manos de unos amateurs. Uso aquí el término en el sentido literal y no despectivo del término: estos amateurs son académicos de reconocido prestigio que hacen en su tiempo libre y con equipos informáticos que vete tú a saber a quién pertenecerán y dónde se alojarán las cosas que pretende hacer el consorcio.

Mejores mensajes de error con deparse + substitute

R
foo <- function(df, column.name){
    if (!column.name %in% colnames(df))
      stop("Column ", column.name, " not found in ", deparse(substitute(df)))

    mean(df$column.name)  # por ejemplo
  }

  foo(iris, "petal.area")

Lanza el error

Error in foo(iris, "petal.area") : Column petal.area not found in iris

mucho más informativo gracias a deparse + substitute.

Diferencia de medias a la bayesiana con salsa de stan

El habitual problema de la diferencia de medias suele formularse de la siguiente manera: hay observaciones $latex y_{1i}$ e $latex y_{2i}$ donde

$$ y_{ji} \sim N(\mu_j, \sigma)$$

e interesa saber si $latex \mu_1 = \mu_2$. Obviamente, se desconoce $latex \sigma$. De cómo resolvió Gosset el problema están los libros de estadística llenos. En R,

set.seed(1234)
N1 <- 50
N2 <- 50
mu1 <- 1
mu2 <- -0.5
sig1 <- 1
sig2 <- 1
 
y1 <- rnorm(N1, mu1, sig1)
y2 <- rnorm(N2, mu2, sig2)
 
t.test(y1, y2)
# Welch Two Sample t-test
#
# data:  y1 and y2
# t = 4.7059, df = 95.633, p-value = 8.522e-06
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#   0.5246427 1.2901923
# sample estimates:
#   mean of x  mean of y
# 0.5469470 -0.3604705

En rstan,

Busco viñetista (para MicroDatosEs)

R

Las viñetas son complementos importantes para un paquete, para que un usuario circunstancial pruebe y use un paquete. Uno de los míos, MicroDatosEs carece de ellas.

Me gustaría poder añadirle una o más que ilustraran cómo usarlo. Por ejemplo, para reproducir algunos de los números que ofrece el INE en sus notas de prensa.

Por eso te ofrezco la posibilidad de que te conviertas en viñetista. Eso te convertiría en colaborador del paquete (que es algo que cabe en un currículo). Eso sí, no subiría cualquier cosa.

SparkR 1.4: carga de ficheros CSV

R

He instalado Spark 1.4 recientemente y he comenzado a cacharrear. Antes de nada, quiero cargar datos.

Advierto que ha cambiado sustancialmente la API de SparkR. Entre otras novedades, desapareció (o más bien, se escondió) la función textFile, que permitía leer ficheros línea a línea. Ahora está pero no se exporta. La verás solo si haces SparkR:::textFile. ¿Signo de deprecación?

Se pueden crear un DataFrame (tablas distribuidas de Spark) a partir de un data.frame de R:

rPython & Anaconda

R

Nota: publico hoy en inglés en atención al público potencial de la entrada.

rPython lets R users call Python code. Anaconda is a completely free enterprise-ready Python distribution for large-scale data processing, predictive analytics, and scientific computing. Not surprisingly, some users want to call Anaconda Python rather than their system’s default Python.

However, Anaconda is a very particular package: unlike most other packages, whose files are scattered in a diversity of locations, it is self contained in a single directory. This helps Anaconda solve some problems, like the library hell. It is intended to provide the same experience regardless of the specifics of the host system.

La encuesta de presupuestos familiares, en MicroDatosEs

R

Hoy he subido una nueva versión del paquete MicroDatosEs a r-forge que incluye herramientas para cargar los datos de la Encuesta de Presupuestos Familiares.

Aún no está en CRAN, pero estáis invitados a probarla instalando la versión de desarrollo mediante

install.packages("MicroDatosEs", repos="http://R-Forge.R-project.org")

La parte del paquete que se encarga de la EPF es obra de Diego Paniagua, que es uno de los estudiantes del Experto en Data Science de la UTAD. La aportación a este paquete es, de hecho, parte de su proyecto final.

Liberado Spark 1.4

R

El anuncio de la liberación de la versión 1.4 de Spark se ha materializado. Está aquí.

¿Qué trae de novedad la versión 1.4? La integración con SparkR —antes había que instalarlo con algo de dolor independientemente— y, aparentemente, data.frames distribuidos y, cuentan, una sintaxis similar a la de dplyr —honestamente, hubiera preferido otra— para manipularlos.

Iré desgranando por aquí novedades. Y estoy pensando organizar una install & tutorial party un día de estos en Madrid. ¿Alguien se animaría?

Paralelismo en R: memo[rándum]

R

Esta es una nota que me dejo a mí mismo sobre paralelización en R para no tener que ir buscándola en otras partes:

library(parallel)

foo <- function(i){
  Sys.sleep(i)
}

cl <- makeCluster(4)

system.time(parSapply(cl, 1:4, foo))
# user  system elapsed
# 0.025   0.006   4.007

system.time(sapply(1:4, foo))
# user  system elapsed
# 0.039   0.033  10.001

stopCluster(cl)