2015-11-29

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

Title 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

Referencias

2015-11-25

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

Title

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

    2015-11-23

    Diagramas de dispersión con el paquete car en R

    Title

    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.
    2. 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)")
      
    3. Añadiendo elipses donde se concentran los datos
    4. scatterplot(prestige ~ income, data = Prestige, ellipse = TRUE)
      

    5. 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.
    6. scatterplot(prestige ~ income|type, data = Prestige, legend.coords = "topleft")

    7. Etiquetar un número determinado de puntos.

    8. 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)

    9. Etiquetar los datos interactivamente.

    10. 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

    Referencias

    2015-11-18

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

    Title

    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

    2015-11-14

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

    Title

    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

    2015-11-11

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

    Title

    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

    Referencias

    Nube de datos