Matriz de correlación y valores p con Hmisc en R

Para examinar la asociación entre variables en nuestros datos vamos a ver cómo calcular una matriz de correlación.

Datos

Como ejemplo, usamos el conjunto de datos swiss sobre la fertilidad y datos socioeconómicos de 47 provincias Suiza.

?swiss

head(swiss)

             Fertility Agriculture Examination Education Catholic
Courtelary        80.2        17.0          15        12     9.96
Delemont          83.1        45.1           6         9    84.84
Franches-Mnt      92.5        39.7           5         5    93.40
Moutier           85.8        36.5          12         7    33.77
Neuveville        76.9        43.5          17        15     5.16
Porrentruy        76.1        35.3           9         7    90.57
             Infant.Mortality
Courtelary               22.2
Delemont                 22.2
Franches-Mnt             20.2
Moutier                  20.3
Neuveville               20.6
Porrentruy               26.6

Matriz de correlación

Representamos la matriz de correlación entre las 6 variables.

cor(swiss)

                  Fertility Agriculture Examination   Education   Catholic
Fertility         1.0000000  0.35307918  -0.6458827 -0.66378886  0.4636847
Agriculture       0.3530792  1.00000000  -0.6865422 -0.63952252  0.4010951
Examination      -0.6458827 -0.68654221   1.0000000  0.69841530 -0.5727418
Education        -0.6637889 -0.63952252   0.6984153  1.00000000 -0.1538589
Catholic          0.4636847  0.40109505  -0.5727418 -0.15385892  1.0000000
Infant.Mortality  0.4165560 -0.06085861  -0.1140216 -0.09932185  0.1754959
                 Infant.Mortality
Fertility              0.41655603
Agriculture           -0.06085861
Examination           -0.11402160
Education             -0.09932185
Catholic               0.17549591
Infant.Mortality       1.00000000

Redondeamos los resultados a dos decimales, para despejar la presentación. Se ve con claridad como la diagonal principal son unos que dividen a la matriz en dos partes simétricas.

round(cor(swiss), 2)

            Fertility Agriculture Examination Education Catholic Infant.Mortality
Fertility           1.00      0.35       -0.65     -0.66     0.46           0.42
Agriculture         0.35      1.00       -0.69     -0.64     0.40          -0.06
Examination        -0.65     -0.69        1.00      0.70    -0.57          -0.11
Education          -0.66     -0.64        0.70      1.00    -0.15          -0.10
Catholic            0.46      0.40       -0.57     -0.15     1.00           0.18
Infant.Mortality    0.42     -0.06       -0.11     -0.10     0.18           1.00

Se puede observar como la fertilidad está negativamente asociada con Education y Examination. En cambio, está positivamente asociada con Agriculture, Catholic y Infant Mortality.

Contraste de hipótesis para un par de variables

cor.test(swiss$Fertility, swiss$Education)

 Pearson's product-moment correlation

data:  swiss$Fertility and swiss$Education
t = -5.9536, df = 45, p-value = 3.659e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.7987075 -0.4653206
sample estimates:
       cor 
-0.6637889

La salida de datos incluye la t de Student del contraste de hipótesis, los grados de libertad (df), el valor p (probabilidad), los intervalos de confianza al 95% y el coeficiente de correlación.

Valores p entre todas las variables

Necesitamos emplear el paquete Hmisc. Es necesario convertir previamente el data frame en una matriz.

library(Hmisc)
rcorr(as.matrix(swiss))

            Fertility Agriculture Examination Education Catholic Infant.Mortality
Fertility           1.00      0.35       -0.65     -0.66     0.46           0.42
Agriculture         0.35      1.00       -0.69     -0.64     0.40          -0.06
Examination        -0.65     -0.69        1.00      0.70    -0.57          -0.11
Education          -0.66     -0.64        0.70      1.00    -0.15          -0.10
Catholic            0.46      0.40       -0.57     -0.15     1.00           0.18
Infant.Mortality    0.42     -0.06       -0.11     -0.10     0.18           1.00

n= 47 


P
            Fertility Agriculture Examination Education Catholic Infant.Mortality
Fertility                  0.0149      0.0000      0.0000    0.0010   0.0036          
Agriculture      0.0149                0.0000      0.0000    0.0052   0.6845          
Examination      0.0000    0.0000                  0.0000    0.0000   0.4454          
Education        0.0000    0.0000      0.0000                0.3018   0.5065          
Catholic         0.0010    0.0052      0.0000      0.3018             0.2380          
Infant.Mortality 0.0036    0.6845      0.4454      0.5065    0.2380         

La salida de datos incluye la matriz de correlaciones (automáticamente redondea a dos decimales) y la matriz de valores p. Se puede observar que casi todos los valores p son estadísticamente significativos salvo algunos relacionados con Infant.Mortality, pero el coeficiente de correlación era casi cero en cualquier caso.

Entradas relacionadas

• Regresión lineal simple para principiantes en R
• Regresión lineal simple
• Regresión lineal simple en R con ggplot2
• Regresión lineal simple para principiantes en Excel

Referencias

• Hmisc: Harrell Miscellaneous
• Coeficiente de correlación de Pearson
• Matriz de correlación
• Valor p
• Prueba t de Student

Dividir un data frame en los grupos definidos por una columna en R

Problema

Queremos dividir un data frame por una de las columnas, creando tantos data frames como grupos contiene la columna. Por ejemplo, dividiremos el data frame ChickWeight por la columna Chick que identifica 50 tipos de polluelos.

str(ChickWeight$Chick)

Ord.factor w/ 50 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 15 15 15 15 ...

Solución

Crear una lista de data frames.

ldf <- split(ChickWeight, ChickWeight$Chick)

Para acceder a un data frame de la lista de 50, por ejemplo el 20:

ldf$`20`

    weight Time Chick Diet
209     41    0    20    1
210     47    2    20    1
211     54    4    20    1
212     58    6    20    1
213     65    8    20    1
214     73   10    20    1
215     77   12    20    1
216     89   14    20    1
217     98   16    20    1
218    107   18    20    1
219    115   20    20    1
220    117   21    20    1

Sin embargo, para prescindir de los números como nombres de data frames, creamos una nueva columna en la que anteponemos a cada número la cadena de texto Chick_.

ChickWeight$id <- paste0("Chick_", ChickWeight$Chick)
# Creamos la lista de nuevo
ldf <- split(ChickWeight, ChickWeight$id)
# Accedemos a un data frame
ldf$Chick_20

Crear data frames independientes

Si queremos crear un objeto separado para cada data frame en nuestro entorno de trabajo.

list2env(ldf, .GlobalEnv)

En nuestro entorno habremos creado 50 data frames.

Referencias

• stackoverflow

Diagramas de dispersión con el paquete car en R

Introducción

La función scatterplot del paquete car nos permite crear diagramas de dispersión mejorados. Incluye diagramas de caja en los márgenes, rectas de regresión suavizadas, identificación de valores atípicos etc. Veremos varios ejemplos usando datos (data frames) incluidos en el propio paquete car: Prestige y UN.

library(car)
head(Prestige)

                     education income women prestige census type
gov.administrators      13.11  12351 11.16     68.8   1113 prof
general.managers        12.26  25879  4.02     69.1   1130 prof
accountants             12.77   9271 15.70     63.4   1171 prof
purchasing.officers     11.42   8865  9.11     56.8   1175 prof
chemists                14.62   8403 11.68     73.5   2111 prof
physicists              15.64  11030  5.13     77.6   2113 prof

Ejemplos

1. Gráfico relacionando la renta y la educación, incluyendo diagramas de caja en los ejes y regresión local (línea LOESS). Empleamos el formato fórmula y ~ x para indicar las variables a representar.

scatterplot(education ~ income,
            data = Prestige,
            pch = 16,
            col = "darkblue",
            main = "Educación y Renta\n del conjunto de datos \"Prestige\"",
            xlab = "Renta (dollares)",
            ylab = "Educación (años)")


2. Añadiendo elipses donde se concentran los datos

scatterplot(prestige ~ income, data = Prestige, ellipse = TRUE)


3. Representamos grupos siguiendo el formato y ~ x | z, en el que z evalúa como un factor otra variable pare dividir los datos en grupos.

scatterplot(prestige ~ income|type, data = Prestige, legend.coords = "topleft")


4. Etiquetar un número determinado de puntos.

Por defecto id.n = 0, lo especificamos para indicar el número de puntos que deseamos etiquetar.

scatterplot(infant.mortality ~ gdp, log = "xy", data = UN, id.n = 5)


5. Etiquetar los datos interactivamente.

Utilizamos el argumento id.method para etiquetar de la misma manera que hicimos con la función identify.

scatterplot(infant.mortality ~ gdp, id.method = "identify", data = UN)

Entradas relacionadas

• Regresión lineal simple para principiantes en R
• Diagramas de dispersión

Referencias

• Manual del paquete car

Ajustar el espacio entre el gráfico y la leyenda en ggplot2

Problema

Queremos aumentar la distancia de separación entre el gráfico y la legenda usando el paquete ggplot2.

library(ggplot2)
xy <- data.frame(x = 1:10, y = 10:1, type = rep(LETTERS[1:2], each=5))
plot <- ggplot(data = xy) +
        geom_point(aes(x = x, y = y, colo r= type))

Solución

plot <- plot + theme(legend.box.margin = margin(0, 0, 0, 20))

Referencias

• stackoverflow

Contar el número de ocurrencias basado en dos columnas en R

Problema

Deseamos añadir un contador con el número de ocurrencias basado en las columnas timey site. Con la particularidad de que la agrupación de time sea por día, si no, cada fila sería una ocurrencia distinta. Es decir, cada vez que la combinación de día y site diferente, reiniciará el contador.

                  time site val
1  2014-09-01 00:00:00 2001   1
2  2014-09-01 00:15:00 2001   0
3  2014-09-01 00:30:00 2001   2
4  2014-09-01 00:45:00 2001   0
5  2014-09-01 00:00:00 2002   1
6  2014-09-01 00:15:00 2002   0
7  2014-09-01 00:30:00 2002   2
8  2014-09-02 00:45:00 2001   0
9  2014-09-02 00:00:00 2001   1
10 2014-09-02 00:15:00 2001   0
11 2014-09-02 00:30:00 2001   2
12 2014-09-02 00:45:00 2001   0
13 2014-09-02 00:00:00 2002   1
14 2014-09-02 00:15:00 2002   0
15 2014-09-02 00:30:00 2002   2
16 2014-09-02 00:45:00 2001   0

df <- structure(list(time = structure(c(1409522400, 1409523300, 1409524200, 1409525100, 1409522400, 1409523300, 1409524200, 1409611500, 1409608800, 1409609700, 1409610600, 1409611500, 1409608800, 1409609700, 1409610600, 1409611500), class = c("POSIXct", "POSIXt"), tzone = ""), site = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L), .Label = c("2001", "2002"), class = "factor"), val = c(1L, 0L, 2L, 0L, 1L, 0L, 2L, 0L, 1L, 0L, 2L, 0L, 1L, 0L, 2L, 0L)), .Names = c("time", "site", "val"), row.names = c(NA, -16L), class = "data.frame")

Solución

Utilizamos dplyr. Primero creamos una columna id, extrayendo el día de la fecha (columna time). Después agrupamos por site e id, y añadimos una nueva variable countercon el número de ocurrencias para esos dos grupos.

library(dplyr)
df$id <- as.factor(format(df$time,'%d'))
library(dplyr)
df %>% group_by(site, id) %>% mutate(counter = row_number()) 

Ahora tenemos un contador que empieza desde 1 cada vez que el día y el site cambian. Si se vuelven a repetir, como en la última fila, continua con la numeración.

Source: local data frame [16 x 5]
Groups: site, id

                  time site val id counter
1  2014-09-01 00:00:00 2001   1 01       1
2  2014-09-01 00:15:00 2001   0 01       2
3  2014-09-01 00:30:00 2001   2 01       3
4  2014-09-01 00:45:00 2001   0 01       4
5  2014-09-01 00:00:00 2002   1 01       1
6  2014-09-01 00:15:00 2002   0 01       2
7  2014-09-01 00:30:00 2002   2 01       3
8  2014-09-02 00:45:00 2001   0 02       1
9  2014-09-02 00:00:00 2001   1 02       2
10 2014-09-02 00:15:00 2001   0 02       3
11 2014-09-02 00:30:00 2001   2 02       4
12 2014-09-02 00:45:00 2001   0 02       5
13 2014-09-02 00:00:00 2002   1 02       1
14 2014-09-02 00:15:00 2002   0 02       2
15 2014-09-02 00:30:00 2002   2 02       3
16 2014-09-02 00:45:00 2001   0 02       6

Referencias

• stackoverflow

Convertir de formato largo a ancho usando una función de agregación en R

Problema

Deseamos transformar un data frame de formato largo a ancho agrupando por la suma de los valores de una de las variables. Es decir pasar de esta tabla,

 V1 V3     V2
1  5 10    low
2  5  3    low
3  5  6 medium
4 45 10    low
5 45  3    low
6 77  1   high

A esta otra, donde mantenemos la columna V1, la columna V2 contiene los nombres de las nuevas columnas, y la columna V3 será la variable que almacenará los valores y que agregaremos sumando.

  V1 low medium high
1  5  13      6    0
2 45  13      0    0
3 77   0      0    1

Datos originales

V1 <- c(5,5,5,45,45,77)  
V2 <- c("low", "low", "medium", "low", "low", "high")  
V3 <- c(10,3,6,10,3,1)  
df <- as.data.frame(cbind(V1,V3,V2))

Soluciones

Primero examinamos el data frame que nos dieron.

str(df)

'data.frame': 6 obs. of  3 variables:
 $ V1: Factor w/ 3 levels "45","5","77": 2 2 2 1 1 3
 $ V3: Factor w/ 4 levels "1","10","3","6": 2 3 4 2 3 1
 $ V2: Factor w/ 3 levels "high","low","medium": 2 2 3 2 2 1

Las columnas V1 y V3 se han convertido en factores, cuando deberían ser númericas. Si utilizáramos cualquier código obtendríamos el error: Error: ‘sum’ not meaningful for factors. Para solucionarlo, empleamos el siguiente código:

df <- cbind.data.frame(V1,V3,V2)
str(df)

'data.frame': 6 obs. of  3 variables:
 $ V1: num  5 5 5 45 45 77
 $ V3: num  10 3 6 10 3 1
 $ V2: Factor w/ 3 levels "high","low","medium": 2 2 3 2 2 1

reshape2

Empleamos dcast para crear un data frame, usando V3 como la variable valores, e indicamos como función de agregación la suma (sum).

library(reshape2)
dcast(df, V1  ~ V2, value.var="V3", fun = sum)

  V1 high low medium
1  5    0  13      6
2 45    0  13      0
3 77    1   0      0

Como se puede observar, ordenado las columnas alfabéticamente (high, low y medium) en lugar de respetar el orden de ocurrencia original: low, medium, high.

dcast(df, V1  ~ factor(V2, levels = unique(V2)), value.var = "V3", sum)

'data.frame': 6 obs. of  3 variables:
 $ V1: num  5 5 5 45 45 77
 $ V3: num  10 3 6 10 3 1
 $ V2: Factor w/ 3 levels "high","low","medium": 2 2 3 2 2 1

paquete base

Empleamos la función xtabs para crear una tabla de contingencia mediante el formato fórmula.

xtabs(V3 ~ V1 + V2, df)

    V2
V1   high low medium
  5     0  13      6
  45    0  13      0
  77    1   0      0

xtabs(V3 ~ V1 + factor(V2, c("low", "medium", "high")), df)

   factor(V2, c("low", "medium", "high"))
V1   low medium high
  5   13      6    0
  45  13      0    0
  77   0      0    1

dplyr y tidyr

Con dplyr con la función spread, el equivalente a dcast en reshape2.

df$V2 <- factor(df$V2, levels = c("low", "medium", "high"))

library(tidyr)# Función spread
library(dplyr)

df %>%
  group_by(V1, V2) %>%
  summarise(sum = sum(V3)) %>%
  spread(V2, sum, fill = 0)

Source: local data frame [3 x 4]

     V1   low medium  high
  (dbl) (dbl)  (dbl) (dbl)
1     5    13      6     0
2    45    13      0     0
3    77     0      0     1

Entradas relacionadas

• Tablas de contingencia en R
• Transformar tabla de contingencia en tabla de de frecuencias

Referencias

• stackoveflow