Undécima sesión con Rcmdr: la distribución normal

Densidad de probabilidad y función de distribución (acumulada) de probabilidad

En las clases hemos visto aparecer la distribución normal de forma natural al considerar distribuciones binomiales B(n,p) para valores de n cada vez más grandes. Esta distribución continua ocupa el lugar más destacado entre todas las que estudiaremos en Estadística. Y en esta entrada vamos a aprender los comandos básicos para trabajar en R con la distribución normal.

La distribución normal es en realidad una familia de distribuciones que dependen de dos parámetros, \mu y \sigma, que se corresponden respectivamente con la media y desviación típica de una variable aleatoria X que siga esta distribución. Decimos entonces que X es de tipo N(\mu,\sigma).

Su función de densidad de probabilidad es

f_{\mu,\sigma}(x)=\dfrac{1}{\sigma\sqrt{2\pi}} e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}

En R podemos calcular los valores de esta función mediante el comando dnorm. Por ejemplo, para calcular el valor de f_{3,2}(4) (es decir, $\latex \mu=3, \sigma=2$ y X=4) usaríamos un código como este:

# Introducimos los parametros de la normal X de tipo  N(mu,sigma)
mu = 3
sigma = 2
# Calculamos la densidad en un valor x
X = 4
dnorm(X,mean=mu,sd=sigma)

y se obtiene como resultado 0.1760327. Antes de seguir, observemos que el comando dnorm sigue el esquema de nombres que ya anunciamos al hablar de la distribución binomial. Todas las funciones relacionadas con la distribución normal usan el sufijo norm, precedido de una letra que nos indica el tipo de resultado. Así, en esta entrada, además de dnorm, y por orden, vamos a conocer las funciones pnorm, qnorm y rnorm.

Aunque está bien saber calcular los valores de la función de densidad mediante dnorm, lo cierto es que, tratándose de una distribución continua, no vamos a usarlos prácticamente nunca. Pero, antes de despedirnos de dnorm,  vamos a aprovechar esta función para hacer una breve incursión en las capacidades gráficas de R. Para dibujar la función de densidad f, por ejemplo en el intervalo que va desde \mu-4\sigma hasta  \mu+4\sigma, usaremos el comando curve de esta manera

curve(dnorm(x,mean=mu,sd=sigma),from=mu-4*sigma,to=mu+4*sigma)

Este es el comando curve más básico, pero podemos añadir parámetros opcionales para controlar al detalle el aspecto final del gráfico. Proponemos aquí sólo un ejemplo para que puedas experimentar cambiando los valores de algunos parámetros.

curve(dnorm(x,mu,sigma),from=mu-4*sigma,to=mu+4*sigma,col="red",ylim=range(0, 0.5),lwd=6)

(Pistas: col viene de “color”, lwd de “line width” (ancho de línea) e ylim es “y limit”). El resultado es Como ya hemos dicho, aparte de los gráficos le vamos a dar poco uso a la función dnorm. En cambio, la función pnorm es muy interesante. Si tenemos una variable X de tipo N(\mu,\sigma) y queremos calcular el valor

P(X<k)

(es decir, lo que llamamos la cola izquierda de la distribución normal en el punto k) entonces podemos ejecutar estas instrucciones:

# Introducimos los parametros de la normal X de tipo  N(mu,sigma)
mu = 3
sigma = 2
# Calculamos la probabilidad de la cola izquierda de un valor k
k = 4
pnorm(k,mean=mu,sd=sigma)

Podemos calcular cualquier valor de probabilidad de la distribución normal combinando la función pnorm con las propiedades básicas de la probabilidad, y teniendo en cuenta que (tratándose de una distribución continua) no hay ninguna diferencia entre

P(X<k)    y     P(X\leq k)

Por ejemplo, si X es de tipo N(10,2) y queremos calcular la probabilidad

P(X>11)

usaríamos

1-pnorm(11,mean=10,sd=2)

Y si queremos calcular

P(7<X<12)

usaríamos

pnorm(12,mean=10,sd=2)-pnorm(7,mean=10,sd=2)

Estos ejemplos cubren, en realidad, todas las situaciones que se presentan en la práctica con problemas de probabilidad directos para la distribución normal. Recordemos que hemos llamado problemas directos a aquellos en los que la pregunta tiene la forma

P(dato)=??

Terminamos esta visita a la función pnorm con un encargo para el lector. Sean cuales sean \mu y \sigma, ¿cuál es el resultado de este comando?

pnorm(mu,mu,sigma)

La normal estándar

De entre todas las distribuciones normales, la más simple, y a la vez más destacada, es la que tiene media o y desviación típica 1, a la que llamamos normal estándar y que representamos con la letra Z. La normal estándar ocupa un lugar destacado en la Estadística, y en R no podía suceder otra cosa. De hecho, si no le damos ninguna información sobre la media y la desviación típica, R siempre asume que la normal de la que hablamos es la estándar. Así pues, si escribimos

pnorm(1.5)

entonces R asume que estamos interesados en conocer el valor de

P(Z<1.5)

donde, insistimos, Z representa a una normal estándar de tipo N(0,1) (y la respuesta es 0.9331928). En toda las funciones que vamos a ver en esta entrada se aplica el mismo principio: una normal, por defecto, es estándar para R, salvo que indiquemos su media y desviación típica explícitamente.

Problemas inversos. Cuantiles.

Un problema inverso de probabilidad (problema de cuantiles) es aquel en el que se la pregunta tiene la forma general

P(??)=dato

Por ejemplo,  en una distribución normal de tipo N(5,2), ¿cuál es el valor k para el que se cumple esta ecuación?

P(X< k)=1/3

Muy pronto veremos (en clase y en próximas entradas del blog) que este tipo de preguntas son extremadamente importantes cuando tenemos que calcular intervalos de confianza basados en la distribución normal (los primeros que encontraremos al adentrarnos en la inferencia).

En el caso de la distribución normal, y para responder a este tipo de preguntas, en R disponemos de la función qnorm. Para el ejemplo anterior la respuesta se obtiene ejecutando:

qnorm(1/3,mean=5,sd=2)

El  resultado es k=4.138545. Podemos verificarlo usando pnorm. Si se calcula

pnorm(4.138545,mean=5,sd=2)

Se obtiene como resultado 0.3333333, que es aproximadamente 1/3. Para otros problemas inversos de probabilidad tendremos que modificar un poco la llamada a la función. Por ejemplo, si en X, que sigue una distribución normal de tipo N(5,2), queremos saber cuál es el valor k para el que se cumple

P(X>k)=0.05

entonces debemos tener en cuenta que

P(X>k)=1-P(X<k)

y que por lo tanto estamos buscando el valor para el que se cumple

1-P(X<k)=0.05,

es decir

P(X<k)=1-0.05=0.95.

Eso nos lleva a usar el comando

qnorm(1-0.05,mean=5,sd=2)

para calcular el cuantil que deseábamos obtener. Ya hemos dicho que este tipo de cálculos jugarán un papel muy importante en la construcción de intervalos de confianza, así que volveremos a hablar de ellos cuando tratemos ese asunto en una próxima entrada del blog.

Finalmente, señalemos que la función qnorm sigue, como ya habíamos anunciado, el convenio de notación de R para las distribuciones de probabilidad: la q es por quantile, y norm identifica a la distribución normal.

Muestras aleatorias de distribuciones normales

Ya vimos como usar la función rbinom para generar valores aleatorios de una distribución binomial. Así que, siguiendo el convenio de notación, no debería sorprendernos que la función rnorm nos proporcione valores aleatorios de una distribución normal. Por ejemplo

rnorm(50,mean=100,sd=5)

nos proporciona 50 valores aleatorios de una distribución de tipo N(100,5). Obtendremos una salida similar a esta (distinta cada vez, claro, son muestras aleatorias):

>rnorm(50,mean=100,sd=5)
[1] 101.33299  98.01174 101.60873  94.73994 101.00582 103.88962 109.29085 103.30387 100.65213  95.60944 100.78252 106.23386 102.42284  99.72567  87.68830 101.69162 100.32440  99.91586
[19]  99.08630 101.33887  95.65277  96.94155  95.68294 103.15907 102.56853 101.74651  94.19673 101.46866  93.44603 103.69243 108.92536  95.41126 102.06478 109.04288  98.17461  98.70458
[37] 107.13320 105.30013 104.27446 101.50243  99.57164  97.45771  96.61668  96.50402  97.46888  92.25941 101.06347  97.29899  95.18880  88.48130

Gracias por la atención.

Décima sesión con Rcmdr: la distribución binomial.

Densidad de probabilidad y función de distribución (acumulada) de probabilidad

La distribución binomial es, sin duda, la más importante de las distribuciones discretas de probabilidad. Una variable X de tipo binomial B(n,p) (n es el número de ensayos y p la probabilidad de éxito en cada uno de esos ensayos) toma valores k de 0 a n, y se cumple que

P(X=k)=\binom{n}{k}p^kq^{n-k},

siendo q=1-p la probabilidad de fallo.

La media, varianza y desviación típica de una distribución binomials se obtienen fácilmente en R con estos comandos (aquí se muestra un ejemplo para una binomial de tipo B(25,0.5)):

# Introducimos los parametros de la binomial X de tipo  B(n,p)
n = 25
p = 0.5
q = 1-p
# Calculamos la media, varianza  y desviacion tipica
mu = n*p
mu
varianza= n*p*q
varianza
sigma = sqrt(n*p*(1-p))
sigma

Además, en R disponemos de varias funciones para calcular valores de probabilidad asociados a las distribuciones binomiales. Las vamos a explicar con detalle por dos razones. Primero, porque, como hemos dicho, la distribución binomial es muy importante. Y segundo,  porque una de las ventajas de R es que el manejo de todas las distribuciones que vamos a ver (y veremos unas cuantas) es homogéneo. Así que buena parte de lo que aprendamos aquí con la binomial nos será de utilidad cuando lleguemos a la normal, la de Poisson, la t de Student, etc.

La más fácil de entender es la función dbinom. Esta función sirve para calcular la probabilidad de que una variable binomial X de tipo B(n,p) tome un cierto valor k. Por ejemplo, con n=20, k=5 y p=3/7 usaríamos este comando:

dbinom(5,size=20,prob=3/7)

que produce como resultado 0.0506978.

Es importante desde el principio fijarse en los nombres que se usan en R para las funciones ligadas a distribuciones, porque el esquema de nombres será el mismo en todas las distribuciones que veamos. En la función dbinom, por ejemplo, la d indica la densidad de probabilidad. Más adelante, cuando veamos un nombre como dpois este convenio de nombres nos servirá para reconocer que es la densidad de probabilidad de una distribución tipo Poisson,

Siguiendo con la binomial, la siguiente función que vamos a  ver nos permite calcular su función de distribución (acumulada) de probabilidad F(k). Recordemos que la función de distribución de probabilidad  se define así:

F(k)=P(X\leq k)

Es decir, es la probabilidad de que X tome un valor igual o menor que k. Para calcular este valor en R se utiliza la función pbinom. Para entender la relación entre dbinom y pbinom vamos a usar el ejemplo anterior. Si ejecutamos este comando:

pbinom(5,size=20,prob=3/7)

estamos calculando lo mismo que si hiciéramos:

dbinom(0,size=20,prob=3/7)+dbinom(1,size=20,prob=3/7)+dbinom(2,size=20,prob=3/7)+dbinom(3,size=20,prob=3/7)+dbinom(4,size=20,prob=3/7)+dbinom(5,size=20,prob=3/7)

Y en ambos casos obtenemos 0.08014222 como resultado. De nuevo subrayamos que siempre que veamos una d delante del nombre de una distribución sabremos que se trata de su función de distribución acumulada (siempre con \leq).

Un comentario muy importante para evitar errores: cuando se trabaja con distribuciones discretas, como la binomial que nos ocupa, es esencial darse cuenta de que

P(X\leq 5)\mbox{ y } P(X<5)

son preguntas distintas.  Y la función pbinom siempre contesta a la pregunta con menor o igual (como veremos, esto es igual para todas las funciones de R que empiezan por p).  Si queremos obtener la respuesta a estas preguntas en R (con n=20, p=3/7, por ejemplo) debemos usar estos dos comandos:

#probabilidad P(X<=5)
pbinom(5,size=20,prob=3/7)
#probabilidad P(X<5)
pbinom(4,size=20,prob=3/7)

Así que, como se ve, para obtener P(X<k) usamos pbinom con k-1. Se debe tener mucho cuidado al responder a preguntas como P(X>=7) utilizando pbinom, para no confundirse al elegir los argumentos de la función.  En este fichero de comandos R tienes una recopilación de ejemplos de cálculo de probabilidades en muchas situaciones distintas.

Problemas directos y problemas inversos. Cuantiles.

Combinando las funciones dbinom y pbinom pueden resolverse todos los problemas directos de probabilidad asociados con la distribución binomial. ¿Qué quiere decir problemas directos? Nos referimos a los problemas de la forma

P(dato)=??

en los que la pregunta tiene la forma general “¿cuánto vale la probabilidad de que…?“. En cambio un problema inverso de probabilidad es aquel en el que se la pregunta tiene la forma general

P(??)=dato

Es decir, “¿cuál es el valor para el que la probabilidad vale….?”. Veamos un ejemplo de problema inverso: en una distribución binomial de tipo B(17,1/6), ¿cuál es el primer valor k (de 0 a 17) para el que se cumple P(X\leq k)>1/2? Para responder a este tipo de preguntas en R disponemos de la función qbinom. Concretamente, obtenemos la respuesta ejecutando:

qbinom(0.5,size=17,prob=1/6)

El  resultado es k=3. Podemos comprobarlo usando pbinom. Vamos a calcular varios valores de la función de distribución:

pbinom(0,size=17,prob=1/6)
pbinom(1,size=17,prob=1/6)
pbinom(2,size=17,prob=1/6)
pbinom(3,size=17,prob=1/6)

Y se obtiene como resultado:

> pbinom(0,size=17,prob=1/6)
[1] 0.04507324
> pbinom(1,size=17,prob=1/6)
[1] 0.1983223
> pbinom(2,size=17,prob=1/6)
[1] 0.4435207
> pbinom(3,size=17,prob=1/6)
[1] 0.6887192

Puedes ver que, en efecto, k=3 es el primer valor para el que la probabilidad acumulada supera el valor 1/2. Es muy posible que esta discusión te haya recordado la noción de mediana, o cuartiles y percentiles en Estadística Descriptiva. Allí se trataba de frecuencias, y aquí estamos hablando de probabilidades, pero la idea es similar: ¿cuál es el primer valor que deja a su izquierda el 50% (o el porcentaje que sea) de la probabilidad?  Por esa razón a estos valores que estamos calculando se les llama los cuantiles de la distribución. En inglés quantiles, y de ahí el nombre qbinom. Un cuantil es el valor que deja a su izquierda una cierta probabilidad. Y en R las funciones que empiezan por q siempre calculan quantiles y siempre a la izquierda (lo destacamos porque será muy importante recordarlo cuando tratemos con distribuciones continuas).

Un ejercicio muy saludable es tratar de obtener la respuesta a esta pregunta usando qbinom.  En una distribución de tipo B(17,1/6), ¿cuál es el menor valor k que cumple esta desigualdad?

P(X\geq k)=0.25

¡Cuidado, se trata de un mayor o igual!

Muestras aleatorias de distribuciones binomiales

En la anterior sesión vimos como usar la función sample de R para generar muestras aleatorias de los valores de una variable aleatoria discreta genérica, a partir de su densidad de probabilidad. En el caso de la distribución binomial R nos brinda una función, rbinom (la r viene de random values, valores aleatorios), para hacer esto mismo sin complicaciones. Si por ejemplo queremos generar 200 valores de una distribución binomial de tipo B(40,4/5) ejecutamos el comando:

rbinom(200,size=40,prob=4/5)

Y obtendremos una salida similar a esta (distinta cada vez, claro, son muestras aleatorias):

> rbinom(200,size=40,prob=4/5)
[1] 32 33 35 33 32 30 35 33 30 28 34 32 32 31 31 33 32 33 33 34 32 36 28 31 34 31 32 32 33 33 30 29 29 39 30 29 31 35 31 33 30 33 31 31 35 27 31 34 33 29 32 32 32 33 30 30 30 33 34 35 37 28 31 33 25 30 32 33 34 28 33 32 33 28 31 34
[77] 33 35 37 33 27 32 33 32 32 35 35 32 28 27 29 33 31 32 31 34 31 33 29 31 28 35 28 38 35 33 33 33 30 32 32 30 29 36 34 34 30 29 37 36 29 33 33 32 32 31 28 33 32 34 35 32 33 28 30 35 30 33 36 32 28 33 32 29 35 34 37 29 29 37 31 32
[153] 30 30 29 30 30 35 36 30 27 36 29 31 34 35 27 32 31 34 33 33 35 26 33 29 34 36 33 33 30 33 34 28 30 34 31 33 31 34 31 29 32 35 32 36 29 33 31 35

Gracias por la atención.

Novena sesión con Rcmdr: variables aleatorias discretas en R. Sampleando en R.

En esta entrada vamos a ocuparnos de los cálculos necesarios para trabajar con una
variable aleatoria discreta X, dada mediante una tabla de valores y probabilidades (la tabla de densidad de probabilidad) como esta:

\begin{array}{c|c|c|c|c|}valor&x_1&x_2&\cdots&x_k \\\hline probabilidad&p_1&p_2&\cdots&p_k\end{array}

En esta situación, queremos calcular \mu_X (la media de X) y \sigma_X (la desviación típica de X), y vamos a aprender a utilizar R para calcularlos. Para fijar ideas vamos a pensar en un ejemplo concreto. Supongamos que la densidad de probabilidad de la variable X es esta:

\begin{array}{c|c|c|c|c|c|}valor&2&4&7&8&11\\\hline probabilidad&1/5&1/10&1/10&2/5&1/5\end{array}

Entonces, para calcular la media en R usaríamos estos comandos:

#Tabla de densidad de la variable aleatoria discreta
valores=c(2,4,7,8,11)
probabilidades=c(1/5,1/10,1/10,2/5,1/5)

#Cálculo de la media
mu=sum(valores*probabilidades)
mu

#Cálculo de la varianza y la desviación típica
varianza=sum((valores-mu)^2*probabilidades)
varianza
sigma=sqrt(varianza)
sigma

Y se obtienen estos resultados:

> #Tabla de densidad de la variable aleatoria discreta
> valores=c(2,4,7,8,11)
> probabilidades=c(1/5,1/10,1/10,2/5,1/5)
> #Cálculo de la media
> mu=sum(valores*probabilidades)
> mu
[1] 6.9
> #Cálculo de la varianza y la desviación típica
> varianza=sum((valores-mu)^2*probabilidades)
> varianza
[1] 9.49
> sigma=sqrt(varianza)
> sigma
[1] 3.080584

También podemos representar gráficamente la distribución de probabilidad, mediante este comando:

barplot(probabilidades,names.arg=valores,col="lightseagreen")

El gráfico que se obtiene (en la ventana R, no en la de Rcmdr) esPara entender el comando barplot, empecemos por lo más sencillo: el último argumento col=”lightseagreen” es, simplemente, para dar color (verdoso, en lugar del gris por defecto) al gráfico.  SI quieres ver (una parte de) los colores disponibles en R, ejecuta el comando colors().  El primer argumento determina la altura de las columnas (que vienen dadas por las probabilidades, claro). Y el segundo argumento, names.arg=valores ,sirve para pedirle a R que rotule las columnas con los valores de la variable X, que están almacenados en ese vector valores.

Función de distribución (probabilidad acumulada)

La función de distribución de la variable aleatoria X es, recordémoslo:

F(k)=P(X\leq k)

Es decir, que dado el valor de k, debemos sumar todos los valores de la densidad de probabilidad para valores menores o iguales que k. En el ejemplo que estamos utilizando, si queremos calcular F(7), debemos hacer:

F(7)=P(X=2)+P(X=4)+P(X=7)=\frac{1}{5}+\frac{1}{10}+\frac{1}{10}

Se trata de sumas acumuladas. ¿Cómo hacemos esto en R? Usaríamos la función cumsum, que ya encontramos en la quinta sesión al hacer tablas de frecuencias acumuladas. Es decir,

#Función de distribución
cumsum(probabilidades)

que produce un vector con los valores de F(k) para cada $k$.

> cumsum(probabilidades)
[1] 0.2 0.3 0.4 0.8 1.0

Seguramente preferiríamos ver estos resultados en forma de tabla, para poder localizar el valor de F que corresponde a cada valor de k. Con lo que hemos aprendido de matrices en la anterior sesión, es fácil conseguirlo. Pondríamos:

#Función de distribución
rbind(valores,cumsum(probabilidades))

y obtenemos una tabla de valores de F:

> rbind(valores,cumsum(probabilidades))
[,1] [,2] [,3] [,4] [,5]
valores  2.0  4.0  7.0  8.0   11
0.2  0.3  0.4  0.8    1

Sampleando

Muchas veces queremos hacer experimentos o, mejor dicho, simulaciones, utilizando el ordenador para imitar lo que sucedería al obtener valores de una determinada variable aleatoria, de la que conocemos su densidad de probabilidad. Es decir, que queremos simular la obtención de una muestra de dicha variable. Como apuntamos en la sexta sesión, puede usarse la función sample (sample es muestra en inglés) de R para hacer esto. En este ejemplo, para obtener una muestra de 100 valores de X haríamos:

muestra=sample(valores,1000,prob=probabilidades,replace=T)

Como se ve, debemos explicarle a sample cuáles son los valores, y cuáles son sus probabilidades (mediante prob=). Además, en el segundo argumento indicamos cuántos elementos queremos (1000). Finalmente, para que sea posible obtener 1000 valores de ua variable que sólo toma cinco valores distintos, debemos permitir que sample repita valores. Por eso aparece el argumento replace=T. Mediante ese argumento R repetirá los valores, con frecuencias que se corresponden con la densidad de probabilidad de X.

Conclusión

Para facilitarte el trabajo con las variables aleatorias discretas, aquí tienes un resumen de todos los comandos que hemos visto en esta entrada. Puedes usar este fichero como plantilla cuando tengas que trabajar con este tipo de problemas, y sólo tendrás que cambiar los vectores de valores y probabilidades.

Y muchas gracias por la atención.

#Tabla de densidad de la variable aleatoria discreta
valores=c(2,4,7,8,11)
probabilidades=c(1/5,1/10,1/10,2/5,1/5)

#Cálculo de la media
mu=sum(valores*probabilidades)
mu

#Cálculo de la varianza y la desviación típica
varianza=sum((valores-mu)^2*probabilidades)
varianza
sigma=sqrt(varianza)
sigma

#Grafico de columnas de la distribucion de probabilidad
barplot(probabilidades,names.arg=valores,col="lightseagreen")

#Función de distribución
rbind(valores,cumsum(probabilidades))

#Muestra. Elegirprimero el numero de elementos de la muestra (y si hay remplazamiento).
numeroElementosMuestra=10000
muestra=sample(valores,numeroElementosMuestra,prob=probabilidades,replace=T)
muestra
mediaMuestral=mean(muestra)
mediaMuestral

Octava sesión con Rcmdr. Datos cualitativos. Matrices y dataframes.

Manejando variables cualitativas en R.

En la anterior sesión con Rcmdr aprendimos a utilizar (de forma básica) las matrices de R. Para refrescar la memoria puedes empezar ejecutando este código, en el que la única novedad es la última fila.

library("gtools") #recuerda que has debido instalar este paquete (ver septima sesion del blog)
combinaciones=combinations(4,3)
combinaciones
# de que tipo es el objeto combinaciones que hemos creado?
class(combinaciones)

Cuando lo ejecutes obtendrás como última linea de la salida:

> class(combinaciones)
[1] "matrix"

Es decir, que el objeto combinaciones es de tipo matrix. Puedes utilizar la función class de R para averiguar el tipo de cualquiera de los objetos que aparecen en una sesión de trabajo con R. Por ejemplo, prueba a usarla con un vector de números:

datos=c(2,5,17,21,9,7)
class(datos)

y observa el resultado. Verás que el vector datos es de tipo numeric.

Todo nuestro trabajo, en las anteriores sesiones del blog se ha centrado en números, en valores de variables cuantitativas. Pero sabemos que a veces es necesario trabajar con variables cualitativas.Recuerda que una variable cualitativa se usa para establecer clasificaciones nominales en los datos. Por ejemplo, la clasificación en hombre o mujer entre los pacientes que siguen un cierto tratamiento. O la clasificación taxonómica por especies.

Es cierto que siempre podríamos codificarlas mediante números.  Pero no es menos cierto que lo más cómodo es poder utilizar nombres como valores de las variables. Para hacer esto posible, R nos permite crear un cierto tipo de valores, que  son simplemente palabras o frases entrecomilladas (en general las llamamos cadenas de texto).  Por ejemplo, este vector podría representar el género de los diez pacientes que se han sometido a un cierto tratamiento:

pacientesPorGenero=c("mujer","mujer","hombre","mujer","hombre","hombre","hombre","mujer","mujer","mujer")
class(pacientesPorGenero)

Cuando ejecutes este código,obtendrás como último resultado

> class(pacientesPorGenero)
[1] "character"

Es decir, que los elementos del vector pacientesPorGenero son, para R, de tipo character. Ese es el tipo de dato que se utiliza en R para representar los valores de una variable cualitativa. Naturalmente, sin intentas  realizar operaciones numéricas con ese vector, como estas,

pacientesPorGenero^2
mean(pacientesPorGenero)

R te obsequiará con un surtido de insultos más o menos ofensivos en la ventana de mensajes de Rcmdr. En el segundo caso (el intento fallido de calcular la media), el valor que devuelve R es interesante: se obtiene NA, que es el valor que R devuelve cuando un resultado numérico es imposible de calcular, o no está disponible (Not Available, de ahí el nombre).

Una última observación: los vectores de tipo character se pueden ordenar, alfabéticamente, pero es muy posible que esa ordenación no tenga sentido en nuestro problema. Volveremos sobre esto en otra entrada, cuando tratemos el tema de los factores en R, que nos permitirán una gestión mucho más eficiente de la información que contienen las variables cualitativas.

Matrices a partir de un vector

En la sétima sesión obtuvimos varias matrices como resultado de las operaciones combinatorias que hacíamos. Pero ¿y si yo quiero rellenar una matriz con mis datos, escribiéndolos directamente?

Podemos hacer esto usando la función matrix de R. Como argumentos de esta función escribiremos los elementos de la matriz como un vector, y el número de filas mediante la opción nrow. Por ejemplo:

unaMatriz=matrix(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15),nrow=3)
unaMatriz
class(unaMatriz)

producirá este (seguramente sorprendente) resultado:

> unaMatriz=matrix(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15),nrow=3)
> unaMatriz
[,1] [,2] [,3] [,4] [,5]
[1,]    1    4    7   10   13
[2,]    2    5    8   11   14
[3,]    3    6    9   12   15
> class(unaMatriz)
[1] "matrix"

Obsérvese que R ha ido colocando los elementos por columnas, hasta formar una matriz de tres filas y cinco columnas (de acuerdo con el argumento nrow=3 que hemos usado).

Naturalmente también podríamos haber fabricado el vector usando cualquiera de los métodos que vimos en la sexta sesión, escribiendo por ejemplo:

unaMatriz=matrix(1:15,nrow=3)

y el resultado sería el mismo. Asimismo podemos emplear un vector guardado previamente mediante su nombre (y de esa manera puedes, por ejemplo, con lo que ya sabes de anteriores sesiones, fabricar una matriz a partir de los datos de un fichero csv).

¿.Y si quiero rellenar por filas (porque me resulta más natural, probablemente)? Pues basta con añadir una opción byrow=T a la  función matrix:

unaMatriz=matrix(1:15,nrow=3,byrow=T)
unaMatriz

Y ahora se obtiene

> unaMatriz=matrix(1:15,nrow=3,byrow=T)
> unaMatriz
[,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15

como queríamos.

Otra opción es usar la función t (simplemente la letra t), que traspone una matriz. Prueba a ejecutar este código para ver como funciona:

unaMatriz=matrix(1:15,nrow=3)
unaMatriz
matrizTranspuesta=t(unaMatriz)
matrizTranspuesta

Una pregunta natural cuando se usa la función matrix es ¿y qué pasa si las dimensiones de la matriz no casan con la longitud del vector? Por ejemplo, si tratamos de meter 13 datos en una matriz de  3 filas:

unaMatriz=matrix(1:13,nrow=3)
unaMatriz

El resultado es que se obtiene:

> unaMatriz
[,1] [,2] [,3] [,4] [,5]
[1,]    1    4    7   10   13
[2,]    2    5    8   11    1
[3,]    3    6    9   12    2

y una advertencia en la ventana de mensajes. R nos avisa de que el vector no era suficiente para rellenar la matriz y por esa razón ha vuelto al principio del vector, utilizando 1 y 2 para rellenar la parte restante de la matriz. Esta forma de trabajar de R nos puede resultar más o menos conveniente o convincente, pero en cualquier caso hay que conocerla.

Matrices por filas y columnas: cbind y rbind.

Otra manera de fabricar matrices a partir de vectores es utilizar las funciones rbind y cbind (donde r y c proceden de row, fila y column, columna, respectivamente) para escribir las filas o columnas directamente como vectores. Por ejemplo, el código

vectorFila1=2:5
vectorFila2=1:4
vectorFila3=c(-2,1,4,7)
matrizPorFilas=rbind(vectorFila1,vectorFila2,vectorFila3)
matrizPorFilas

produce este resultado:

> matrizPorFilas
[,1] [,2] [,3] [,4]
vectorFila1    2    3    4    5
vectorFila2    1    2    3    4
vectorFila3   -2    1    4    7

Y para ilustrar la función cbind vamos a añadirle una columna más a esta matriz, escribiendo:

quintaColumna=c(11,17,-4)
matrizModificada=cbind(matrizPorFilas,quintaColumna)
matrizModificada

que produce:

> matrizModificada
                          quintaColumna
vectorFila1  2 3 4 5            11
vectorFila2  1 2 3 4            17
vectorFila3 -2 1 4 7            -4

Obsérvese que cuando construimos la matriz mediante cbind y rbind, R usa los nombres de los vectores como nombres de las filas y columnas correspondientes. Volveremos pronto sobre esto (en una próxima entrada).

Dataframes

Después de conocer la forma en que R maneja las variables cualitativas, y de aprender un poco más sobre matrices, podemos estar tentados de combinar ambas ideas. Por ejemplo, puede resultar útil disponer de una tabla en la que guardamos, para cada uno de los pacientes ingresados en una planta de un hospital, su género (hombre o mujer),  su edad y el número de días que llevan ingresados. Algo como el contenido de este fichero pacientes.csv .

Puesto que cada columna contiene unos datos, podemos introducir los datos mediante vectores por columnas (usando cbind). El código sería:

vectorColumna1=c("hombre","mujer","mujer","mujer","hombre","hombre","mujer","mujer","hombre")
vectorColumna2=c(53,67,31,56,25,43,41,76,17)
vectorColumna3=c(3,1,4,2,1,7,2,4,8)
matrizPorColumnas=cbind(vectorColumna1,vectorColumna2,vectorColumna3)
matrizPorColumnas

y el resultado es:

> matrizPorColumnas
vectorColumna1 vectorColumna2 vectorColumna3
[1,] "hombre"       "53"           "3"
[2,] "mujer"        "67"           "1"
[3,] "mujer"        "31"           "4"
[4,] "mujer"        "56"           "2"
[5,] "hombre"       "25"           "1"
[6,] "hombre"       "43"           "7"
[7,] "mujer"        "41"           "2"
[8,] "mujer"        "76"           "4"
[9,] "hombre"       "17"           "8"

¿Perfecto, verdad? Espera un momento, algo va mal; esas comillas en las columnas de números no deberían estar ahí…

Para comprobar que en efecto tenemos un problema, tratemos de calcular la media de esa columna. Primero vemos el contenido de la columna (con matrizPorColumnas[,2]):

matrizPorColumnas[,2]
mean(matrizPorColumnas[,2])

Y al tratar de calcular la media obtenemos NA y un mensaje de advertencia en la ventana de mensajes.:

[18] NOTA: Aviso en mean.default(matrizPorColumnas[, 2]) :
argument is not numeric or logical: returning NA

¿Qué ha sucedido? Pues que nos hemos topado con la principal limitación de las matrices de R (que afecta también a los vectores, aunque hasta ahora no se nos había planteado el problema). Una matriz está formada por elementos que son todos del mismo tipo. Todos numeric, o todos de tipo character, pero no se admiten mezclas. Por eso, cuando le hemos pedido a R que usara los tres vectores columna para hacer una matriz, ha hecho lo que ha podido: ha convertido las columnas de números a tipo character, y las ha usado como tales (la conversión contraria, convertir character a numeric casi nunca tiene sentido, por eso R ha hecho lo que ha hecho).

Vaya decepción. Entonces, ¿en R no puedo escribir una tabla en la que se mezclen números y nombres? ¡Los dataframes al rescate! Un dataframe está pensado exactamente para eso, para albergar una tabla con tipos mixtos de datos. Veamos como funcionaría en este ejemplo:

pacientes=data.frame(vectorColumna1,vectorColumna2,vectorColumna3)
pacientes

Como puede verse, ahora no hay problema, el resultado es:

> pacientes
vectorColumna1 vectorColumna2 vectorColumna3
1         hombre             53              3
2          mujer             67              1
3          mujer             31              4
4          mujer             56              2
5         hombre             25              1
6         hombre             43              7
7          mujer             41              2
8          mujer             76              4
9         hombre             17              8

y si queremos calcular la media de la segunda columna, podemos hacer:

mean(pacientes[,2])

y se obtiene el resultado deseado.

> mean(pacientes[,2])
[1] 45.44444

Puesto que la columna dos tiene nombre (vectorColumna2), podemos conseguir lo mismo con esta otra notación  con el símbolo $ (volveremos a encontrarla en futuras entradas):

mean(pacientes$vectorColumna2)

Y ya que estamos en ello, es el momento de darse cuenta de que los nombres de las columnas en nuestro flamante data.frame no son especialmente afortunados. Podríamos haber sido más cuidadosos al elegir los nombres de los vectores, pero incluso si no lo hemos sido, todavía tiene remedio. Podemos hacer:

colnames(pacientes)=c("genero","edad","diasHospital")
pacientes

Y ahora se obtiene:

> pacientes
genero edad diasHospital
1 hombre   53            3
2  mujer   67            1
3  mujer   31            4
4  mujer   56            2
5 hombre   25            1
6 hombre   43            7
7  mujer   41            2
8  mujer   76            4
9 hombre   17            8

como queríamos.

En estos primeros ejemplos ya hemos visto que los dataframes se pueden manejar de manera muy parecida a como se hace con las matrices. En próximas entradas iremos aprendiendo más sobre ellos, por el procedimiento de usarlos. De momento, si quieres ver más ejemplos, te conviene saber que R incorpora una buena cantidad de dataframes de ejemplo. Para ver la lista escribe el comando

data()

y aparecerá (en otra ventana). Uno de los ejemplos es el dataframe llamado crabs. Puedes ver su contenido tecleando en R simplemente

crabs

y verás aparecer una tabla de datos (de 200 filas), con datos morfológicos de cangrejos (crabs, en inglés) del género Leptograpsus. Naturalmente, sin más información se hace difícil interpretar esta tabla. Para averiguar más sobre el fichero (o sobre casi cualquier otra cosa en R) puedes escribir su nombre precedido de una interrogación y ejecutar esa línea:

?crabs

Cuando lo hagas se abrirá una ventana de tu navegador (no te preocupes, aunque se abra el navegador no se necesita estar conectado a internet para esto), con información sobre el fichero, como en esta figura:

Finalmente, dataframes y ficheros csv. Funciones write.table y read.table

De la misma forma que hemos aprendido a usar scan y cat para leer y escribir vectores en ficheros csv, podemos leer y escribir dataframes a esos ficheros. Por ejemplo, para escribir el dataframe crabs a un fichero crabs.csv (que puedes abrir en Calc, por ejemplo), usaríamos la función write.table así (acuérdate de fijar el directorio de trabajo primero):

write.table(crabs,"crabs.csv",sep=" ",dec=",")

Un par de comentarios para explicar las opciones:

  1. la opción dec=”,” sirve para pedirle a R que al escribir el fichero reemplace los puntos decimales que emplea por comas. Esto es necesario si vas a usar el fichero en una versión española de Calc, que usa comas para los decimales.
  2. la opción sep=” ” le dice  a R que use espacios para separar los valores dentro de una línea del fichero csv. En este caso es necesario usar espacios, y no comas,  porque ya hemos usado comas para los decimales.

Si lo que queremos es leer un fichero de datos como el fichero pacientes.csv que vimos antes en esta misma entrada, podemos utilizar la función read.table así (¡recuerda fijar el directorio de trabajo!)

datosPacientes=read.table("pacientes.csv")
datosPacientes

Estos comandos producen esta salida:

> datosPacientes=read.table("pacientes.csv")
> datosPacientes
V1    V2
1 hombre ,53,3
2  mujer ,67,1
3  mujer ,31,4
4  mujer ,56,2
5 hombre ,25,1
6 hombre ,43,7
7  mujer ,41,2
8  mujer ,76,4
9 hombre ,17,8

donde vemos que, en efecto, R ha creado un nuevo data.frame datosPacientes, y le ha asignado automáticamente nombres V1 y V2 a las variables columna. En futuros encuentros con la función read.table aprenderemos más detalles sobre su funcionamiento.

Muchas gracias por la atención, especialmente en esta entrada tan extensa.

Séptima sesión con Rcmdr: paquetes y combinatoria. Datos en tablas

En la anterior sesión dejamos pendiente el problema de simular los lanzamientos de un par de dados. Queremos obtener resultados como (1,3) , (5,2), (4,3), etc. Vamos a utilizar este problema para aprender algunas cosas más sobre R.

Instalando paquetes adicionales

Más adelante, cuando veamos las estructuras básicas de bucles e iteración en R, dispondremos de otras maneras de construir los resultados. Pero, por el momento, para generar los valores que queremos, es mejor que nos aprovechemos del trabajo que otros han hecho por nosotros.

Una de las ventajas que se desprenden del hecho de que R sea un programa de código abierto esla gran cantidad de paquetes de código R que se han escrito, para resolver los problemas más variados. Ya hemos visto cómo guardar nuestros propios fragmentos de código en R. Otros usuarios, a menudo expertos en R y en alguna de sus aplicaciones, escriben sus propios fragmentos de código, para resolver problemas (a veces muy sofisticados). Y después ponen esos fragmentos de código, llamados a veces paquetes o librerías,  a disposición de la comunidad de usuarios de R, de forma libre y gratuita. A fecha de hoy, R-cran, que es el mayor almacén de código R,  contiene más de 4000 de esos paquetes. Para que te hagas una idea de lo que aportan, y lo elaborados que pueden lleagr a ser, la interfaz R-Commander que estamos usando es uno de esos paquetes.

Cuando instalas  R en tu ordenador, normalmente sólo instalas una pequeña parte de esos paquetes, para no hacer la instalación innecesariamente larga y complicada.Si quieres ver la lista de paquetes instalados en tu máquina, en la ventana Rgui (no es Rcmdr) utiliza el menú PaquetesCargar paquete…, y verás esa lista. Si en algún momento necesitas uno de los paquetes  que no llegaste a instalar, puedes conectar con un almacén en Internet, descargarlo e instalarlo de forma casi automática. Veamos cómo.

Para la sesión de hoy quiero usar un paquete que se llama gtools, que contiene, entre otras cosas, código para hacer Combinatoria. Cuando R arranca, lo normal es que sólo algunos de los paquetes instalados en nuestra máquina están disponibles (de nuevo para ahorrar recursos).  Para saber si el paquete gtools es uno de los que están instalados puedo buscarlo en la lista que hemos visto antes, o puedo teclear este comando, con el que le digo a R que quiero usar ese paquete:

library("gtools")

En este caso, el color rojo de la ventana de mensajes nos indica que algo ha ido mal, porque el paquete gtools no está instalado en nuestra máquina. Para instalarlo (conectados a internet, claro) tenemos que dirigirnos a la ventana  titulada Rgui, y usar el menú Paquetes → Instalar paquete(s)… Se abrirá una ventana, titulada CRAN mirror, en la que elegimos un mirror o almacén de paquetes R  Es conveniente elegir uno cercano geográficamente, así que yo voy a elegir Spain (Madrid). Tras unos instantes, en los que R descarga la lista de todos los paquetes disponibles en ese almacén, aparece una larga lista alfabética de paquetes. Ya casi hemos acabado, basta con localizar en esta lista el paquete gtools y pulsar en Ok

Si es el primer paquete que instalamos aparecerá una pregunta, en la que nos piden permiso para crear un “almacén local” de paquetes descargados en nuestro ordenador. Aceptamos, y comienza la descarga e instalación del paquete.

Cuando la instalación acabe, vuelve a Rcmdr, limpia la ventana de mensajes, y ejecuta de nuevo el comando library(“gtools”).  La ventana de mensajes ahora debería aparecer en verde. Ya podemos empezar el trabajo.

Permutaciones y combinaciones.

Los problemas típicos de la  Combinatoria más sencilla consisten en formar las combinaciones o permutaciones de unos objetos dados. En España distinguimos también variaciones, pero la tradición anglosajona considera las variaciones como otro tipo de permutaciones. Recordemos que la diferencia entre unas y otras estriba en que, en las combinaciones, el orden en que se colocan los objetos no importa. En cambio en las permutaciones-variaciones, dos listas se consideran distintas si los objetos aparecen en ellas en un orden distinto.

El paquete gtools añade dos funciones nuevas a R, permutations y combinations, que nos van a servir precisamente para obtener esos resultados.  No se trata de saber cuántos hay, sino, de hecho, de construirlos y enumerarlos. Para verlas en acción vamos a utilizar estos comandos para fabricar todas las listas posibles de tres elementos, elegidos de entre cuatro posibles. Para empezar, no admitimos repeticiones de los elementos. Entonces, si el orden no importa, estamos formando las combinaciones de cuatro elementos tomados de tres en tres. Y si el orden importa, entonces se trata de las variaciones (un inglés diría permutaciones) de cuatro elementos tomados de tres en tres.

Los correspondientes comandos son (recuerda que has debido utilizar primero library(“gtools”)):

combinations(4,3)
permutations(4,3)

Y la salida que se obtiene es, para el primer comando (combinaciones):

> combinations(4,3)
[,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    2    4
[3,]    1    3    4
[4,]    2    3    4

mientras que para el segundo comando (permutaciones), la salida es

> permutations(4,3)
[,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    2    4
[3,]    1    3    2
[4,]    1    3    4
[5,]    1    4    2
[6,]    1    4    3
[7,]    2    1    3
[8,]    2    1    4
[9,]    2    3    1
[10,]    2    3    4
[11,]    2    4    1
[12,]    2    4    3
[13,]    3    1    2
[14,]    3    1    4
[15,]    3    2    1
[16,]    3    2    4
[17,]    3    4    1
[18,]    3    4    2
[19,]    4    1    2
[20,]    4    1    3
[21,]    4    2    1
[22,]    4    2    3
[23,]    4    3    1
[24,]    4    3    2

que, evidentemente es más larga, puesto que ahora tenemos en cuenta la ordenación.

¿Y para lanzar dos dados? Bueno, si queremos obtener resultados equiprobables, debemos distinguir el orden (como si cada dado fuera de un color), y además, a diferencia de los dos ejemplos anteriores, debemos permitir resultados repetidos. Con esas premisas, hacemos

dosDados=permutations(6,2,repeats.allowed=TRUE)
dosDados

y obtenemos  los 36 resultados esperados:

[,1] [,2]
[1,]    1    1
[2,]    1    2
[3,]    1    3
[4,]    1    4
[5,]    1    5
[6,]    1    6
[7,]    2    1
[8,]    2    2
[9,]    2    3
[10,]    2    4
[11,]    2    5
[12,]    2    6
[13,]    3    1
[14,]    3    2
[15,]    3    3
[16,]    3    4
[17,]    3    5
[18,]    3    6
[19,]    4    1
[20,]    4    2
[21,]    4    3
[22,]    4    4
[23,]    4    5
[24,]    4    6
[25,]    5    1
[26,]    5    2
[27,]    5    3
[28,]    5    4
[29,]    5    5
[30,]    5    6
[31,]    6    1
[32,]    6    2
[33,]    6    3
[34,]    6    4
[35,]    6    5
[36,]    6    6

Ahora que ya vemos que esto funciona, podemos pararnos un momento y preguntarnos: estas respuestas que estamos obteniendo, que evidentemente no son vectores, ¿qué són? Y la respuesta es que los objetos como este dosDados que acabamos de fabricar son matrces, una estructura adecuada para almacenar datos en tablas.

Tablas de datos en R, un primer contacto

Las matrices son una de las formas que utiliza R para almacenar y organizar la información. (A este tipo de objetos -vectores, matrices, tablas, etc.- se los conoce en Informática como Estructuras de Datos). Su manejo tiene muchos parecido con lo que ya hemos aprendido en el caso de los vectores. Si vuelves a mirar la anterior matriz (la de los dos dados), verás que R la ha adornado con corchetes, que son de hecho una forma de indicarnos cómo se accede a los distintos elementos de la matriz. Por ejemplo, si escribes:

dosDados[23,]

obtendrás  como respuesta

> dosDados[23,]
[1] 4 5

Es decir, la fila número 23 de la matriz. Si lo que quieres es obtener su segunda columna, escribe

dosDados[,2]

y la respuesta de R será un vector con los elementos de esa columna:

> dosDados[,2]
[1] 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6

Como puede verse, los números entre corchetes que rodean la matriz nos ayudan a recordar  todo esto.  Si queremos, por poner un ejemplo, calcular la suma de los dados en cada tirada, podemos hacer:

sumaTiradas=dosDados[,1]+dosDados[,2]

y veremos aparecer un vector con esas sumas. Y si necesitas saber las dimensiones de la matriz (cuántas filas y cuuántas columnas tiene), puedes usar la función dim:

dim(dosDados[,1])

para averiguar esas dimensiones.Por su parte, la función length te dirá cuántos elementos tiene la matriz (en este caso, 72).

También podemos usar la notación de corchetes para seleccionar los elementos de la matriz que cumplen alguna condición, como hacíamos con los vectores. Por ejemplo, las filas en las que el segundo dado (la segunda columna) es 4 o más, se seleccionan con:

dosDados[dosDados[,2]>=4,]

que produce como resultado una nueva matriz con los resultados que queremos. Atención: los números de fila de esta nueva matriz no se corresponden con los de la original. Antes de cerrar esta entrada, te propongo un pequeño ejercicio: trata de adivinar lo que se obtiene con este comando:

dosDados[(dosDados[,1]+dosDados[,2])>7,]

En esta entrada hemos avanzado bastante, y hemos conseguido nuestro objetivo de representar las 36 tiradas diferentes de dos dados. Todavía podemos ir un poco más allá, pero para eso necesitamos aprender sobre estructuras de control. Con eso estaremos en condiciones de escribir auténticos programas en R. Eso lo veremos en las próximas entradas.

Gracias por la atención.

Sexta sesión con Rcmdr: más sobre vectores y combinatoria.

Cualquier curso de introducción a la Estadística incluye, a su vez, una introducción a la Probabilidad. Y los ejemplos iniciales de cálculo de probabilidades recurren necesariamente (por simplicidad, y por razones históricas) a juegos de azar, como los dados, lanzamientos de monedas, loterías, ruletas, etc. Estas situaciones simplificadas, e idealizadas, sirven de modelo a muchas otras que encontramos después en aplicaciones “más serias”. Y en muchas de ellas juega un papel la Regla de Laplace, que en esencia dice que, en un experimento en el que hay una familia de sucesos elementales equiprobables (SEE), la probabilidad del suceso A es:

P(A)=\dfrac{\mbox{n\'umero de SEE en los que ocurre A}}{\mbox{n\'umero total de SEE posibles}}

Para utilizar esta regla efizcamente hay que contar los sucesos que cumplen una cierta condición. Y contar, por elemental que pueda parecer, es una actividad extremadamente complicada. La rama de las matemáticas que se encarga de este tipo de problemas es la Combinatoria. Y en esta entrada vamos a empezar a utilizar la combinatoria como excusa para aprender un poco más sobre cómo se fabrican vectores a medida con R, y a mejorar nuestras habilidades con esos vectores.

Fabricando vectores a medida

Ya sabemos que, en R, podemos crear un vector con la función c, simplemente  enumerando sus elementos, como en

unVector=c(2,7,2,3,1,5)

Esto está muy bien, pero si queremos fabricar así un vector que contenga los números del 1 al 1000, es muy probable que antes de acabar nos invada la melancolía (por no mencionar el hecho de que, casi con seguridad, cometeremos algún error). Afortunadamente, R dispone de mecanismos para simplificar este tipo de tareas repetitivas. En este ejemplo, basta con usar:

milNumeros=1:1000
milNumeros

La segunda líneas, ya lo sabemos, es para ver el resultado. Y se obtiene:

(Por cierto, el resultado puede servir para entender porque R no muestra ningún resultado a menos que se lo pidamos explícitamente).

Es un avance sin duda, pero aún no resuelve muchos problemas. ¿Y si queremos los números del 1 al 1000, pero saltando de tres en tres? El lector mañoso tal vez habrá pensado “hago un vector con los primeros 333, y lo multiplico por 3…”. Pero en R hay una manera mucho más sencilla, usando la función seq. Hacemos:

deTresEnTres=seq(1,1000,3)
deTresEnTres

Y ahora la salida es:exactamente lo que queríamos.

La función seq permite fabricar vectores cuyos elementos se diferencian en una cantidad fija; que puede ser cualquier número, entero, fraccionario, expresado con decimales, etc.  Hay otra función, la función rep, estrechamente emparentada con est, cuyo principal uso es fabricar vectores siguiendo un patrón. El nombre rep viene de repeat, repetir. Y eso ayuda a adivinar lo que sucede si escribimos estos comandos:

datos=c(1,2,3,4,5)
rep(datos,3)

La salida correspondiente es, como cabría esperar, el vector datos repetido tres veces:

> datos=c(1,2,3,4,5)
> rep(datos,3)
[1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

Quizá no sea tan fácil adivinar lo que sucede con este otro comando:

datos=c(1,2,3,4,5)
rep(datos,c(2,9,1,4,3))

Si lo ejecutamos,  la salida es:

> datos=c(1,2,3,4,5)
> rep(datos,c(2,9,1,4,3))
[1] 1 1 2 2 2 2 2 2 2 2 2 3 4 4 4 4 5 5 5

Y ahora está más claro lo que sucede: el segundo argumento nos dice cuántas veces debemos repetir cada uno de los correspondientes elementos del vector datos (el primer argumento).

Este funcionamiento de rep es especialmente interesante cuando tenemos una tabla de frecuencias.  Todavía no hemos visto, en este blog, como crear tablas en R (aunque nos hemos encontrado con las primeras al hacer Estadística Descriptiva). Pero vamos a suponer que tenemos dos vectores:

valores=c(2,3,5,8,13)
frecuencias=c(5,7,12,2,14)

de manera que el primero contiene los valores (distintos) que aparecen en nuestros datos, y el segundo contiene las frecuencias de esos valores (por supuesto, los dos vectores tendrán la misma longitud). Entonces , para recuperar el conjunto inicial de datos (sin agrupar) partiendo de esta “tabla de frecuencias”, podemos usar así la función rep:

valores=c(2,3,5,8,13)
frecuencias=c(5,7,12,2,14)
datosSinAgrupar=rep(valores,frecuencias)
datosSinAgrupar

Y se obtiene este resultado:

> valores=c(2,3,5,8,13)
> frecuencias=c(5,7,12,2,14)
> datosSinAgrupar=rep(valores,frecuencias)
> datosSinAgrupar
[1]  2  2  2  2  2  3  3  3  3  3  3  3  5  5  5  5  5  5  5  5  5  5  5  5  8  8 13 13
[29] 13 13 13 13 13 13 13 13 13 13 13 13

(Recuerda que la colocación del vector salida depende de la anchura de tu ventana de Rcmdr. Aquí aparecen [1] y [29], pero en tu máquina aparecerán seguramente otros valores, o sólo uno de ellos.)

Finalmente, si lo que queremos es ordenar un vector de números, la función sort viene en nuestra ayuda. Por ejemplo, el resultado de estos comandos:

datos=c(12,13,93,4,42,77,24,15,79,49,20,70,46,33,10,15,62,48,56,88,74,68,8,83,94,48,95,61,15,32)
sort(datos)

es este, en el que aparece el vector ordenado:

> datos=c(12,13,93,4,42,77,24,15,79,49,20,70,46,33,10,15,62,48,56,88,74,68,8,83,94,48,95,61,15,32)
> sort(datos)
[1]  4  8 10 12 13 15 15 15 20 24 32 33 42 46 48 48 49 56 61 62 68 70 74 77 79 83 88 93 94 95

Una aclaración importante. Esta ordenación no afecta en absoluto al vector datos original. Si quieres ordenar el vector y guardar el vector ordenado con el mismo nombre, entonces puedes usar el comando

datos = sort(datos)

pero es importante entender que este es un viaje sin retorno: el vector datos original, sin ordenar, se habrá perdido, y sólo tendremos la versión ordenada. ¡Cuidado con la posible pérdida de información!

Lanzando los dados: la función sample

Disponer de ordenadores, que hacen el trabajo duro por nosotros,  nos permite realizar muchos experimentos en los que poner a prueba nuestras ideas . En los primeros pasos dentro del mundo de la probabilidad, pensamos a menudo en ejemplos que involucran el lanzamiento de un dado.¿Cómo podemos usar R para “lanzar un dado virtual”? La forma más sencilla de hacerlo es utilizar otra función que fabrica vectores, la función sample.El código que usa esta función para simular 100 lanzamientos de un dado es este:

sample(1:6,100,replace=TRUE)

Los tres argumentos de la función sample, en este ejemplo simple, son:

  • el rango de resultados deseados. Por ejemplo, para un dado, del 1 al 6.
  • el número de valores que queremos obtener. Es decir, cuántas veces tiramos el dado, que aquí son 100.
  • en principio, R no admite repeticiones de valores que ya han aparecido (son muestras sin remplazamiento). Si queremos valores repetidos (remplazamiento), debemos indicarlo con el argumento replace=TRUE.

Y el resultado es algo como esto:

> sample(1:6,100,replace=TRUE)
[1] 1 3 1 6 4 6 5 5 5 4 2 2 1 1 4 5 1 3 3 6 1 2 6 6 6 6 1 2 2 5 5 2 3 5 1 2 6 4 6 1 6 2 4 2 5 5 4 1
[49] 3 1 6 5 4 5 5 3 5 2 4 2 6 6 3 3 1 1 1 1 2 1 2 6 4 4 3 5 4 1 3 1 1 3 5 3 6 2 4 6 5 1 5 6 2 3 5 5
[97] 6 3 5 4

En tu máquina los números serán distintos, claro está, porque la función sample está eligiendo aleatoriamente los valores (y en este caso, todos son equiprobables).

Si quieres saber cuántas veces ha aparecido cada uno de los valores, puedes utilizar

unDadoCienVeces=sample(1:6,100,replace=TRUE)
table(unDadoCienVeces)

Le hemos dado un nombre al vector para facilitar las cosas, y para que puedas hacer un experimento: en lugar de 100, escribe 100000, y compara las frecuencias con las anteriores. Además, si no le pusiéramos nombre, al ejecutar el comando verías desfilar ante tus ojos 100000 números en la ventana de resultados…

La función sample permite hacer ajustes mucho más finos que los que hemos visto en este primer ejemplo. Volveremos sobre esto en próximas entradas, pero por el momento nos despedimos de ella con un ejemplo que pretende preparar el terreno y avivar la curiosidad. Ejecuta unas cuantas veces este comando, observa los resultados  y trata de adivinar lo que está pasando (como pista, prob viene de probabilidad):

sample(c(1,2,4,7),100,prob=c(0.2,0.1,0.1,0.6),replace=TRUE)

Seleccionando elementos de vectores

Hemos lanzado 100 veces un dado, y hemos hecho la tabla de frecuencias del vector unDadoCienVeces resultante de ese experimento. Eso nos permite saber cuántas veces ha aparecido, por ejemplo, el número 5. Pero no nos permite saber cuál fue el resultado precisamente en la tirada número 23.  Para averiguar esto, en R escribimos

unDadoCienVeces[23]

De esa forma, usando los corchetes, podemos acceder directamente a cualquier posición del vector. Por ejemplo, el primer y el último elemento del vector son:

unDadoCienVeces[1]
unDadoCienVeces[length(unDadoCienVeces)]

En este caso sabemos que el vector tiene 100 elementos, y para obtener el último elemento podríamos haber escrito

unDadoCienVeces[100]

pero la otra versión (usando length) tiene la ventaja de que seguirá funcionando aunque el vector cambie.

Este recurso de los corchetes es mucho más potente de lo que parece a primera vista. Podemos seleccionar varios elementos a la vez indicando sus posiciones. Para seleccionar los elementos de unDadoCienVeces en las posiciones de la 12 a la 27 hacemos:

unDadoCienVeces[12:27]

y obtenemos un vector con los números que ocupan exactamente esas posiciones. Y si las posiciones que queremos no son consecutivas, basta con indicarlas mediante un vector:

unDadoCienVeces[c(3,8,21,43,56)]

nos proporciona los elementos en esas posiciones del vector. Cuando se combina esto con otras ideas que estamos viendo en esta entrada del blog se pueden conseguir resultados muy útiles. Por ejemplo, ¿qué crees que sucede al hacer esto? (Haz el experimento).

datos=1:100
datos
datos[seq(1,length(datos),3)]

Pero aún podemos ir más lejos. Podemos utilizarlo para seleccionar los elementos del vector que cumplen cierta condición. Por ejemplo, queremos saber cuántos de los resultados en unDadoCienVeces son mayores que 3. En primer lugar, igual que aprendimos (en la segunda sesión con Rcmdr) a sumar una cantidad a un vector, o a elevar un vector al cuadrado, etc., podemos también aplicar una condición lógica, como una desigualdad, a un vector:

unDadoCienVeces > 3

(Añado espacios por claridad; no cambia nada, los añadas o no). La respuesta es un nuevo vector como este (te recuerdo que en tu caso los resultados concretos variarán):

> unDadoCienVeces > 3
[1] FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE
[39] FALSE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
[77]  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE

Cada posición de este vector nos indica si la condición (o suceso) “ser mayor que 3” se cumple en la correspondiente posición del vector unDadoCienVeces. Si lo que queremos es obtener los valores originales del vector, entonces combinamos esta idea con la selección por corchetes:

unDadoCienVeces[ unDadoCienVeces > 3 ]

Y obtendremos una salida parecida a esta

> unDadoCienVeces[ unDadoCienVeces > 3 ]
[1] 6 5 6 6 5 5 4 6 5 4 6 4 6 4 5 5 4 4 5 5 6 4 4 4 4 4 5 4 6 5 5 5 6 4 5 6 5 4 4 5 5 5 5 5 6

con un vector en el que se ha seleccionado todos los elementos (y sólo esos elementos) del vector unDadoCienVeces que cumplen la condición “ser mayor que 3”.

Resumen final

Digerir esta entrada, en la que hemos visto muchos conceptos nuevos, puede llevar algo de tiempo y requiere de un rato de práctica con R. Haz todos los experimentos que quieras: son gratis y no vas a romper nada, claro. Con las técnicas que hemos aprendido, podemos dar los primeros pasos en el mundo de la probabilidad. Pero todavía, si queremos hacer un experimento tan sencillo como lanzar dos dados a la vez, y repetirlo 1000 veces, no disponemos del lenguaje necesario para representar esto en R de la forma más eficaz. Trataremos ese tema en la próxima sesión, en la que conoceremos las tablas de R (ya las hemos intuido, al pedirle a R que haga tablas de frecuencias).

Gracias por la atención.

Quinta sesión con Rcmdr: estadística descriptica con un vector (variable cuantitativa, datos no agrupados)

En esta sesión (quizá un poco más larga que las anteriores), vamos a empezar por fin el trabajo propiamente estadístico. Vamos a comenzar con un objetivo bastante modesto, pero para eso sirven los comienzos. Nuestro objetivo es trabajar con una colección de datos, concretamente un vector de datos, almacenados en un fichero csv. Para fijar ideas, vamos a trabajar con el mismo fichero edades01.csv, que ya usamos en la anterior sesión.

Hacemos esto para que nuestro viaje hacia el mundo de R sea lo más cómodo posible, partiendo de una cierta experiencia trabajando con hojas de cálculo. En una próxima sesión, una vez que hayamos introducido la noción de dataframe en R, volveremos a hacer este estudio, aprovechando las facilidades que proporciona para ello Rcmdr. Pero creo que es mejor empezar por lo más elemental, para que después podamos ver con más claridad donde estamos.

Así, que si no lo recuerdas, este es un buen momento para volver a la anterior sesión con Rcmdr, y ver cómo allí se utilizaba la función scan para leer el fichero csv y almacenarlo en un vector llamado vectorEdades. (Presta especial atención al directorio de trabajo en el que está ubicado el fichero). Ese vector contiene una lista de valores de una variable cuantitativa (de tipo discreto),  no agrupados por intervalos.  El tipo de análisis que vamos a hacer aquí se aplica a este tipo de variables; en otras entradas del blog comentaremos las diferencias que aparecen cuando el análisis se lleva a cabo con otros tipos de variables, o cuando trabajamos con datos agrupados.

Una vez cargados los datos, vamos a utilizar R para ir obteniendo, paso a paso, la descripción estadística de este conjunto de datos. Esta lista es como la lista que uno utiliza para hacer la maleta cuando se va de viaje. Algunas personas prefieren viajar con más equipaje, otras prefieren viajar más ligeros, y en cualquier caso la decisión siempre depende del tipo de viaje. Pero, a medida que uno va viajando, es muy posible que la lista sea cada vez más completa (sin llegar a ser enciclopédica).  La lista que vamos a presentar aquí es, seguramente, un compromiso razonable, un buen punto de partida, que el lector podrá ir completando a medida que gane experiencia viajando por la Estadísitica.

Iremos presentando los valores uno por uno, y al final los resumiremos todos en un fichero de instrucciones R que el lector puede usar como plantilla para este tipo de análisis.

1) Número y rango de los datos

Lo más básico es saber con cuántos datos estamos trabajando. Y para definir el rango de los datos necesitamos el valor mínimo y el máximo. Las instrucciones en R son:

n=length(vectorEdades)
n
min(vectorEdades)
max(vectorEdades)

2) Valores centrales (“medias”)

El siguiente paso puede ser calcular los valores centrales. Los mas sencillos son la media aritmética y la mediana. Los calculamos así:

mean(vectorEdades)
median(vectorEdades)

3) Medidas de dispersión: varianza, desviación típica.
Como siempre, al hablar de medidas de dispersión y software estadístico, se hace necesario hacer una advertencia. Es fundamental que el usuario comprenda si la varianza que está obteniendo es la varianza poblacional, definida así:

S^2=\dfrac{\sum_{j=1}^n (x_j-\bar x)^2}{n}

con n en el denominador, o si se trata de la varianza muestral (o cuasivarianza), definida así:

s^2=\dfrac{\sum_{j=1}^n (x_j-\bar x)^2}{n-1}

con n-1 en el denominador. En el caso de R, la función varsiempre calcula la varianza muestral (la cuasivarianza). ¿Y la varianza poblacional? Bueno, el caso es que, sorprendentemente, R no tiene una función para calcular directamente la varianza poblacional. Pero si se recuerda que la relación entre ambas varianzas es:

S^2=\dfrac{n-1}{n}\, s^2

entonces es muy fácil recuperar la varianza poblacional a partir de la función var de R (y del valor de n, que ya hemos obtenido en el primer punto).

La desviación típica muestral (respectivamente poblacional) es la raíz cuadrada de la varianza muestral (resp. poblacional). Como en el caso de las varianzas, R incluye una función, concretamente la función sd (de standard deviation),  para calcular directamente la desviación típica muestral de los elementos de un vector. Si lo que queremos es calcular la desviación típica poblacional,  calcularemos la varianza poblacional y haremos su raíz cuadrada. Y para eso es necesario saber que la raíz cuadrada se calcula con la función sqrt (de square root).

Con estas ideas, el código adecuado para calcular (y mostrar) ambas varianzas y desviaciones típicas es:

varianzaMuestral=var(vectorEdades)
varianzaMuestral
desvTipicaMuestral=sd(vectorEdades)

varianzaPoblacional=((n-1)/n)*varianzaMuestral
varianzaPoblacional
desvTipicaPoblacional=sqrt(varianzaPoblacional)
desvTipicaPoblacional

4) Medidas de posición: cuartiles y percentiles.

El cálculo de los cuartiles (o percentiles) en R, como en cualquier otro programa estadístico, requiere una discusión más detallada que los  valores anteriores que hemos discutido. Por un lado, si lo único que se desea es obtener la información por defecto que proporciona R, dejándo que sea el programa el que decida la mejor forma de proceder, entonces las cosas son muy sencillas. Basta con escribir

summary(vectorEdades)

Y, para el ejemplo con el que estamos trabajando, obtenemos (en la ventana de resultados) una tabla como esta:

> summary(vectorEdades)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
3.00    6.00    7.00    6.75    8.00   10.00

En la que aparece la información que buscábamos. Por otra parte, en lugar de utilizar summary, si sólo queremos uno de los cuantiles, o si lo que queremos son percentiles, podemos conseguir un control más fino de la respuesta utilizando la función quantile. Esta función necesita dos argumentos: el vector de datos, y un segundo vector en el que aparecen los porcentajes (expresados como tantos por uno) de los percentiles que queremos calcular. Por ejemplo, para pedirle a R que calcule (con el método que más le gusta) los percentiles 5%, 15% y 75% de vectorEdades escribimos este comando:

quantile(vectorEdades,c(0.05,0.15,0.75))
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
3.00    6.00    7.00    6.75    8.00   10.00

y obtenemos

> quantile(vectorEdades,c(0.05,0.15,0.75))
5%  15%  75%
4.95 5.00 8.00

Aquí, como queremos mantener la discusión lo más simple posible, no vamos a entrar en más detalles técnicos sobre las decisiones que R toma para calcular estos números, pero queremos prevenir al lector de que las cosas se complican enseguida: en R hay varias (concretamente, nueve) formas distintas de calcular los cuartiles. Para no entretenernos aquí, dejaremos la información más detallada para otra entrada del blog.

5) Medidas de forma: simetría y apuntamiento.

En muchos cursos de introducción a la Estadística, al presentar la Estadística Descriptiva (típicamente uno de los primeros, si no el primer tema del curso), se describen también las llamadas medidas de forma de un vector de datos. Concretamente, se describen a menudo la simetría y apuntamiento (skewness y kurtosis, en inglés), como formas de saber cuánto se alejan nuestros datos de una distribución normal. Yo prefiero dejar que la distribución normal aparezca de una manera más natural, y más adelante en el curso,  pero en cualquier caso, añado aquí los comandos necesarios para calcular los coeficientes de simetría y apuntamiento  con R. UNa observación, el comando kurtosis no forma parte del núcleo central de comandos de R, y por eso es necesario cargar una librería para aplicarlo (de ahí el comando library(e1071) ):

skewness(vectorEdades)
library(e1071)
kurtosis(vectorEdades)

6) Tablas de frecuencias (absolutas, acumuladas, relativas)

A veces necesitaremos analizar los datos de un vector en forma de tabla de frecuencias. La tabla de frecuencias puede definirse agrupando los valores por intervalos, o simplemente haciendo un recuento de frecuencias para cada uno de los valores diferentes que forman el vector. La decisión depende de la naturaleza de los datos que forman el vector. Para el caso de vectorEdades (el ejemplo que estamos usando), la tabla de frecuencias (absolutas) , sin agrupar por intervalos se obtiene con el sencillo comando:

table(vectorEdades)

Y la respuesta es:

> table(vectorEdades)
vectorEdades
3  4  5  6  7  8  9 10
1  4 16 21 25 23  9  1

Si lo que queremos es una tabla de frecuencias acumuladas, entonces debemos aplicar la función cumsum al resultado de table. Es decir:

cumsum(table(vectorEdades))

que produce el resultado deseado:

> cumsum(table(vectorEdades))
3   4   5   6   7   8   9  10
1   5  21  42  67  90  99 100

Antes de seguir, queremos detenernos en un punto que puede pasar inadvertido. Hasta ahora, en este blog, y desde que empezamos a aprender sobre R, sólo hemos usado variables y vectores. Pero ahora estamos empezando a obtener tablas como resultado de nuestras operaciones. Las tablas son otra forma de almacenar, mostrar  y manipular datos en R. Y aún nos quedan muchos otros objetos en el camino: matrices, listas, datasets, factores… La lista es larga. Pero iremos aprendiendo sobre todos estos recursos sólo a medida que los necesitemos, y sólo los más sencillos ser necesarios en un primer encuentro con la Estadística.

Las tablas, que acabamos de conocer, se comportan en algunos sentidos como los vectores de R. Por ejemplo, si queremos multiplicar por dos todos los elementos de la tabla de frecuencias de vectorEdades, basta con hacer:

2*table(vectorEdades)

Y la respuesta es:

> 2*table(vectorEdades)
vectorEdades
3  4  5  6  7  8  9 10
2  8 32 42 50 46 18  2

Es decir, que lo que aprendimos de aritmética vectorial se extiende a las tablas. Esto tiene una consecuencia inmediata: la tabla de frecuencias relativas se obtiene simplemente dividiendo la tabla de frecuencias absolutas por el número n (el total de datos, es decir length(vectorEdades)). Y para la tabla de frecuencias acumuladas relativas aplicamos cumsum y luego dividimos. Esto es:

table(vectorEdades)/n
cumsum(table(vectorEdades))/n

7) Gráficos (de columnas y boxplot)

En Estadística el tipo de representación gráfica adecuada depende de la naturaleza de los datos. Para los datos que forman el vectorEdades, un gráfico de barras (no confundir con un histograma) puede ser la opción más sencilla y razonable, al tratarse de una variable aleatoria discreta, que ya hemos decidido antes no agrupar en intervalos. El gráfico se obtiene muy fácilmente mediante la función barplot de R, pero aplicada a la tabla de frecuencias de vectorEdades. Es decir:

barplot(table(vectorEdades))

Y el gráfico que se obtiene es este (cuidado; al menos en Windows, el gráfico no aparece en la ventana de Rcmdr, sino en la ventana llamada Rgui, que puede haber quedado cubierta por la de Rcmdr. A primera vista, puede parecer que el comando barplot no ha tenido efecto):

Dos observaciones:

  1. En R disponemos de muchísimas opciones para cambiar los colores, forma, rótulos, etc. de los gráficos. Nos vamos a conformar por el momento con las opciones más simples, pero conviene que el lector sepa que esto no es ni mucho menos el final de la historia con los gráficos de R.
  2. Naturalmente se pueden obtener gráficos similares a partir de las frecuencias relativas, acumuladas, etc.

El gráfico de caja y bigotes es otro tipo habitual de gráficos para datos como los que estamos manejando en este ejemplo. Y obtenerlo con R es también inmediato:

boxplot(table(vectorEdades))

produce como resultado: 

Resumen final

Con esto hemos completado un conjunto de análisis descriptivos del vector de datos bastante completo. Como resumen, y para faciltar el trabajo, agrupamos en este fichero todos los comandos que hemos utilizado. Hemos remplazado vectorEdades por un vectorDatos más genérico, para que el lector no tenga dificultad en reuitlizar este fichero para sus propios datos. Y hemos intercalado comentarios en el listado de comandos para ayudar a recordar el tipo de análisis que se está efectuando. En próximas entradas veremos análisis similares para valores agrupados por intervalos, variables cualitativas, etc.

Gracias por la atención