Guía R

Cátedra de Probabilidad y Estadística

R es un lenguaje de programación libre para el análisis estadístico y gráfico de datos. R Está

disponible para todos los sistemas operativos y puede descargarse libremente desde el sitio web
http://www.r-project.org/. De la misma forma, RStudio puede descargarse libremente desde el
sitio web https://www.rstudio.com/.

Índice General
1 La ayuda de R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2 Definir el directorio de trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

3 Cargar y leer paquetes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

4 Leer datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

5 Visualización de datos y selección de variables . . . . . . . . . . . . . . . . . . . . . 4

6 Matrices y vectores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
6.1 Vectores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
6.2 Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

7 Funciones básicas y estadísticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

8 Resumen numérico y gráfico de los datos . . . . . . . . . . . . . . . . . . . . . . . . 11

8.1 Una variable cualitativa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
8.2 Dos variables cualitativas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
8.3 Una variable cuantitativa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
8.4 Una variable cuantitativa vs. una cualitativa . . . . . . . . . . . . . . . . . . . . . . . 16
8.5 Dos variables cuantitativas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
8.6 Una variable cuantitativa vs. dos cualitativas (factores) . . . . . . . . . . . . . . . . . 21
8.7 Más sobre gráficos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

1
9 Funciones de distribución . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
9.1 Variables aleatorias discretas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
9.2 Variables aleatorias continuas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

10 Intervalos de confianza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
10.1 Una muestra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
10.2 Dos muestras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
10.3 Intervalos de confianza para una proporción . . . . . . . . . . . . . . . . . . . . . . . 31

11 Test de hipótesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

12 ANOVA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

13 Regresión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

2
Instalación de R y RStudio
Para comenzar a utilizar R, el primer paso es descargarlo de la página http://www.r-project.org/
y luego proceder a su instalación.
Una vez instalado R se podrá instalar RStudio, un entorno de desarrollo integrado para el lenguaje
R que lo vuelve más fácil de utilizar. Se lo debe descargar de la página https://www.rstudio.com/.

1. La ayuda de R
Dentro de RStudio, una de las ventanas contiene la pestaña "Help", en la que se podrá ingresar el
nombre de cualquier comando que se quiera utilizar. Es muy recomendable acostumbrarse a usarla,
ya que cada explicación viene acompañada de ejemplos de uso.
En las ventanas de Consola y Fuente, también se podrá utilizar la ayuda, con los comandos help()
o ?. Por ejemplo si queremos buscar la ayuda para instalar paquetes, podremos hacerlo con los
comandos
help(install.packages)
ó
?install.packages

2. Definir el directorio de trabajo

Por defecto RStudio trabaja en un directorio específico que se crea al instalar el programa, y podemos
saber cuál es a través del comando getwd(). Por ejemplo,

> getwd()
[1] "C:/Documents and Settings/R"

Para cambiar el directorio de trabajo a través de la barra de menú podemos hacerlo en Session > > Set Wor
e ir a la carpeta deseada.
También podemos hacerlo desde la ventana de comando utilizando el comando setwd especificando
el directorio donde deseamos trabajar. Por ejemplo,

> setwd(C:/Usuario/PyE)

3. Cargar y leer paquetes

Para poder utilizar un paquete específico es necesario leerlo. Para hacerlo utilizamos el comando
library. Por ejemplo, para leer el paquete Sleuth2 utilizamos: library(Sleuth2).

3
La primera vez que se lo quiera utilizar, será necesario instalarlo. Para hacerlo en RStudio debemos
ir a: Tools > > Install Packages... allí seleccionar si se lo desea instalar desde internet o desde
un archivo (en general lo hacemos desde internet a través de un repositorio (CRAN)), escribir el
nombre del paquete, por ejemplo, Sleuth2 e instalarlo presionando Install.
Si queremos instalarlo desde las ventanas de Consola y Fuente, se puede utilizar el comando
install.packages, con el ejemplo anterior podemos hacerlo escribiendo install.packages(’’Sleuth2’’)

La lista de paquetes que será necesario instalar para el desarrollo de esta materia es:

readr, car, MASS, mosaic, multcomp, Sleuth2, tseries, lattice

4. Leer datos
Para leer y cargar los datos tenemos que fijarnos que tipo de extensión tiene el archivo. De acuerdo
a la extensión es la función que usamos.

Archivos desde librerías de R: Algunas librerías de R contienen varios conjuntos de datos. Por
ejemplo, la librería Sleuth2 contiene, además de otros, un conjunto de datos llamado case0801,
para leer estos datos primero leemos la librería como en 3. y luego leemos los datos como sigue:

library(Sleuth2)
datos = case0801

Importar archivos: En RStudio se puede ir a Files > > Import Dataset y elegir el tipo de
archivo. Hay que tener en cuenta

• Tipo de separación de los datos (coma, punto y coma, espacio, tabulación)

• Si tiene encabezado
• Cómo tratar los datos faltantes (si es que los hay)
• El nombre que tendrán los datos

Los datos también se pueden cargar desde la Consola o la Fuente.

• Archivos de extensión .R o .txt:

> data1 = read.table(’salchichas.txt’, header = TRUE)

> data2 = read.table(’dieta.R’, header = TRUE)

4
 lee los archivos .R y .txt y los guarda en data1 y data2 respectivamente. Con el argumento
header estamos indicando que queremos que también se lean y guarden los títulos que
identifican cada variable. Este comando se utiliza siempre y cuando los datos tengan
encabezado, si no los tiene se debe poner header = FALSE.
• Archivos de valores separados por comas (.csv), utilizar el siguiente comando:

> data3 = read.csv(’alumnos.csv’, header = TRUE, sep = ’;’)

> data4 = read.csv(’animales.csv’, header = TRUE, sep = ’,’)

> data5 = read.csv(’pennstate.csv’, header = TRUE, sep = ’,’)

Cuando los datos están ordenados usando las columnas de excel, como en el caso de
alumnos.csv, debemos declarar sep = ’;’ para que [R] entienda que cada columna es
una variable. En cambio, si solamente se separan utilizando las comas, como en el archivo
pennstate.csv, declaramos sep = ’,’. Si los datos están separados por una tabulación
debemos usar sep = ’\t’.
• Archivos de Excel (.xls, .xlsx), utilizando el paquete readxl y el comando read_excel().
• Datos faltantes: Hay veces en que los datos tienen “datos faltantes” por lo que en ese
caso debemos usar sep = ’\t’. Además si el dato faltante corresponde a una variable
cualitativa, debemos agregar la instrucción na.string = ’’ como sigue,

> data3 = read.csv(’alumnos.csv’, header = TRUE, sep = ’\t’, na.string = ’’

5. Visualización de datos y selección de variables

Para poder visualizar los datos, una vez cargados sólo escribimos en [R] su nombre.

> data1
CALORIAS SODIO TIPO
1 186 495 A
2 181 477 A
3 176 425 A
4 149 322 A
......................

Para poder saber el nombre de las columnas hacemos:

> names(data1)
[1] "CALORIAS" "SODIO" "TIPO"
5
Para poder trabajar con una columna específica se utiliza:

> data1$CALORIAS
[1] 186 181 176 149 184 190 158 139 175 148 152 111 141 153 190 157 131 149 135 132
[21] 173 191 182 190 172 147 146 139 175 136 179 153 107 195 135 140 138 129 132 102
[41] 106 94 102 87 99 170 113 135 142 86 143 152 146 144

pero si queremos directamente utilizar el nombre de la variable luego de leer los datos debemos hacer:

> attach(data1)

y entonces ya podemos llamar a cada variable directamente por su nombre, por ejemplo:

> CALORIAS
[1] 186 181 176 149 184 190 158 139 175 148 152 111 141 153 190 157 131 149 135 132
[21] 173 191 182 190 172 147 146 139 175 136 179 153 107 195 135 140 138 129 132 102
[41] 106 94 102 87 99 170 113 135 142 86 143 152 146 144

Para ver el comienzo de los datos (títulos y primeras filas):

> head(data1)
CALORIAS SODIO TIPO
1 186 495 A
2 181 477 A
3 176 425 A
4 149 322 A
5 184 482 A
6 190 587 A

Para ver la descripción de los mismos:

> str(data1)
’data.frame’: 54 obs. of 3 variables:
$ CALORIAS: int 186 181 176 149 184 190 158 139 175 148 ...
$ SODIO : int 495 477 425 322 482 587 370 322 479 375 ...
$ TIPO : Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...

Para saber si una variable tiene datos faltantes usamos la función is.na que es una función lógica que
devuelve TRUE en la posición donde falta el dato y FALSE en las restantes. Por ejemplo:

> a = c(1,2,NA,4,5,NA)
> a
[1] 1 2 NA 4 5 NA

6
> is.na(a)
[1] FALSE FALSE TRUE FALSE FALSE TRUE

Si directamente queremos saber si una variable tiene algún dato faltante podemos hacer:

> sum(is.na(a) == FALSE)

[1] 4

Para saber el índice de un elemento utilizamos el comando which como sigue:

> a = c(-1,2,4,-5)
> which(a>0)
[1] 2 3

Si por ejemplo nos interesa quedarnos con aquellos valores del sodio que corresponde a las salchichas tipo A,
lo podemos hacer como sigue y lo guardamos en la nueva variable SOA:

> SOA = SODIO[TIPO == ’A’]

> SOA
[1] 495 477 425 322 482 587 370 322 479 375 330 300 386 401 645 440 317 319 298
[20] 253

Otro ejemplo más,

> CALA = CALORIAS[SODIO >= 400]

Si nos interesa quedarnos con aquellos valores del sodio que corresponde a las salchichas que no son de tipo
A, lo podemos hacer como sigue:

> SOA = SODIO[TIPO != ’A’]

Si nos interesa quedarnos con aquellos valores de calorías correspondientes a las salchichas de tipo A con un
nivel de sodio mayor o igual a 400, lo podemos hacer como sigue:

> CALORIAS[SODIO > = 400 & TIPO == ’A’]

Si deseamos quedarnos con aquellos valores de calorías correspondientes a las salchichas de tipo C o con un
nivel de sodio menor a 400, lo podemos hacer como sigue:

> CALORIAS[SODIO < 400 | TIPO == ’C’]

7
6. Matrices y vectores
6.1. Vectores
Vamos a ver algunas formas de generar vectores. Se puede construir un vector de tipo numérico, lógico o
caracter. Algunos ejemplos:

> vector1 = c(74, 122, 235, 111, 292, 111, 211, 133, 156, 79)

> vector2 = c(’Siempre’, ’casi siempre’, ’a veces’, ’raramente’,’Nunca’)

con la expresión c() concatenamos los valores que le ingresamos para crear el vector. Otras funciones para
la generación de vectores:

seq(from = a,to = b,by = c): genera una secuencia de números desde a hasta b con incrementos
de c.

a:b: lo mismo que antes con incremento 1.

rev(a:b): genera una secuencia de números desde b hasta a.

rep(x,v): repite el elemento x un número v de veces.

Para acceder a alguna posición del vector utilizamos la expresión vector[]. Por ejemplo, si queremos acceder
a la posición 3 del vector1 hacemos:

> vector1[3]
[1] 235

Para saber la longitud que tiene un vector utilizamos el comando length como sigue:

> length(vector1)
[1] 10

6.2. Matrices
Para crear una matriz la forma más simple en R, es mediante la función matrix(). A continuación, se muestra
como crear una matriz llamada M de dimensión 3x2:

> M = matrix(nrow = 2, ncol = 3)

> M
[,1] [,2] [,3]
[1,] NA NA NA
[2,] NA NA NA

8
Ya que no se ha definido ningún elemento, R nos lo informa con NA, para crear una matriz con elementos,
por ejemplo, mediante columnas: (1, 2, 3, 4, 5, 6):

> M = matrix(c(1,2,3,4,5,6), nrow = 2, ncol = 3)

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

En caso que queramos crearlo por filas:

> M = matrix(c(1,2,3,4,5,6), nrow = 2, ncol = 3, byrow = TRUE)

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

Y por último, si queremos añadirles nombres a los índices de columnas: Naranjas, Platanos y Melon, y en
las filas: Supermercado y Tienda, se muestra a continuación:

> M = matrix(c(1,2,3,4,5,6), nrow = 2, ncol = 3, byrow = TRUE,

+dimnames = list(c(’Supermercado’, ’Tienda’),
+c(’Naranjas’,’Platanos’, ’Melon’)))

> M
Naranjas Platanos Melon
Supermercado 1 2 3
Tienda 4 5 6

Para acceder a alguna posición de la matriz utilizamos la expresión M[,]. Por ejemplo, si queremos acceder
a la posición (2,3) de M hacemos:

> M[2,3]
[1] 6

Importante: Observar la diferencia entre [] y (). Los primeros sirven para acceder a un elemento del vector
(o de una matriz) y los segundos sirven sólo para funciones, como por ejemplo para crear vectores utilizamos
la función concatenar c() o para crear una matriz matrix().

Veamos algunas funciones sobre las matrices:

dim(M): devuelve las dimensiones de la matriz M.

colnames(M): devuelve el nombre de las columnas de la matriz M.

9
 rownames(M): devuelve el nombre de las filas de la matriz M.

M[ , ]: accede a elementos dentro de la matriz M.

t(M): devuelve la matriz M transpuesta.

apply: aplica una función sobre las filas o columnas de una matriz.

cbind(a,b,...): crea una matriz con los vectores o matrices a, b, ... como columnas.

rbind(a,b,...): crea una matriz con los vectores o matrices a, b, ... como filas.

Veamos algunos ejemplos, utilizando la matriz M creada antes:

> M[1,]
Naranjas Platanos Melon
1 2 3

> M[,1:2]
Naranjas Platanos
Supermercado 1 2
Tienda 4 5

> cbind(M, c(0,0))

Naranjas Platanos Melon
Supermercado 1 2 3 0
Tienda 4 5 6 0

> a = c(5, 27, 32)

> b = 4:6
> c = c(15, 95, 4)
> S = rbind(a,b,c)
> S
[,1] [,2] [,3]
a 5 27 32
b 4 5 6
c 15 95 4

Una ventaja que posee este software es que podemos extraer de manera muy sencilla filas o columnas que
no son contiguas. Por ejemplo, si queremos extraer las columnas 1 y 3 de M hacemos:

> M[,c(1,3)]
Naranjas Melon
Supermercado 1 3
Tienda 4 6
10
O si queremos todas las columnas menos la 1, hacemos:

> M[,-1]
Platanos Melon
Supermercado 2 3
Tienda 5 6

O mejor aún, si queremos todas las filas de M menos la 2 y 3 hacemos:

> M[-c(2,3),]
Naranjas Platanos Melon
1 2 3

por supuesto que en este caso sencillo esto es lo mismo que hacer M[1,].

7. Funciones básicas y estadísticas

A continuación se enumeran algunas funciones útiles. Se pueden buscar ejemplos de utilización con la Ayuda
de R.

Las funciones min(x) y max(x) calculan el mínimo y el máximo, respectivamente, del elemento x. En
caso de existir datos faltantes o perdidos es necesario agregar la opción na.rm = TRUE.

La función range(x) devuelve un vector que contiene el mínimo y el máximo de x.

La función sum(x) devuelve la suma de todos los elementos de x. En caso de existir datos perdidos es
necesario agregar la opción na.rm = TRUE como antes.

La función mean(x) devuelve el promedio de todos los elementos de x.:

La función median(x) devuelve la mediana de x.

La función var(x) devuelve la varianza de x.

La función sd(x) devuelve el desvío estándar de x.

La función quantile(x) devuelve el mínimo, el primer, segundo y tercer cuantil, y el máximo de x.

La función IQR(x) el rango intercuartílico de x, es decir, la diferencia entre el tercer y primer cuantil.

Obervarción: Para las funciones estadísticas, en caso de existir datos perdidos es necesario agregar la opción
na.rm = TRUE como antes.

11
8. Resumen numérico y gráfico de los datos
8.1. Una variable cualitativa
Resumen numérico: Tabla de frecuencias absolutas y relativas. Utilizamos la función table y prop.table
respectivamente. Por ejemplo,

> table(TIPO)
TIPO
A B C
20 17 17

> prop.table(table(TIPO))
TIPO
A B C
0.3703704 0.3148148 0.3148148

Resumen gráfico: Gráficos de tortas y barras.

Para realizar un graáfico de torta utilizamos la función pie como sigue:

> pie(table(TIPO), clockwise = T, main = ’Tipo de Salchichas’)

Tipo de Salchichas

C
A

Para realizar un graáfico de barra de frecuencias absolutas utilizamos la función barplot: como sigue:
y . Veamos ejemplos,
> barplot(table(TIPO), main = ’Tipo de Salchichas’, col = ’blue’)

12
 Tipo de Salchichas

20
15
10
5
0
A B C

Si en lugar del gráfico de barras de frecuencias absolutas estamos interesados en el de frecuencias

relativas, debemos utilizar la función prop.table en el primer argumento de la función barplot como
sigue:
> barplot(prop.table(table(TIPO)), main = ’Tipo de Salchichas’, col = ’blue’)

Tipo de Salchichas
0.30
0.20
0.10
0.00

A B C

8.2. Dos variables cualitativas

Tabla de contigencia: Para realizar una tabla de contingencia utilizamos la función table como sigue:

> table(NumeroAleatorio, Sexo)

Sexo
NumeroAleatorio Female Male
1 1 1
2 5 4
3 11 11
4 15 6
5 10 8
6 13 10
7 28 28
8 11 8
13
 9 8 6
10 1 5

Si en lugar de las frecuencias absolutas queremos calcular las relativas anteponemos a table la función
prop.table especificando con un 1 si queremos calcular las proporciones por fila o con un 2 si queremos
las proporciones por columna (sin esa especificación calcula el porcentaje sobre el total de la tabla).
Veamos un ejemplo:

> prop.table(table(NumeroAleatorio,Sexo),2)

Sexo
NumeroAleatorio Female Male
1 0.009708738 0.011494253
2 0.048543689 0.045977011
3 0.106796117 0.126436782
4 0.145631068 0.068965517
5 0.097087379 0.091954023
6 0.126213592 0.114942529
7 0.271844660 0.321839080
8 0.106796117 0.091954023
9 0.077669903 0.068965517
10 0.009708738 0.057471264

Gráfico de barras: Para realizar un gráfico de barras apilado utilizamos la función barplot:
> a = prop.table(table(Sexo,NumeroAleatorio),1)
> barplot(a)
0.4
0.2
0.0

1 2 3 4 5 6 7 8 9 10

Para realizar un gráfico de barras adosado es necesario agregar la opción beside = TRUE:
> a = prop.table(table(Sexo,NumeroAleatorio),1)
> barplot(a, beside = TRUE)
14
 0.30
0.20
0.10
0.00
1 2 3 4 5 6 7 8 9 10

Observar que para realizar un barplot es importante el orden en que se ponen las variables al armar la
tabla: en segundo lugar debemos poner la variable que queremos que esté en el eje horizontal, teniendo
mucho cuidado al agregar la condición para que las proporciones se calculen por filas o columnas según
corresponda (recordar que deben sumar 1 las proporciones para cada nivel de la variable predictora).

8.3. Una variable cuantitativa

Medidas Resúmenes:

• Medidas de tendencia central

◦ mean(x) devuelve el promedio del argumento x
◦ mean(x,trim = a) devuelve la media recortada del elemento x al a % de cada lado (a ∈
[0, 0,5])
◦ median(x) devuelve la mediana del argumento x
• Medidas de variabilidad
◦ sd(x) devuelve el desvío estandar del argumento x
◦ var(x) devuelve la varianza del vector o matriz x
◦ IQR(x) devuelve el rango intercuartil del argumento x
◦ range(x) devuelve el valor mínimo y el valor máximo en x
◦ quantile(x) devuelve el mínimo, el primer, segundo y tercer quantil, y el máximo de x
◦ quantile(x, a) devuelve el quantil a de x (a ∈ [0, 1])
◦ fivenum(x) devuelve el mínimo, el primer, segundo y tercer quantil, y el máximo de x
◦ summary(x) devuelve el mínimo, el primer quantil, la mediana (o segundo quantil), la media,
el tercer quantil, y el máximo.

En todos los casos es necesario agregar el argumento na.rm = TRUE si x contiene datos perdidos, para
que la función los ignore.

Gráficos:
15
• Histograma: Veamos un ejemplo con la variable CALORIAS:
> par(mfrow = c(2,2))
> hist(CALORIAS, main = ’Histograma’) ] Histograma en escala de frecuencias absolutas
> hist(CALORIAS, main = ’Histograma’, freq = FALSE) ] Histograma en escala de densi-
dad
> hist(CALORIAS, br = seq(80,200,10), main = ’Histograma’) ] Con los puntos de corte
deseados
> hist(CALORIAS, br = c(86,100,125,140,160,180,200),
+ main = ’Histograma’, freq = FALSE)] Con otros puntos de corte

Histograma Histograma
Frequency

0.010
Density
5 10

0.000
0

80 120 160 200 80 120 160 200

CALORIAS CALORIAS

Histograma Histograma
Frequency

0.010
Density
8
4

0.000
0

80 120 160 200 100 140 180

CALORIAS CALORIAS

• Boxplot: Veamos un ejemplo con la variable CALORIAS:

> boxplot(CALORIAS, main = ’Boxplot para la variable CALORIAS’)

16
 Boxplot para la variable CALORIAS

180
160
140
120
100
• Frecuencias acumuladas: Veamos un ejemplo con la variable Velmax:
> plot.ecdf(Velmax)

ecdf(x)
1.0

●●●
●●
●●
●
0.8

●●

●●
●
●●
0.6
Fn(x)

●●●
0.4

●
●●
0.2

●
●●

●
●
●●●
0.0

20 40 60 80 100 120 140 160

8.4. Una variable cuantitativa vs. una cualitativa

Estudiar la relación entre una variable cuantitativa (respuesta) y una cualitativa (predictora) involucra
separar los datos en grupos basados en las “categorías” o “niveles” de la variable predictora.

Medidas Resúmenes: El resumen numérico es dar una medida de centro y de variabilidad de la varia-
ble respuesta por cada nivel de la variable preditora. La función “favstats” (abreviatura de “favorite
statistic”) proporcionas las medidas resúmenes más comunes para variables numéricas (i.e. la medias,
desviación estándar, mediana, etc).
Por ejemplo, si queremos comparar las calorías contenidas en cada tipo de salchichas, podemos
utilizar el siguiente comando:
17
library(mosaic)
favstats(CALORIAS~TIPO)

La respuesta de [R] será:

TIPO min Q1 median Q3 max mean sd n missing

1 A 111 140.5 152.5 177.25 190 156.8500 22.64201 20 0
2 B 107 139.0 153.0 179.00 195 158.7059 25.23580 17 0
3 C 86 102.0 129.0 143.00 170 122.4706 25.48313 17 0

Si sólo nos interesan algunas medidas resúmenes de los que nos proporciona favstats(), podemos
indecarlo de la siguiente manera:

favstats(CALORIAS~TIPO)[c(’TIPO’,’mean’,’sd’)]

En tal caso, [R] nos mostrará:

TIPO mean sd
1 A 156.8500 22.64201
2 B 158.7059 25.23580
3 C 122.4706 25.48313

Gráficos

• Histograma por categoría: Cuando realizamos histogramas para comparar la variable para cada
categoría, es importante que escalas vertical y horizantal sean las mismas para cada gráfico:

> par(mfrow = c(1,3))

> hist(CALORIAS[TIPO == ’A’], br = seq(80,200,20))
> hist(CALORIAS[TIPO == ’B’], br = seq(80,200,20))
> hist(CALORIAS[TIPO == ’C’], br = seq(80,200,20))

18
 Histogram of CALORIAS[TIPO
Histogram
== "A"]
of CALORIAS[TIPO
Histogram
== "B"]
of CALORIAS[TIPO == "C"]

5
4

4
6

3
Frequency

Frequency

Frequency
4

2
2

1
0

0
80 140 200 80 140 200 80 140 200

CALORIAS[TIPO == "A"] CALORIAS[TIPO == "B"] CALORIAS[TIPO == "C"]

• Boxplot por categoría: Veamos un ejemplo:

> boxplot(CALORIAS~TIPO, main = ’Boxplot de CALORIAS para cada TIPO’)

Boxplot de CALORIAS para cada TIPO

180
160
140
120
100

A B C

Otra forma de realizar este gráfico es utilizando la librería lattice:

> bwplot(CALORIAS~TIPO)

También puede observarse la distribución de la variable cuantitativa por categoría de la cualitativa

graficando las densidades como sigue:

> densityplot(~CALORIAS, groups = TIPO, auto.key = TRUE)

19
 A
B
C

0.015

0.010

Density
0.005

0.000

50 100 150 200

CALORIAS

• Diagrama de dispersión: Veamos un ejemplo donde comparamos las velocidades máximas (Vel-
Max) para cada sexo (Sexo). La función as.numeric transforma la variable Sexo en un vector de
1 y 2 (pues son dos categorías) asignandole 1 a Female y 2 a Male.
> plot(as.numeric(Sexo), Velmax, xlab = ’Sexo’, xlim = c(0.5,2.5))

●
●
100 120 140

● ●
●
●
● ●
●
● ●
●
●
● ●
● ●
●
Velmax

● ●
● ●
●
● ●
●
●
●
●
● ●
80

●
● ●
●
●
60

●
●
●
40

0.5 1.0 1.5 2.0 2.5

Sexo

Dos o más gráficos superpuestos: Para superponer gráficos es necesario agregar la opción add = TRUE al
partir del segundo gráfico. Por ejemplo, si queremos superponer dos histogramas hacemos lo siguiente:

> hist(CALORIAS[TIPO == ’A’], br = seq(80,200,20), main = ’Histograma

+ de CALORIAS’, xlab = ’’, ylab = ’’)
> hist(CALORIAS[TIPO == ’B’], br = seq(80,200,20), add = TRUE,
+col = rgb(1, 0, 0, 1/5))
> hist(CALORIAS[TIPO == ’C’], br = seq(80,200,20), add = TRUE,
+col = rgb(0, 1, 0, 1/5))

20
 Histograma de CALORIAS

8
6
4
2
0

80 100 120 140 160 180 200

8.5. Dos variables cuantitativas

Podemos realizar un gráfico de dispersión para observar la relación que existe entre dos variables cuantitativas.
Existen varias formas de realizar esto, la más común es utilizar la función plot, pero existen otras alternativas
que se muestran a continuación:

plot(CALORIAS, SODIO)
600
500
SODIO

400
300
200

100 120 140 160 180

CALORIAS

Si queremos discriminar por categorías, cambiar colores, etc., se puede ver en la Ayuda de RStudio cómo
hacerlo
Otra forma más directa de hacerlo es utilizando el paquete car:
21
library(car)
scatterplot(Calorias, Sodio, groups = Tipo, legend.title = ’Tipo’,
legend.coords = ’bottomright’, smooth = FALSE, reg.line = FALSE, boxplots = FALSE)

600
500
SODIO

400
300

Tipo
200

A
B
C

100 120 140 160 180

CALORIAS

8.6. Una variable cuantitativa vs. dos cualitativas (factores)

Un resúmen numérico de la variable cuantitativa para las diferentes combinaciones de categorías (niveles) de
las dos variables cualitativas (factores) se puede obtener de la siguiente manera:

> datos = read.table(’Contraccion.txt’, header = TRUE)

> attach(datos)

> library(mosaic)
> favstats(Contraccion~Velocidad+Temperatura)

min Q1 median Q3 max mean sd n missing

Alta.Alta 91.62 93.0325 93.270 93.4450 93.89 93.168 0.62668085 10 0
Baja.Alta 75.82 76.0075 76.135 76.1775 76.25 76.090 0.13507200 10 0
Alta.Bajo 71.48 71.5275 71.590 71.6950 71.74 71.606 0.09605554 10 0
Baja.Bajo 72.42 72.5575 72.600 72.8000 73.07 72.674 0.20630075 10 0

Para realizar el gráfico de las medias de una variable cuantitativa para dos combinaciones de factores y así
poder visualizar si existen interacciones entre los diferentes factores hacemos

> interaction.plot(Velocidad, Temperatura, Contraccion)

22
 Temperatura

Alta

90
Bajo

mean of Contraccion

85
80
75

Alta Baja

Velocidad

8.7. Más sobre gráficos

Para realizar más de un gráfico en una misma página, antes de graficar debemos agregar la opción
par(mfrow = c(#1,#2)) con #1 la cantidad de gráficos que queremos por fila y #2 la cantidad que
queremos por columna.

Para realizar un nuevo gráfico sin borrar el anterior, antes del nuevo gráfico debemos agregar la opción
dev.new().

Para agregar una leyenda al gráfico, dentro de la función gráfica (plot, barplot, etc.) debemos agregar la
opción legend indicando en un vector los nombres de las variables. Por ejemplo, legend = c(’Hombres’,
’Mujeres’).

Para agregar/modificar las etiquetas de los ejes del gráfico, dentro de la función gráfica (plot, barplot,
etc.) debemos agregar la opción xlab o ylab con el nombre que deseamos ponerle al eje. Por ejemplo,
xlab = ’Peso en kg.’, ylab = ’Altura en cm.’.

Para modificar los ejes del gráfico, dentro de la función gráfica (plot, barplot, etc.) debemos agregar
la opción xlim o ylim indicando en un vector los limites queremos que tengan los ejes. Por ejemplo,
xlim = c(0,100) o ylim = c(0,200).

Si deseamos cambiar las marcas de los ejes, primero debemos borrar las que hace por defecto agre-
gando la opción xaxt = ’n’ o yaxt = ’n’ dentro de la función gráfica y luego agregando la sentencia
axis(side = eje, at = marcas). Por ejemplo, si queremos poner marcas en 1, 2 y 3 del eje x, de-
bemos agregar axis(side = 1, at = c(1,2,3)).

Para cambiar el tipo de gráfico, dentro de la función gráfica (plot, barplot, etc.) debemos agregar la
opción type indicando el tipo que gráfico que deseamos realizar. Por ejemplo, si deseamos que grafique

23
 líneas, debemos escribir type = ’l’, si deseamos que grafique puntos, debemos escribir type = ’p’,
etc. Recurrir a la ayuda de R para ver más tipos de gráficos (help(plot)).

Para cambiar el caracter del gráfico, dentro de la función gráfica (plot, barplot, etc.) debemos agregar la
opción pch indicando el tipo de caracter que deseamos realizar. Por ejemplo, si deseamos que grafique
cuadraditos, debemos escribir pch = 0, si deseamos círculos, debemos escribir pch = 1, triángulos
debemos escribir pch = 2, etc. Recurrir a la ayuda de R para ver más tipos de gráficos (help(points)).

Para cambiar el color del gráfico, dentro de la función gráfica (plot, barplot, etc.) debemos agregar la
opción col indicando el color que deseamos. Por ejemplo, si deseamos que el gráfico sea rojo, debemos
escribir col = ’red’ o col = 2. Recurrir a la ayuda de R para ver más tipos de gráficos (help(plot)).
Para poner colores personalizados, podemos utilizar la opción col = rgb(#1, #2, #3, #4) con #1
la cantidad de rojo que deseamos, #2 la cantidad de verde y #3 la cantidad de azul. El argumen-
to #4 es para especificar la transparencia de los colores. Por ejemplo, para graficar en rojo escribimos
col = rgb(1, 0, 0, 1), para graficar en rojo pero más transparente escribimos col = rgb(1, 0, 0, 1/5).

Ejemplo:

par(mfrow = c(1,2))
plot(SODIO, CALORIAS, main = ’Antes’)

plot(SODIO, CALORIAS, col = rgb(0.6, 0, 0.5, 1), pch = 6, xaxt = ’n’,

+ yaxt = ’n’, xlim = c(100,700), ylim = c(0,200), xlab = ’Sodio’,
+ ylab = ’Calorias’, main = ’Despues’)
axis(side = 1, at = c(100,300,500,700))
axis(side = 2, at = c(0,100,200))

Antes Despues
200
180
160
CALORIAS

Calorias

100
140
120
100

200 400 600 100 300 500 700

SODIO Sodio

24
9. Funciones de distribución
9.1. Variables aleatorias discretas
En esta sección veremos como graficar funciones de probabilidad y funciones de distribución acumulada tanto
empíricas como teóricas.

1. Comencemos con un ejemplo de una variable aleatoria discreta: sea X una variable aleatoria que
toma valores 1, 3 y 4 con probabilidades 0.5, 0.2 y 0.3, respectivamente. Definimos los elementos que
necesitaremos:

x = c(1, 3, 4)
p = c(0.5, 0.2, 0.3)

a) Distribución teórica: Para graficar la función de probabilidad puntual de X utilizamos la opción

type = ’h’ dentro de la función plot como sigue:

plot(x, p, type = ’h’, xaxt = "n", ylim = c(0,max(p)),

+ main = ’Funcion de probabilidad puntual’)
points(x, p , pch = 16)
axis(side = 1, at = c(1,3,4))

Para calcular la función de distribución acumulada utilizamos la función stepfun y para graficarla
utilizamos la función plot.stepfun como sigue:

f = stepfun(x, c(0, cumsum(p)))

plot.stepfun(f, vertical = FALSE, pch = 16,
+ ylab = ’F(X)’, main = ’Funcion de distribucion acumulada’)

La función cumsum(x) devuelve un vector con las sumas acumuladas de los elementos de x.

Funcion de probabilidad puntual Funcion de distribucion acumulada

0.5

1.0
0.4

0.8
0.3

0.6
F(X)
p

0.2

0.4
0.1

0.2
0.0

0.0

1 3 4 0 1 2 3 4 5

x x

25
 b) Distribución empírica: Para obtener una muestra de tamaño 1000 de la variable aleatoria X,
utilizamos la función sample como sigue:

n = 1000
muestra = sample(x, prob = p, size = n, replace = TRUE)

Para graficar el histograma con estos datos utilizamos la función hist (ver sección 8.3).

hist(muestra, br = seq(0.5, 4.5), prob = TRUE, xlab = ’x’,

+ ylab = ’frecuencias relativas’, main = ’Histograma muestral’ )

Para graficar la función de distribución acumulada empírica utilizamos la función plot.ecdf

como sigue:

plot.ecdf(muestra, main = ’Funcion de distribucion acumulada de X’)

Histograma muestral Funcion de distribucion acumulada de X

0.5

1.0
0.4

0.8
frecuencias relativas

0.3

0.6
Fn(x)
0.2

0.4
0.1

0.2
0.0

0.0

1 2 3 4 0 1 2 3 4 5

x x

Observación: La distribución uniforme discreta se realiza de igual manera simplemente cambiando

el vector de valores y de probabilidades como sigue:

k = 10
x = 1:10
p = rep(1/k, k)

2. Distribución Binomial: sea X una variable aleatoria con distribución Binomial de parámetros m y p.
Definimos los elementos que necesitaremos:

m = 6
p = 0.5
x = 0:m

26
a) Distribución teórica: Para obtener la función de probabilidad puntual y la función de distri-
bución acumulada de X utilizamos, respectivamente, las funciones dbinom(x, size, prob) y
pbinom(x, size, prob) con argumentos x = 0:m (posibles valores de la variable aleatoria X),
size = m (el número de ensayos Bernoulli) y prob = p (la probabilidad de éxito).

par(mfrow = c(1,2))

plot(x, dbinom(x,m,p), type = ’h’, xaxt = "n", ylim = c(0,1),

+ ylab = ’frecuencias relativas’,
+ main = ’Funcion de probabilidad puntual’)
points(x, dbinom(x,m,p) , pch = 16)
axis(side = 1, at = x)

f = stepfun(x, c(0, pbinom(x,m,p)))

plot.stepfun(f, vertical = FALSE, pch = 16, ylab = ’F(X)’,
+ main = ’Funcion de distribucion acumulada’)

Funcion de probabilidad puntual Funcion de distribucion acumulada

1.0

1.0
0.8

0.8
frecuencias relativas

0.6

0.6
F(X)
0.4

0.4
0.2

0.2
0.0

0.0

0 1 2 3 4 5 6 0 2 4 6

x x

b) Distribución empírica: Para obtener una muestra de tamaño 1000 de la variable aleatoria X,
utilizamos la función rbinom como sigue:

n = 1000
muestra = rbinom(n, m, p)

Para graficar el histograma y la función de distribución acumulada empírica de X lo hacemos

como antes:

par(mfrow = c(1,2))

hist(muestra, br = seq(-0.5,m+.5), xaxt = "n",

+ ylim = c(0,1), prob = TRUE, xlab = ’x’, ylab = ’frecuencias relativas’,
27
 + main = ’Histograma muestral’)
axis(side = 1, at = x)

plot.ecdf(muestra, xaxt = "n",

+ main = ’Funcion de distribucion acumulada de X’)
axis(side = 1, at = x)

Histograma muestral Funcion de distribucion acumulada de X

1.0

1.0
0.8

0.8
frecuencias relativas

0.6

0.6
Fn(x)
0.4

0.4
0.2

0.2
0.0

0.0
0 1 2 3 4 5 6 0 1 2 3 4 5 6

x x

Observaciones:

El gráfico de probabilidad puntual (teórico) de X

plot(x, p, type = ’h’)

points(x, p , pch = 16)

puede reemplazarse por un histograma: barplot(prop.table(table(x))).

El histograma muestral hist(muestra) puede reemplazarse por:

barplot(prop.table(table(muestra)))

En general, cuando trabajamos con una variable aleatoria discreta (excepto para la uniforme que difiere
levemente) realizamos los gráficos anteriores de manera análoga, teniendo en cuenta la siguiente tabla:

Generación de Función de Función de

Distribución números aleatorios probabilidad puntual distribución acumulada Cuantiles
Bernoulli rbinom dbinom pbinom qbinom
Binomial rbinom dbinom pbinom qbinom
Poisson rpois dpois ppois qpois
Hipergeométrica rhyper dhyper phyper qhyper

28
9.2. Variables aleatorias continuas
En esta sección veremos como graficar funciones de densidad y funciones de distribución acumulada tanto
empíricas como teóricas.

1. Distribución Normal: sea X una variable aleatoria con distribución normal de parámetro µ y σ. Defi-
nimos los elementos que necesitaremos:

mu = 2
sigma = 2

a) Distribucion teórica: Para obtener la función de densidad y la función de distribución acumulada

de X utilizamos, las funciones dnorm(x,mu,sigma) y pnorm(x,mu,sigma), respectivamente. Para
graficarlas utilizamos el comando curve como sigue:

par(mfrow = c(1,2))

x1 = mu-3*sigma
x2 = mu+3*sigma
curve(dnorm(x, mu, sigma), xlim = c(x1,x2), main = ’Funcion de densidad’,
+ylab = ’f(x)’)

curve(pnorm(x, mu, sigma), xlim = c(-10,15), ylab = ’F(X)’,

+main = ’Funcion de distribucion acumulada’)

Funcion de densidad Funcion de distribucion acumulada

0.20

0.0 0.2 0.4 0.6 0.8 1.0

0.15

F(X)
0.10
f(x)

0.05
0.00

-4 -2 0 2 4 6 8 -10 -5 0 5 10 15

x x

b) Distribución empírica: Para obtener una muestra de tamaño 1000 de la variable aleatoria X,
utilizamos la función rnorm de la siguiente forma:

n = 1000
muestra = rnorm(n, mu, sigma)

29
 Para graficar el histograma y la función de distribución acumulada empírica de X lo hacemos
como sigue:

par(mfrow = c(1,2))

hist(muestra, prob = TRUE, ylab = ’frecuencias relativas’,

+main = ’Histograma muestral’)

plot.ecdf(muestra, xaxt = "n", main = ’Funcion de distribucion

+acumulada de X’)

Histograma muestral Funcion de distribucion acumulada de X

1.0
0.20

0.8
frecuencias relativas

0.15

0.6
Fn(X)
0.10

0.4
0.05

0.2
0.00

0.0
-5 0 5 -5 0 5 10

x x

En general, cuando trabajamos con una variable aleatoria continua realizamos los gráficos anteriores de ma-
nera análoga, teniendo en cuenta la siguiente tabla:

Generación de Función de Función de

Distribución números aleatorios probabilidad puntual distribución acumulada Cuantiles
Normal rnorm dnorm pnorm qnorm
Uniforme runif dunif punif qunif
Chi-cuadrada rchisq dchisq pchisq qchisq
F rf df pf qf
T-student rt dt pt qt
Lognormal rlnorm dlnorm plnorm qlnorm
Exponencial rexp dexp pexp qexp

10. Intervalos de confianza

En esta sección encontrarán un resumen de como costruir en R todos los intervalos de confianza que vimos
en clase.

30
10.1. Una muestra
Intervalo de confianza para µ con varianza conocida: No podemos hacerlo con la compu, lo debemos
calcular a mano.

Intervalo de confinza para µ con varianza desconocida: Vamos a utilizar la función de R t.test de la
siguiente forma:

> t.test(x, alternative = "two.sided", conf.level = 0.95)

• El vector x contiene la muestra.

• En la opción alternative indicamos si queremos un intervalo a bilateral ("two.sided"), unilateral
superior ("less") o unilateral inferior ("greater").
• En conf.level indicamos el nivel de confianza con el que queremos construir el intervalo.

10.2. Dos muestras

Intervalo de confianza para µ1 −µ2 con varianza conocidas: No podemos realizarlo en la compu, tenemos
que calcularlo a mano.

Intervalo de confianza para µ1 − µ2 con varianza desconocidas, pero iguales: En la función t.test
ingresamos lo siguiente

> t.test(x, y, alternative = "two.sided", conf.level = 0.95, var.equal = T)

• En los vectores x e y ingresamos las dos muestras

• En var.equal indicamos que queremos que el intervalo se construya considerando varianzas
iguales.
• Los demás argumentos son como en el caso de una muestra.

Intervalo de confianza para µ1 − µ2 con varianza desconocidas, pero distintas: Realizamos lo mismo
que en el item anterior, pero ingresando en var.equal el valor lógico False:

> t.test(x, y, alternative = "greater", conf.level = 0.99, var.equal = F)

Por ejemplo, en este caso construimos el intervalo unilateral inferior con 99 % de confianza.

Intervalo de confianza para muestras apareadas: Aclaramos esta condición en t.test por medio del
argumento paired = T

> t.test(x, y, alternative = "two.sided", conf.level = 0.95, paired = T)

σ12
Intervalo de confianza para : Lo realizamos mediante la función var.test
σ22

> var.test(x, y, alternative = "two.sided", conf.level = 0.95)

31
10.3. Intervalos de confianza para una proporción
Una muestra: Vamos a utilizar la funcion prop.test de la siguiente forma

> prop.test(X, n, conf.level = 0.95, alternative = "two.sided")

• X es el número de éxitos en n experimentos.

• Los demás argumetos son análogos a lo visto antes.

Dos muestras: Para calcular el intervalo de confianza para la diferencia de proporciones utilizamos
prop.test como se explica a continuación

> prop.test(x, n, conf.level = 0.95, alternative = "two.sided")

• x es un vector de longitud 2, que contiene la cantidad de éxitos para cada muestra.

• n es un vector de longitud 2 que contiene el número de experimentos para cada muestra
• Los demás argumetos son análogos a lo visto antes.

11. Test de hipótesis

Para obtener el p-valor correspondiente a un test de hipótesis en general usamos las mismas funciones que
para la construcción de intervalos de confianza. Supongamos un test para la media de una población en el
que queremos testear la hipótesis nula H0 : µ = 500 contra la altenativa H1 : µ > 500, con un nivel de sig-
nificancia α = 0,01 y disponemos de una muestra de datos guardada en el vector muestra. En este caso usamos

> t.test(muestra,alternative="greater",mu=500,conf.level=0.99)

Si, en cambio, queremos testear la hipótesis nula H0 : µ = 500 contra la altenativa H1 : µ < 500, con un
nivel de significancia α = 0,05, en este caso usamos

> t.test(muestra,alternative="less",mu=500,conf.level=0.95)

Por último, si la hipótesis alternativa es en cambio H1 : µ 6= 500, la orden a utilizar es

> t.test(muestra,alternative="two.sided",mu=500,conf.level=0.95)

R también proporciona funciones para el cálculo directo de la potencia de un t-test, a través de la función
power.t.test(). Para su uso debemos recordar que la potencia (power) de un t-test es función del tamaño
muestral (n), del nivel de significancia (sig.level), del desvío estándar de la población de la cual se obtuvo
la muestra (sd), de la diferencia entre la media verdadera y la postulada por la hipótesis nula (delta) y del
tipo de hipótesis alternativa (alternative); es decir, si se trata de un test unilateral o bilateral. También

32
debemos precisar si se trata de un test para una sola muestra o para dos muestras. La función power.t.test
espera que especifiquemos todos excepto uno de esos datos y con ello computa el aspecto restante. Por ejem-
plo, si nos interesa calcular la potencia del test para la media dado por la hipótesis nula H0 : µ = 500 contra
la altenativa H1 : µ < 500, con un nivel de significancia α = 0,01, usando una muestra de tamaño n = 20,
cuando la media verdadera es µ = 480 y la varianza de la población es σ 2 = 202 , usamos

> power.t.test(type="one.sample",n=20,delta=20,sig.level=0.01,sd=20,alternative="one.sided")

Obtenemos

One-sample t test power calculation

n = 20
delta = 20
sd = 20
sig.level = 0.01
power = 0.9653099
alternative = one.sided

de donde vemos que la potencia del test especificado es 0,9653099. Si para las mismas hipótesis quisiéramos
saber cuál es el n que deberíamos usar para tener una potencia de 0,90, entonces usamos

> power.t.test(power=0.9,delta=20,sig.level=0.01,sd=20,alternative="one.sided")

En este caso obtenemos

One-sample t test power calculation

n = 15.85579
delta = 20
sd = 20
sig.level = 0.01
power = 0.9
alternative = one.sided

El resultado nos muestra que el tamaño de muestra más chico que nos asegura esa potencia es n = 16.

Cuando tenemos un test para la diferencia de las medias de dos poblaciones, también usamos t.test() y
especificamos si consideraremos que las varianzas (poblacionales) son iguales o no del mismo modo que lo
hicimos cuando discutimos la construcción de intervalos de confianza. Para los cálculos de potencia utilizan-
do la función power.t.test() es necesario especificar que se trata de un test para dos muestras a través
de la opción type="two.sample". Es necesario notar también que esta función solamente realiza

33
cómputos relacionados con la potencia para un test en el que las dos muestras tienen el mismo
tamaño y se supone además que la varianza de ambas poblaciones es la misma.

Un último comentario sobre test para dos muestras refiere al uso de formas del tipo

> t.test(CALIFICACIONA~SEXO,alternative="greater")

En este caso estamos queriendo testear si las calificaciones (medias) obtenidas en Matemática A por las mu-
jeres son superiores a las obtenidas por los varones. La orden suministrada a R realiza ese test, en este caso
suponiendo que las varianzas poblacionales son iguales (lo que no implica que sea la opción más adecuada
en este caso). Sin embargo, no es obvio que el test realizado realmente emplea la hipótesis alternativa que
queremos. ?ťCómo nos damos cuenta que la alternativa usada es que la media de las mujeres es mayor que
la de los hombres y no al revés? La respuesta es que el nivel femenino del factor SEXO se considera antes
que el nivel masculino, por el orden alfabético de sus letras iniciales. Para evitar confusiones, se sugiere usar
una forma explícita del tipo

> t.test(CALIFICACIONA[SEXO=="femenino"],CALIFICACIONA[SEXO=="masculino"],
alternative="greater")

12. ANOVA
Para realizar un anális ANOVA en R vamos a necesitar instalar los siguientes paquetes: car, multcomp y
lawstat.

PASO 1: Resumen numérico y gráfico de los datos.

Mediante un gráfico boxplot por niveles del factor podemos elaborar algunas conjeturas.

> boxplot(Respuesta~Factor)

Podemos realizar un resumen numérico por niveles del factor de la siguiente forma:

> ns = tapply(Respuesta, factor, length)

> Promedio = tapply(Respuestas ,factor, mean)
> Desvio = tapply(Respuestas, factor, sd)
> print(cbind(ns, Promedio, Desvio), digits = 3)

PASO 2: Realizar el test ANOVA.

Es muy importante primero declarar la variable Factor de la siguiente manera:

34
> factor = as.factor(factor)

Por medio de la función aov realizamos la tabla anova

> modelo = aov(Respuesta ~ factor)

> summary(modelo)

Df Sum Sq Mean Sq F value Pr(> F)

factor 5 12734 2546.8 57.10 <2.2e-16
Residuals 343 15297 44.6
---
El valor que nos indica si rechazamos o no la hipótesis nula (p-value) esta marcado en color rojo.

PASO 3: Si en el paso 2 rechazamos H0 realizamos comparaciones múltiples. En caso contrario, pasar

al paso siguiente.

Mediante la función glht de R, aplicamos los métodos que estudiamos en clase

• Tukey

> T = glht(modelo, linfct = mcp(Factor = "Tukey"))

> summary(T)

La tabla mostrará todas las comparaciones posibles entre los niveles del factor, y en la última
columna los p-valores. Por ejemplo, para los datos serum.txt se obtuvo lo siguiente

Linear Hypotheses:
Estimate Std. Error t value Pr(> |t|)
G2 - G1 == 0 -24.084 14.483 -1.663 0.3522
G3 - G1 == 0 20.593 14.483 1.422 0.4891
G4 - G1 == 0 28.046 15.845 1.770 0.2986
G3 - G2 == 0 44.677 17.009 2.627 0.0552 .
G4 - G2 == 0 52.130 18.183 2.867 0.0312 *
G4 - G3 == 0 7.453 18.183 0.410 0.9761

Si queremos obtener lo intervalos de confianza, hacemos:

confint(T)

• Bonferroni
Utilizamos lo que ya guardamos en T como sigue:

summary(T, test = adjusted("bonferroni"))

35
 Nos devolverá una tabla similar a la anterior pero con los p-valores correspondientes a este
método.
Si ahora queremos los intervalos de confianza, hacemos:

k = 6
alfa = 0.05
tt = qt(1-alfa/(2*k))
confint(T, calpha = tt)

Observación: El valor de k especificado en el comando anterior, corresponde a la cantidad de

comparaciones múltiples que se realizan.
• Dunnet: Este método lo usamos cuando queremos comparar los niveles del factor con un control.
Es importante reordenar los niveles del factor, de modo que el primero corresponda al nivel
control.

> factor1 = relevel(Factor,’control’)

> modelo1 = aov(Respuesta ËIJ factor1)
> dunnett = glht(modelo1, linfct = mcp(factor1 = ’Dunnett’))
> summary(dunnett)

PASO 4: Análisis de los residuos.

Contamos con los siguientes test y gráficos para analizar los residuos:

• Gráfico de residuos vs valores ajustados y gráfico de normalidad de los residuos.

Son los dos primeros gráficos que se obtienen mediante el comando > plot(modelo).

• Test de igualdad de varianzas.

> levene.test(Respuesta, Factor, location = ’median’)

Recordad que este test no es robusto a la no normalidad.

• Test de presencia de outliers.

> outlierTest(modelo)
rstudent unadjusted p-value Bonferonni p
21 -3.71598 0.00023648 0.002294

• Gráfico de la distancia de Cook y outliers

> influenceIndexPlot(modelo)

36
 PASO 5: Si en el paso 4 no encontramos problemas con los residuos, el análisis terminó. En caso con-
trario, realizamos las modificaciones o transformacines necesarias y retomamos al PASO 1.

Si queremos quitar un dato que resultó atípico según el test de outliers, como por ejemplo el dato 21
en el item anterior, recordar que lo hacemos definiendo nuevas variables sin ese dato:

> Respuesta2 = Respuesta[-21]

> factor2 = factor[-21]

Un gráfico que es útil para analizar si es necesario una tranformación es el siguiente:

> boxcox(modelo)

13. Regresión
En lo que sigue llamaremos Y donde guardamos los datos de la variable respuesta y X para los datos predic-
tores.

Un diagrama de dispersión de las variables lo realizamos mediante

> plot(X,Y)

Este gráfico es muy importante para observar como se relacionan los datos y elaborar hipótesis. Para ajustar
el modelo de regresión lineal utilizamos la función lm.

> modelo = lm(Y ~ X)

> summary(modelo)

Estimate Std. Error t value Pr(> |t|)

(Intercept) 1.769135 0.078664 22.490 <2e-16 ***

X -0.001965 0.003493 -0.563 0.576

Residual standard error: 0.09258 on 49 degrees

of freedom Multiple R-squared: 0.006416,
Adjusted R-squared: -0.01386
F-statistic: 0.3164 on 1 and 49 DF, p-value: 0.5763

Veamos cual es la información que me da esta tabla:

37
 Los estimadores de β0 (ordenada al origen) y β1 (pendiente) en la columna Estimate. En este ejemplo
βb0 = 1,769 y βb1 = −0,001

Un estimador de σ, que en este caso es σ

b = 0,092

En la última linea del summary, el F-test compara: H0: modelo con solamente la media vs H1: modelo
más completo que estoy ajustando. En este caso, el p-value da 0.57 por lo que aceptamos H0.

Si queremos agregrar al gráfico de dispersión la recta ajustada, lo hacemos como sigue:

> plot(X, Y)
> lines(X, modelo$fitted)

También podemos hacer:

> plot(X, Y)
> abline(modelo)

Como en ANOVA, tenemos que estudiar los residuos. Las funciones que utilizaremos son similares:

> plot(modelo)
> influenceIndexPlot(modelo, id.n = 3)
> outlierTest(modelo)

El outlierTest devuelve el p-valor correspondiente al test H0 : No existen outliers vs. H1 : Existen outilies.
Además, en el caso de existir outliers nos dice en que posición se encuentran los mismos.
Para obtener los intervalos de confianza utilizamos la función predic. Si queremos intervalo de confianza
para el valor medio lo aclaramos en el argumento:

> predict(modelo, interval = "confidence")

En cambio, si queremos un intervalo de predicción:

> predict(modelo, interval = "predict")

Estas funciones nos devuelve 3 columnas. En la primera se da el valor ajustado, en la segunda el límite inferior
y luego el límite superior del intervalo para cada valor de la variable predictora que está en el conjunto de
datos. Si en cambio, nos interesa un intervalo para un valor en particular a (que puede no estar en el conjunto
de datos X) procedemos como sigue:

> nuevox = data.frame(X = a)

> predict(modelo, nuevox, interval = "confidence")
> predict(modelo, nuevox, interval = "predict")

En el caso de existir outliers y querer eliminarlos podemos proceder de varias maneras, una de ellas es
eliminar a mano los datos atípicos (como se mostró en la Sección anterior). Supongamos que el outlierTest
nos dice que los datos atípicos se encuentran en las posiciones 2 y 23 de los datos. Por lo tanto hacemos:
38
> Xnuevo = X[-c(2,23)]
> Ynuevo = Y[-c(2,23)]

y volvemos a ajustar el modelo:

> modelo1 = lm(Ynuevo ~ Xnuevo)

Otra forma es utilizar directamente la sentencia subset dentro de la función lm como sigue:

> modelo1 = lm(Y ~ X, subset = -c(2,23))

En caso de realizar una transformación, recordemos de la Sección anterior que podemos ver que transforma-
ción aplicarle a la variable respuesta mediante boxcox(modelo). Por ejemplo, si éste gráfico resulta como el
que sigue:

95%
-500
log-Likelihood

-700
-900

-2 -1 0 1 2

Al estar el intervalo de confianza centrado en cero, esto nos indicaría que debemos aplicar una transformación
logarítmica. Por el contrario, el índice indicará la potencia a la cual deberá ser elevada la variable respuesta
para ser transformada. En este caso procedemos de la siguiente manera:

> modelo = lm(log(Y) ~ X)

Y si, por ejemplo, queremos obtener intervalos de confianza y predicción para un valor en particular a, en
este caso debemos proceder de la siguiente manera:

> nuevox = data.frame(X = a)

> exp(predict(modelo, nuevox, interval = "confidence"))
> exp(predict(modelo, nuevox, interval = "predict"))

Es decir, debemos aplicar la función inversa del logaritmo (exponencial) para obtener los intervalos para los
valores originales.

39