Del prompt al script

Hasta ahora hemos hecho todo a través del intérprete. De hecho, en la práctica habitual haremos muchas cosas a través del intérprete, pero es una buena idea mantener nuestro código en algún sitio más permanente. Es en este caso en el que algo como RStudio se convierte en una buena idea ya que lo que queremos es algo que nos permita editar archivos de texto con herramientas que hagan sencillo trabajar con el código, pero que al mismo tiempo permita enviar instrucciones al intérprete.

Análisis básico de datos

Empecemos leyendo datos de Internet:

affairs <- read.csv("http://koaning.io/theme/data/affairs.csv")

La base de datos contiene información sobre los affairs de 601 políticos y alguna información sociodemográfica. Una descripción de las variables (y algunos resultados curiosos) se pueden encontrar en Fair, R. (1977) “A note on the computation of the tobit estimator”, Econometrica, 45, 1723-1727.

Empecemos echándole un vistazo a los datos. Por ejemplo, podemos imprimir en pantalla las primeras 6 líneas del archivo usando la función head:

head(affairs)

También podemos echarle un vistazo a algunos descriptivos usando la función summary aplicada a variables:

summary(affairs$nbaffairs)

y

summary(affairs$child)

Es importante notar que summary hace cosas diferentes dependiendo el input que reciba. En el primer caso, summary recibe una variable numérica y produce la media y los cuartiles. En el segundo caso, summary recibe un factor (una variable categórica) y produce una tabla de frecuencias. Este es un patrón de comportamiento que veremos muy a menudo en R.

Podemos ser más específicos y pedir una tabla de frecuencias con la función:

table(affairs$child)

Para transformar la tabla a proporciones podemos tomar diferentes rutas. Una de ellas sería hacerlo manualmente, dividiendo las frecuencias por el total de casos, tomando la longitud de la columna:

table(affairs$child)/length(affairs$child)

No hemos creado ninguna variable, sino que hemos realizado las operaciones necesarias en una única llamada. La segunda opción es “componer” dos funciones:

prop.table(table(affairs$child))

En ese caso, el resultado de table pasa a la función prop.table, la cual transforma frecuencias a proporciones.

Hypothesis testing

Podemos empezar ahora a analizar los datos. Por ejemplo, imaginemos que estamos interesados en la diferencia en la media de aventuras que tienen los políticos dependiendo de si tienen hijos o no. La media muestral para cada grupo puede calcularse como:

mean(affairs$nbaffairs[affairs$child == "no"])
mean(affairs$nbaffairs[affairs$child == "yes"])

Un contraste de t puede hacerse de varias formas. La forma más natural es pasando variables a la función t.test. Por ejemplo, si quisiésemos contrastar una variable contra la hipótesis nula habitual, usaríamos:

t.test(affairs$nbaffairs)

También podríamos contrastar la igualdad de medias pasando dos vectores a la función:

t.test(affairs$nbaffairs[affairs$child == "no"], affairs$nbaffairs[affairs$child == "yes"])

Lo que es importante a tener en cuenta es que el segundo vector es un argumento opcional y que la función hace algo diferente cuando el segundo argumento existe. Echemos un vistazo a la documentación de t.test:

?t.test

Vemos que existen dos métodos separados para interactuar con t.test: el que acabamos de usar, pasando los argumentos x y tal vez y, y otro que usa una formula. Las fórmulas son muy imporatntes en R y las usaréis a menudo en la segunda mitad del curso:

my_test <- t.test(nbaffairs ~ child, data=affairs)
my_test

La parte izquierda de la ecuación es la variable que queremos contrastar para cada uno de los valores indicados en el lado derecho. El argumento data indica dónde viven esas variables: son columnas del data.frame affairs. La interfaz de fórmula cobra más sentido si pensamos el contraste como si fuese un modelo lineal. En ese caso caso, el lado izquierdo de la fórmula, separado por el símbolo ~, se corresponde con el lado izquierdo de la ecuación. La parte derecha es una variable ficticia (un factor) que divide la muestra en dos grupos.

Podemos inspeccionar los contenidos de my_test usando la función str:

str(my_test)

Tened en cuenta que my_test es una lista que contiene toda la información relativa al contraste que hemos estimado. Esta información está accesible a través de una lista:

my_test$statistic
my_test$conf.int
my_test$estimate

Quizás sea ahora un buen momento para volver a la documentación y comparar el resultado del contraste con lo que dice la sección “Value” en la ayuda.

LS0tIAp0aXRsZTogIlJ1ZGltZW50b3MgZGUgYW7DoWxpc2lzIGVzdGFkw61zdGljbyIKZGF0ZTogImByIGZvcm1hdChTeXMudGltZSgpLCAnJUIgJWQsICVZJylgIgpvdXRwdXQ6IAogIG1kX2RvY3VtZW50CmV4dGVuc2lvbjogZm9vdG5vdGVzCi0tLQoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0UsIGNhY2hlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQoZXZhbCA9IEZBTFNFKSAKYGBgCgojIyBEZWwgX3Byb21wdF8gYWwgX3NjcmlwdF8KCkhhc3RhIGFob3JhIGhlbW9zIGhlY2hvIHRvZG8gYSB0cmF2w6lzIGRlbCBpbnTDqXJwcmV0ZS4gRGUgaGVjaG8sIGVuIGxhIHByw6FjdGljYQpoYWJpdHVhbCBoYXJlbW9zIG11Y2hhcyBjb3NhcyBhIHRyYXbDqXMgZGVsIGludMOpcnByZXRlLCBwZXJvIGVzIHVuYSBidWVuYSBpZGVhCm1hbnRlbmVyIG51ZXN0cm8gY8OzZGlnbyBlbiBhbGfDum4gc2l0aW8gbcOhcyBwZXJtYW5lbnRlLiBFcyBlbiBlc3RlIGNhc28gZW4gZWwgcXVlCmFsZ28gY29tbyBSU3R1ZGlvIHNlIGNvbnZpZXJ0ZSBlbiB1bmEgYnVlbmEgaWRlYSB5YSBxdWUgbG8gcXVlIHF1ZXJlbW9zIGVzIGFsZ28KcXVlIG5vcyBwZXJtaXRhIGVkaXRhciBhcmNoaXZvcyBkZSB0ZXh0byBjb24gaGVycmFtaWVudGFzIHF1ZSBoYWdhbiBzZW5jaWxsbwp0cmFiYWphciBjb24gZWwgY8OzZGlnbywgcGVybyBxdWUgYWwgbWlzbW8gdGllbXBvIHBlcm1pdGEgZW52aWFyIGluc3RydWNjaW9uZXMgYWwKaW50w6lycHJldGUuCgojIyBBbsOhbGlzaXMgYsOhc2ljbyBkZSBkYXRvcwoKRW1wZWNlbW9zIGxleWVuZG8gZGF0b3MgZGUgSW50ZXJuZXQ6CgpgYGB7cn0KYWZmYWlycyA8LSByZWFkLmNzdigiaHR0cDovL2tvYW5pbmcuaW8vdGhlbWUvZGF0YS9hZmZhaXJzLmNzdiIpCmBgYAoKTGEgYmFzZSBkZSBkYXRvcyBjb250aWVuZSBpbmZvcm1hY2nDs24gc29icmUgbG9zIF9hZmZhaXJzXyBkZSA2MDEgcG9sw610aWNvcyB5CmFsZ3VuYSBpbmZvcm1hY2nDs24gc29jaW9kZW1vZ3LDoWZpY2EuIFVuYSBkZXNjcmlwY2nDs24gZGUgbGFzIHZhcmlhYmxlcyAoeSBhbGd1bm9zCnJlc3VsdGFkb3MgY3VyaW9zb3MpIHNlIHB1ZWRlbiBlbmNvbnRyYXIgZW4gRmFpciwgUi4gKDE5NzcpICJBIG5vdGUgb24gdGhlCmNvbXB1dGF0aW9uIG9mIHRoZSB0b2JpdCBlc3RpbWF0b3IiLCBFY29ub21ldHJpY2EsIDQ1LCAxNzIzLTE3MjcuCgpFbXBlY2Vtb3MgZWNow6FuZG9sZSB1biB2aXN0YXpvIGEgbG9zIGRhdG9zLiBQb3IgZWplbXBsbywgcG9kZW1vcyBpbXByaW1pciBlbgpwYW50YWxsYSBsYXMgcHJpbWVyYXMgNiBsw61uZWFzIGRlbCBhcmNoaXZvIHVzYW5kbyBsYSBmdW5jacOzbiBgaGVhZGA6CgpgYGB7cn0KaGVhZChhZmZhaXJzKQpgYGAKClRhbWJpw6luIHBvZGVtb3MgZWNoYXJsZSB1biB2aXN0YXpvIGEgYWxndW5vcyBkZXNjcmlwdGl2b3MgdXNhbmRvIGxhIGZ1bmNpw7NuCmBzdW1tYXJ5YCBhcGxpY2FkYSBhIHZhcmlhYmxlczoKCmBgYHtyfQpzdW1tYXJ5KGFmZmFpcnMkbmJhZmZhaXJzKQpgYGAKCnkKCmBgYHtyfQpzdW1tYXJ5KGFmZmFpcnMkY2hpbGQpCmBgYAoKRXMgaW1wb3J0YW50ZSBub3RhciBxdWUgYHN1bW1hcnlgIGhhY2UgY29zYXMgZGlmZXJlbnRlcyBkZXBlbmRpZW5kbyBlbCBpbnB1dCBxdWUKcmVjaWJhLiBFbiBlbCBwcmltZXIgY2FzbywgYHN1bW1hcnlgIHJlY2liZSB1bmEgdmFyaWFibGUgbnVtw6lyaWNhIHkgcHJvZHVjZSBsYQptZWRpYSB5IGxvcyBjdWFydGlsZXMuIEVuIGVsIHNlZ3VuZG8gY2FzbywgIGBzdW1tYXJ5YCByZWNpYmUgdW4gZmFjdG9yICh1bmEKdmFyaWFibGUgY2F0ZWfDs3JpY2EpIHkgcHJvZHVjZSB1bmEgdGFibGEgZGUgZnJlY3VlbmNpYXMuIEVzdGUgZXMgdW4gcGF0csOzbiBkZQpjb21wb3J0YW1pZW50byBxdWUgdmVyZW1vcyBtdXkgYSBtZW51ZG8gZW4gYFJgLgoKUG9kZW1vcyBzZXIgbcOhcyBlc3BlY8OtZmljb3MgeSBwZWRpciB1bmEgdGFibGEgZGUgZnJlY3VlbmNpYXMgY29uIGxhIGZ1bmNpw7NuOgoKYGBge3J9CnRhYmxlKGFmZmFpcnMkY2hpbGQpCmBgYAoKUGFyYSB0cmFuc2Zvcm1hciBsYSB0YWJsYSBhIHByb3BvcmNpb25lcyBwb2RlbW9zIHRvbWFyIGRpZmVyZW50ZXMgcnV0YXMuIFVuYSBkZQplbGxhcyBzZXLDrWEgaGFjZXJsbyBtYW51YWxtZW50ZSwgZGl2aWRpZW5kbyBsYXMgZnJlY3VlbmNpYXMgcG9yIGVsIHRvdGFsIGRlCmNhc29zLCB0b21hbmRvIGxhIGxvbmdpdHVkIGRlIGxhIGNvbHVtbmE6CgpgYGB7cn0KdGFibGUoYWZmYWlycyRjaGlsZCkvbGVuZ3RoKGFmZmFpcnMkY2hpbGQpCmBgYAoKTm8gaGVtb3MgY3JlYWRvIG5pbmd1bmEgdmFyaWFibGUsIHNpbm8gcXVlIGhlbW9zIHJlYWxpemFkbyBsYXMgb3BlcmFjaW9uZXMKbmVjZXNhcmlhcyBlbiB1bmEgw7puaWNhIGxsYW1hZGEuIExhIHNlZ3VuZGEgb3BjacOzbiBlcyAiY29tcG9uZXIiIGRvcyBmdW5jaW9uZXM6CgpgYGB7cn0KcHJvcC50YWJsZSh0YWJsZShhZmZhaXJzJGNoaWxkKSkKYGBgCgpFbiBlc2UgY2FzbywgZWwgcmVzdWx0YWRvIGRlIGB0YWJsZWAgcGFzYSBhIGxhIGZ1bmNpw7NuIGBwcm9wLnRhYmxlYCwgbGEgY3VhbAp0cmFuc2Zvcm1hIGZyZWN1ZW5jaWFzIGEgcHJvcG9yY2lvbmVzLgoKIyMgSHlwb3RoZXNpcyB0ZXN0aW5nCgpQb2RlbW9zIGVtcGV6YXIgYWhvcmEgYSBhbmFsaXphciBsb3MgZGF0b3MuIFBvciBlamVtcGxvLCBpbWFnaW5lbW9zIHF1ZSBlc3RhbW9zCiBpbnRlcmVzYWRvcyBlbiBsYSBkaWZlcmVuY2lhIGVuIGxhIG1lZGlhIGRlIGF2ZW50dXJhcyBxdWUgdGllbmVuIGxvcyBwb2zDrXRpY29zCiBkZXBlbmRpZW5kbyBkZSBzaSB0aWVuZW4gaGlqb3MgbyBuby4gTGEgbWVkaWEgbXVlc3RyYWwgcGFyYSBjYWRhIGdydXBvIHB1ZWRlCiBjYWxjdWxhcnNlIGNvbW86CgpgYGB7ciBldmFsPUZBTFNFfQptZWFuKGFmZmFpcnMkbmJhZmZhaXJzW2FmZmFpcnMkY2hpbGQgPT0gIm5vIl0pCm1lYW4oYWZmYWlycyRuYmFmZmFpcnNbYWZmYWlycyRjaGlsZCA9PSAieWVzIl0pCmBgYAoKVW4gY29udHJhc3RlIGRlIF90XyBwdWVkZSBoYWNlcnNlIGRlIHZhcmlhcyBmb3JtYXMuIExhIGZvcm1hIG3DoXMgbmF0dXJhbCBlcwpwYXNhbmRvIHZhcmlhYmxlcyBhIGxhIGZ1bmNpw7NuIGB0LnRlc3RgLiBQb3IgZWplbXBsbywgc2kgcXVpc2nDqXNlbW9zIGNvbnRyYXN0YXIKdW5hIHZhcmlhYmxlIGNvbnRyYSBsYSBoaXDDs3Rlc2lzIG51bGEgaGFiaXR1YWwsIHVzYXLDrWFtb3M6CgpgYGB7cn0KdC50ZXN0KGFmZmFpcnMkbmJhZmZhaXJzKQpgYGAKClRhbWJpw6luIHBvZHLDrWFtb3MgY29udHJhc3RhciBsYSBpZ3VhbGRhZCBkZSBtZWRpYXMgcGFzYW5kbyBfZG9zXyB2ZWN0b3JlcyBhIGxhIGZ1bmNpw7NuOgpgYGB7cn0KdC50ZXN0KGFmZmFpcnMkbmJhZmZhaXJzW2FmZmFpcnMkY2hpbGQgPT0gIm5vIl0sIGFmZmFpcnMkbmJhZmZhaXJzW2FmZmFpcnMkY2hpbGQgPT0gInllcyJdKQpgYGAKCkxvIHF1ZSBlcyBpbXBvcnRhbnRlIGEgdGVuZXIgZW4gY3VlbnRhIGVzIHF1ZSBlbCBzZWd1bmRvIHZlY3RvciBlcyB1biBhcmd1bWVudG8Kb3BjaW9uYWwgeSBxdWUgbGEgZnVuY2nDs24gaGFjZSBhbGdvIGRpZmVyZW50ZSBjdWFuZG8gZWwgc2VndW5kbyBhcmd1bWVudG8KZXhpc3RlLiBFY2hlbW9zIHVuIHZpc3Rhem8gYSBsYSBkb2N1bWVudGFjacOzbiBkZSAgYHQudGVzdGA6CgpgYGB7ciBldmFsPUZBTFNFfQo/dC50ZXN0CmBgYAoKVmVtb3MgcXVlIGV4aXN0ZW4gZG9zIG3DqXRvZG9zIHNlcGFyYWRvcyBwYXJhIGludGVyYWN0dWFyIGNvbiBgdC50ZXN0YDogZWwgcXVlCmFjYWJhbW9zIGRlICB1c2FyLCBwYXNhbmRvIGxvcyBhcmd1bWVudG9zIGB4YCB5IHRhbCB2ZXoKYHlgLCB5IG90cm8gcXVlIHVzYSB1bmEgYGZvcm11bGFgLiBMYXMgZsOzcm11bGFzIHNvbiBtdXkgaW1wb3JhdG50ZXMgZW4gYFJgIHkgbGFzCnVzYXLDqWlzIGEgbWVudWRvIGVuIGxhIHNlZ3VuZGEgbWl0YWQgZGVsIGN1cnNvOgoKYGBge3J9Cm15X3Rlc3QgPC0gdC50ZXN0KG5iYWZmYWlycyB+IGNoaWxkLCBkYXRhPWFmZmFpcnMpCm15X3Rlc3QKYGBgCgpMYSBwYXJ0ZSBpenF1aWVyZGEgZGUgbGEgZWN1YWNpw7NuIGVzIGxhIHZhcmlhYmxlIHF1ZSBxdWVyZW1vcyBjb250cmFzdGFyIHBhcmEKY2FkYSB1bm8gZGUgbG9zIHZhbG9yZXMgaW5kaWNhZG9zIGVuIGVsIGxhZG8gZGVyZWNoby4gRWwgYXJndW1lbnRvIGBkYXRhYCBpbmRpY2EKZMOzbmRlIHZpdmVuIGVzYXMgdmFyaWFibGVzOiBzb24gY29sdW1uYXMgZGVsIGBkYXRhLmZyYW1lYCBgYWZmYWlyc2AuIExhIGludGVyZmF6CmRlIGbDs3JtdWxhIGNvYnJhIG3DoXMgc2VudGlkbyBzaSBwZW5zYW1vcyBlbCBjb250cmFzdGUgY29tbyBzaSBmdWVzZSB1biBtb2RlbG8KbGluZWFsLiBFbiBlc2UgY2FzbyBjYXNvLCBlbCBsYWRvIGl6cXVpZXJkbyBkZSBsYSBmw7NybXVsYSwgc2VwYXJhZG8gcG9yIGVsCnPDrW1ib2xvIGB+YCwgc2UgY29ycmVzcG9uZGUgY29uIGVsIGxhZG8gaXpxdWllcmRvIGRlIGxhIGVjdWFjacOzbi4gTGEgcGFydGUKZGVyZWNoYSBlcyB1bmEgdmFyaWFibGUgZmljdGljaWEgKHVuIGZhY3RvcikgcXVlIGRpdmlkZSBsYSBtdWVzdHJhIGVuIGRvcwpncnVwb3MuIAoKUG9kZW1vcyBpbnNwZWNjaW9uYXIgbG9zIGNvbnRlbmlkb3MgZGUgYG15X3Rlc3RgIHVzYW5kbyBsYSBmdW5jacOzbiBgc3RyYDoKCmBgYHtyfQpzdHIobXlfdGVzdCkKYGBgCgpUZW5lZCBlbiBjdWVudGEgcXVlIGBteV90ZXN0YCBlcyB1bmEgbGlzdGEgcXVlIGNvbnRpZW5lIHRvZGEgbGEgaW5mb3JtYWNpw7NuCnJlbGF0aXZhIGFsIGNvbnRyYXN0ZSBxdWUgaGVtb3MgZXN0aW1hZG8uIEVzdGEgaW5mb3JtYWNpw7NuIGVzdMOhIGFjY2VzaWJsZSBhCnRyYXbDqXMgZGUgdW5hIGxpc3RhOgoKYGBge3IgZXZhbD1GQUxTRX0KbXlfdGVzdCRzdGF0aXN0aWMKbXlfdGVzdCRjb25mLmludApteV90ZXN0JGVzdGltYXRlCmBgYAoKUXVpesOhcyBzZWEgYWhvcmEgdW4gYnVlbiBtb21lbnRvIHBhcmEgdm9sdmVyIGEgbGEgZG9jdW1lbnRhY2nDs24geSBjb21wYXJhciBlbApyZXN1bHRhZG8gZGVsIGNvbnRyYXN0ZSBjb24gbG8gcXVlIGRpY2UgbGEgc2VjY2nDs24gIlZhbHVlIiBlbiBsYSBheXVkYS4KCg==