R para profesionales de la salud
Juan José Lasarte Velillas. Pediatra
06 April, 2023
Introducción
¿Por qué R?
Es un programa gratuito y muy potente para análisis de datos que está siendo utilizado por docentes en la mayoría de Universidades y centros de investigación. También la utiliza Google, Facebook, Microsoft, IBM…
Existen librarías para hacer gran variedad de técnicas estadísticas y gráficas que no pueden hacerse con otros programas como SPSS y que facilitan mucho el trabajo.
El tutorial
Este tutorial es un brevísimo resumen del manejo de R que puede ser de utilidad para profesionales sanitarios que no tenemos grandes conocimientos matemáticos. Pretende ser un minimanual para poder recordar la sintaxis de los cálculos más frecuentes en ciencias de la salud. La mayoría de lo que se recoge aquí está más explicado en este curso: https://isglobal-brge.github.io/curso_R/
En este documento, lo que hay en recuadros grises es el código que hay que escribir. En los recuadros blancos se ven los resultados tras ejecutar el código.
Para seguir este tutorial te recomiendo que vayas copiando, pegando y modificando el código en tu R-Studio para ver los resultados.
Instalar R en Windows y Mac
Primero se descarga e instala R desde estos enlaces:
Windows
https://cran.r-project.org/bin/windows/base/
IOS
https://cran.r-project.org/bin/macosx/
Elegir el archivo pkg según el sistema operativo (Sierra, Big Sur…)
Después se descarga e instala RStudio (¡versión
free!), que es el programa que utilizamos para manejar R
https://www.rstudio.com/products/rstudio/download/
Esta instalación está explicada en los primeros 4 minutos de este vídeo:
Manejo básico de RStudio
En R studio trabajamos escribiendo el código en un fichero script, que tiene una extendión .R (se abre en la ventana superior izquierda). Aquí podemos escribir y guardar el código en el fichero lo que facilita mucho el trabajo.
Luego el código se ejecuta en la consola (inferior izda) tras el símbolo >. Para ejecutar el código puede hacerse línea a línea con el botón Run o con Ctrl+Enter (Windows) o Cmd+Enter (Mac). En ambos casos se ejecuta la línea en la que está el cursor. También puede ejecutarse un bloque de código seleccionándolo.
Cuando realizamos una acción mediante un menú, en esta consola se ve el código que se ha ejecutado y podemos copiarlo para utilizarlo en nuestro script. Siempre es más rápido ejecutar tareas con código que con un menú.
Para abrir un fichero script: File > New file > R Script
R es un lenguaje de programación dirigido a objetos (conforme se crean, pueden verse en la ventana Environment). En la pestaña History se guardan las líneas de código que se han ejecutado.
En la ventana inferior derecha pueden verse los Plots (gráficos que vamos creando), los Packages (paquetes instalados), la ayuda Help y un Viewer en donde se ve lo que se escribe en ficheros RMarkdown. También hay una pestaña Files que es un explorador de archivos.
Un objeto puede contener cualquier cosa, por ejemplo, creamos un objeto con este código:
objeto1 <- 4+5el símbolo <- equivale a = pero en R es mejor escribirlo de esa manera para no confundirnos en otros lugares en donde se emplea el signo = para otras funciones
Para ver el objeto:
objeto1[1] 9Otro ejemplo que crea un vector aleatorio con la función normal
objeto2<- rnorm(10000, mean = 0, sd = 1)Y que podemos representarlo en una gráfica:
hist(objeto2)
Órdenes básicas
Establecer el directorio de trabajo
Puede hacerse con este código o desde el menú Session –> Set Working Directory
setwd("ruta del directorio de trabajo")
Para ver el directorio de trabajo en el que estamos:
getwd()
Instalar y cargar librerías
Hay funciones ya precargadas en R pero otras funciones es necesario instalarlas (esto solo hay que hacerlo una vez) y cargarlas cada vez que las utilizamos.
Pueden instalarse librerías desde el menú Tools o bien con este código:
install.packages()
Por ejemplo:
install.packages("summarytools")
Luego, para cargar la librería cuando tengamos que utilizarla, sin comillas:
library(summarytools)Ayudas en cualquier momento
Para obtener ayuda sobre una función basta con poner un signo ? delante del nombre de la función y luego consultar la información en los links que ofrece.
?meanHelp on topic 'mean' was found in the following packages:
Package Library
base /Library/Frameworks/R.framework/Resources/library
rapportools /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library
Using the first match ...Para obtener ayuda de una librería
library(help="summarytools")También hay webs con mucha información y ayuda:
R-Bloggers
Stack Overflow…
Más adelante, cuando se conocen las librerías, hay hojas resumen
(cheetsheets) de funciones en esta página, las básicas con versión en
español:
https://www.rstudio.com/resources/cheatsheets/
Importar datos
En ciencias de la salud normalmente trabajamos con dataframes que son objetos en los que las filas son las observaciones y las columnas las variables.
La mejor forma de crear un datagrame es leyendo (importando) los datos de un archivo.
Habitualmente tenemos que importar estos datos de una hoja excel o de otro tipo de archivo. Lo que mejor funciona para importar son los archivos .csv o .txt separados por tabulaciones.
Puede importarse un archivo desde la pestaña Environment,
que nos facilita la tarea y podemos ver como quedará el dataframe, o
bien con código:
Recuerda que si la primera vez lo haces con el menú, luego puedes copiar
el código de la cónsola para que la próxima vez sea más rápido.
data <-read.delim(file.choose(),header=T, sep=",")
Para un CSV, con el argumento file.choose() elegimos el archivo. Si sabemos la ruta, podemos ponerla directamente entre comillas (“ruta”), el argumento header especifica que la primera fila son los nombres de las variables y sep la separación entre los valores de cada fila, habitualmente comas en un CSV.
En este documento trabajaremos con el archivo
antro.csv, así que primero hay que importarlo:
Puedes descargarlo en: https://anonfiles.com/Fc51v0Wbxc/antro_csv
antro <- read.csv2("~/Documents/TRABAJO/Curso R/antro.csv")También puede exportarse un dataframe a CSV separado por comas y con los nombres de las variables
write.table(dataframe, file = "fichero.CSV", sep = ",", col.names = TRUE, row.names = FALSE)
Funciones básicas
Ver los nombres de las variables y las primeras observaciones:
head(antro) zona fechanac sexo peso talla fechavisit ip rango ciap_ap ciap_imc
1 Zona4 1/2/04 mujer 50.5 155.5 16/2/17 -1.67 Normopeso SI NO
2 Zona4 12/2/04 mujer 33.9 149.5 16/2/16 -1.67 Normopeso NO NO
3 Zona4 28/2/04 mujer 27.5 132.0 7/2/14 -1.67 Normopeso NO NO
4 Zona4 11/4/04 mujer 33.4 148.0 1/7/16 -1.67 Normopeso NO NO
5 Zona4 11/4/04 mujer 33.9 149.0 1/7/16 -1.67 Normopeso NO NO
6 Zona4 30/4/04 mujer 19.5 110.0 7/5/10 -1.67 Normopeso NO NO
Hb creat edad1sex edad1emb edaddiag edadmadre niveledu tabaco ets acos
1 15.00 1 16 16 40 64 primaria ex-fumador si no
2 15.00 1 14 14 NA 48 ninguno fumador no no
3 14.99 1 23 24 30 41 primaria no fumador no no
4 14.99 1 29 29 45 51 primaria no fumador no si
5 14.99 1 22 22 42 58 primaria no fumador no no
6 14.98 1 25 28 NA 57 primaria ex-fumador no si
puntostest familia
1 9 negativo
2 11 positivo
3 3 negativo
4 1 negativo
5 8 negativo
6 5 negativoVer los nombres de las columnas (variables):
colnames(antro) [1] "zona" "fechanac" "sexo" "peso" "talla"
[6] "fechavisit" "ip" "rango" "ciap_ap" "ciap_imc"
[11] "Hb" "creat" "edad1sex" "edad1emb" "edaddiag"
[16] "edadmadre" "niveledu" "tabaco" "ets" "acos"
[21] "puntostest" "familia" Ver las variables y su tipo (numéricas, carácter…):
str(antro)'data.frame': 10551 obs. of 22 variables:
$ zona : chr "Zona4" "Zona4" "Zona4" "Zona4" ...
$ fechanac : chr "1/2/04" "12/2/04" "28/2/04" "11/4/04" ...
$ sexo : chr "mujer" "mujer" "mujer" "mujer" ...
$ peso : num 50.5 33.9 27.5 33.4 33.9 19.5 45.4 17.7 51.9 50.1 ...
$ talla : num 156 150 132 148 149 ...
$ fechavisit: chr "16/2/17" "16/2/16" "7/2/14" "1/7/16" ...
$ ip : num -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 ...
$ rango : chr "Normopeso" "Normopeso" "Normopeso" "Normopeso" ...
$ ciap_ap : chr "SI" "NO" "NO" "NO" ...
$ ciap_imc : chr "NO" "NO" "NO" "NO" ...
$ Hb : num 15 15 15 15 15 ...
$ creat : num 1 1 1 1 1 1 1 1 1 1 ...
$ edad1sex : int 16 14 23 29 22 25 33 36 18 17 ...
$ edad1emb : int 16 14 24 29 22 28 33 36 18 17 ...
$ edaddiag : int 40 NA 30 45 42 NA NA NA 23 45 ...
$ edadmadre : int 64 48 41 51 58 57 64 56 41 52 ...
$ niveledu : chr "primaria" "ninguno" "primaria" "primaria" ...
$ tabaco : chr "ex-fumador" "fumador" "no fumador" "no fumador" ...
$ ets : chr "si" "no" "no" "no" ...
$ acos : chr "no" "no" "no" "si" ...
$ puntostest: int 9 11 3 1 8 5 2 6 3 13 ...
$ familia : chr "negativo" "positivo" "negativo" "negativo" ...Número de columnas (variables)
ncol(antro)[1] 22Número de filas (observaciones)
nrow(antro)[1] 10551Número de filas y columnas
dim(antro)[1] 10551 22Las variables se escriben con el nombre del data frame, el símbolo $ y el nombre de la variable, por ejemplo: antro$talla
Contar el número de registros (filas) que cumplen una condición (ver apartado operadores lógicos)
sum(antro$sexo=="mujer")[1] 5081Descritivo para una variable categórica table
table(antro$sexo)
mujer varon
5081 5470 Descriptivo para una variable continua summary
summary(antro$creat) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.060 0.550 1.500 1.407 2.240 3.000 8 Media de una variable continua
mean(antro$creat,na.rm = T)[1] 1.407486(na.rm=T se pone para evitar errores por ausencia de datos en alguna casilla)
Redondear a un número de decimales, por ejemplo a 1
mean(antro$creat, na.rm = T)[1] 1.407486round(mean(antro$creat, na.rm = T),1)[1] 1.4Truncar (muestra solo el valor entero)
mean(antro$talla, na.rm = T)[1] 123.0786trunc(mean(antro$talla, na.rm = T))[1] 123Listar un número de filas, las 2 primeras, por ejemplo
antro [1:2, ] zona fechanac sexo peso talla fechavisit ip rango ciap_ap ciap_imc
1 Zona4 1/2/04 mujer 50.5 155.5 16/2/17 -1.67 Normopeso SI NO
2 Zona4 12/2/04 mujer 33.9 149.5 16/2/16 -1.67 Normopeso NO NO
Hb creat edad1sex edad1emb edaddiag edadmadre niveledu tabaco ets acos
1 15 1 16 16 40 64 primaria ex-fumador si no
2 15 1 14 14 NA 48 ninguno fumador no no
puntostest familia
1 9 negativo
2 11 positivoEsta nomenclatura con corchetes representa [filas, columnas]
Símbolos y atajos del teclado
Consultar más en el menú Tools ->Keyboard Shortcuts help
El símbolo # sirve para añadir texto que no se ejecutará como código y es útil para añadir explicaciones al código. En un script, cuando detrás de ese texto se añade #### se crean títulos que pueden contraerse y expandirse.
| Teclas Windows | Teclas MAC | Símbolo |
|---|---|---|
| Alt-Ñ | Option-Ñ | ~ |
| Ctr+May+M | Cmd+May+M | %>% |
| Ctrl-L | Ctrl-L | Borra la consola |
| Ctrl+Enter | Cmd+Return | Run línea o selección |
| Ctrl+Alt+R | Cmd+Option+R | Run todo el script actual |
| Ctrl+Alt+B | Cmd+Option+B | Run desde el inicio hasta la línea actual |
| Ctrl+Alt+E | Cmd+Option+E | Run desde la línea actual hasta el final del documento |
| Ctrl+Alt+i | Cmd+Option+i | Insertar chunk |
Operadores lógicos: Y = &, O = |, NO = !
| Operadores relacionales | Mayor | Menor | mayor o igual | menor o igual | igual | distinto |
|---|---|---|---|---|---|---|
| . | > | < | >= | <= | == | != |
Creación y recodificación de variables
Crear variables
Crear variable “NA”
antro$variable1 <- NA
Crear variable con casillas con un “1”
antro$variable2 <- 1
Crear variable con numeración ascendente desde 1
antro$variable3 <- 1:nrow(antro)
Crear variable con condiciones
Si la edad es = 10, pondrá 1 y si no 0
antro$variable4 <- ifelse(antro$edad==10,1,0)
Si el sexo es M pondrá 1 y si no 0
antro$variable5 <- ifelse(antro$sexo=="M",1,0)
Crear variable a partir de otras variables
Por ejemplo, para crear la varible imc
antro$imc <- antro$peso/((antro$talla/100)^2)Esto también puede hacerse con la librería tidyverse (instalarla y cargarla antes).
library(tidyverse)Por ejemplo,
mutate(antro, imc = peso/((talla/100)^2))
Recodificar variable en cuartiles
Por ejemplo, la variable talla,
antro <- mutate(antro, talla_q = cut(talla, quantile(talla, na.rm=TRUE)))
table(antro$talla_q)
(80,108] (108,119] (119,137] (137,189]
2660 2732 2542 2605 NOTA: podemos averiguar los percentiles X de una variable:
quantile(antro$peso, probs=c(0.5, 0.97, 0.98, 0.99), na.rm = T) 50% 97% 98% 99%
23.100 58.000 61.548 68.000 Recodificar variable según puntos de corte y ponerles etiqueta.
Ejemplo: talla menor de 100, entre los 100 y 130 (ambos incluidos) y más de 130,
antro <- mutate(antro, talla_groups = cut(talla, c(-Inf, 100, 130, Inf),))
table(antro$talla_groups)
(-Inf,100] (100,130] (130, Inf]
1044 5970 3526 Lo mismo, poniendo etiquetas,
antro <- mutate(antro, talla_etiq = cut(talla, c(-Inf, 100, 130, Inf),
labels=c("<100", "100-130", ">130")))
table(antro$talla_etiq)
<100 100-130 >130
1044 5970 3526 Recodificar variable con categorías en otras categorías diferentes.
Ejemplo: sexo hombre o mujer por M / F
antro <- mutate(antro,
sex = recode(sexo,
"varon" = "M",
"mujer" = "F"))
table(antro$sex)
F M
5081 5470 Transformar el tipo de variable
Cambiar a tipo fecha (fijarse que en as.Date es D mayúscula) (“Y” mayúscula para años en 4 cifras)
antro$fechanac <- as.Date(antro$fechanac, format="%d/%m/%y")| Símbolo | Significado |
|---|---|
| %d | día (numérico, de 0 a 31) |
| %a | día de la semana abreviado a tres letras |
| %A | día de la semana (nombre completo) |
| %m | mes (numérico de 0 a 12) |
| %b | mes (nombre abreviado a tres letras) |
| %B | mes (nombre completo) |
| %y | año (con dos dígitos) |
| %Y | año (con cuatro dígitos) |
Cambiar a tipo factor
antro$sexo <- as.factor(antro$sexo)Cambiar a tipo numérica
antro$edaddiag <- as.numeric(antro$edaddiag)Cambiar a tipo carácter
dataframe$variable <- as.character(dataframe$variable)
Cambiar a tipo número entero
dataframe$variable <- as.integer(dataframe$variable)
Añadir etiquetas a las variables
Puede hacerse con la librería labelled
library(labelled)var_label(antro$fechanac) <- "Fecha de nacimiento"
View(antro)También pueden añadirse varias etiquetas al mismo tiempo
var_label(antro) <- list(zona="zona básica de salud", Hb = "Hemoglobina", creat = "Creatinina")
View(antro)Para ver las etiquetas (las variables sin etiquetas se muestran como “NULL”)
var_label(antro)$zona
[1] "zona básica de salud"
$fechanac
[1] "Fecha de nacimiento"
$sexo
NULL
$peso
NULL
$talla
NULL
$fechavisit
NULL
$ip
NULL
$rango
NULL
$ciap_ap
NULL
$ciap_imc
NULL
$Hb
[1] "Hemoglobina"
$creat
[1] "Creatinina"
$edad1sex
NULL
$edad1emb
NULL
$edaddiag
NULL
$edadmadre
NULL
$niveledu
NULL
$tabaco
NULL
$ets
NULL
$acos
NULL
$puntostest
NULL
$familia
NULL
$imc
NULL
$talla_q
NULL
$talla_groups
NULL
$talla_etiq
NULL
$sex
NULLSeleccionar registros para hacer subgrupos
Pueden utilizarse la función filter o subset. La función filter (de la librería tidyverse) no funciona si hay NAs, mejor utilizar subset()
varones <- subset(antro, sexo=="varon")
dim(varones)[1] 5470 27varones_menores_100 <- filter(antro, sexo=="varon" & talla < 100)
dim(varones_menores_100)[1] 425 27Cálculos con edades y fechas
Para cálculo de edades utilizar la librería eeptools
library(eeptools)Para el manejo de fechas, utilizar la librería lubridate
library(lubridate)Ejemplo para crear variables de edad en días, meses y años a partir de la fecha de nacimiento y de la fecha de la visita:
Como las variables fechanac y fechavisit son carácter (chr), primero hay que convertirlas en fecha (Date)
antro$fechanac<-as.Date(antro$fechanac, format="%d/%m/%Y")
antro$fechavisit<-as.Date(antro$fechavisit, format="%d/%m/%y")y luego utilizar la función interval de lubridate
antro$dias<-interval(antro$fechanac,antro$fechavisit)%/% days(1)
antro$meses<- interval(antro$fechanac,antro$fechavisit)%/% months(1)
antro$anos<-interval(antro$fechanac,antro$fechavisit)%/% years(1)
str(antro)'data.frame': 10551 obs. of 30 variables:
$ zona : chr "Zona4" "Zona4" "Zona4" "Zona4" ...
..- attr(*, "label")= chr "zona básica de salud"
$ fechanac : Date, format: "2004-02-01" "2004-02-12" ...
$ sexo : Factor w/ 2 levels "mujer","varon": 1 1 1 1 1 1 1 1 1 1 ...
$ peso : num 50.5 33.9 27.5 33.4 33.9 19.5 45.4 17.7 51.9 50.1 ...
$ talla : num 156 150 132 148 149 ...
$ fechavisit : Date, format: "2017-02-16" "2016-02-16" ...
$ ip : num -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 ...
$ rango : chr "Normopeso" "Normopeso" "Normopeso" "Normopeso" ...
$ ciap_ap : chr "SI" "NO" "NO" "NO" ...
$ ciap_imc : chr "NO" "NO" "NO" "NO" ...
$ Hb : num 15 15 15 15 15 ...
..- attr(*, "label")= chr "Hemoglobina"
$ creat : num 1 1 1 1 1 1 1 1 1 1 ...
..- attr(*, "label")= chr "Creatinina"
$ edad1sex : int 16 14 23 29 22 25 33 36 18 17 ...
$ edad1emb : int 16 14 24 29 22 28 33 36 18 17 ...
$ edaddiag : num 40 NA 30 45 42 NA NA NA 23 45 ...
$ edadmadre : int 64 48 41 51 58 57 64 56 41 52 ...
$ niveledu : chr "primaria" "ninguno" "primaria" "primaria" ...
$ tabaco : chr "ex-fumador" "fumador" "no fumador" "no fumador" ...
$ ets : chr "si" "no" "no" "no" ...
$ acos : chr "no" "no" "no" "si" ...
$ puntostest : int 9 11 3 1 8 5 2 6 3 13 ...
$ familia : chr "negativo" "positivo" "negativo" "negativo" ...
$ imc : num 20.9 15.2 15.8 15.2 15.3 ...
$ talla_q : Factor w/ 4 levels "(80,108]","(108,119]",..: 4 4 3 4 4 2 4 2 4 4 ...
$ talla_groups: Factor w/ 3 levels "(-Inf,100]","(100,130]",..: 3 3 3 3 3 2 3 2 3 3 ...
$ talla_etiq : Factor w/ 3 levels "<100","100-130",..: 3 3 3 3 3 2 3 2 3 3 ...
$ sex : chr "F" "F" "F" "F" ...
$ dias : num 4764 4387 3632 4464 4464 ...
$ meses : num 156 144 119 146 146 72 154 72 145 154 ...
$ anos : num 13 12 9 12 12 6 12 6 12 12 ...Función “age_calc” calcula la edad entre 2 fechas, pero requiere que ninguna variable tenga datos “NA”
antro <- na.omit(antro)Para calcular la edad decimal en años, a día de hoy:
antro$anoshoy <- age_calc(antro$fechanac, today(), units="years")Para calcular la edad decimal en años, a una fecha concreta:
antro$anosdec <- age_calc(antro$fechanac, as.Date("2022-05-04"), units="years")Algunos trucos NO imprescindibles
Evitar escribir el nombre del dataframe antes del nombre de la variable
attach(dataframe)
Para volver a dirigir la salida a la pantalla
detach()
Guardar información NO grafica en un fichero de texto
sink("fichero.txt")
Para redirigir otra vez la salida a la pantalla
sink()
Ejecutar un script dentro de otro script
El scrpt tiene que estar en el directorio de trabajo o poner toda la ruta completa
source("script.R", echo = T)
Borrar variables (columnas)
dataframe$variable = NULL
Cambiar nombres de las columnas (variables)
dataframe <- rename(dataframe, variablenew= variableold)
Eliminar filas con “NA”
A veces al importar aparecen filas sin contenido y con “NA”
Para eliminar estas filas:
dataframe <- na.omit(dataframe)
En otras ocasiones nos interesará dejar en blanco las casillas “NA”
dataframe[is.na(dataframe)] <- ""
Añadir variables (columnas) mediante una variable común
La variable común en ambos archivos tiene que ser la misma con el mismo nombre y tipo. Por ejemplo, es útil para añadir el valor de los percentiles por edad (en días o meses, que será la variable de unión) a otro archivo que contiene los registros con la edad.
archivofinal <- left_join(archivo1, archivo2, by = "nombre variable de unión")
Unir registros (filas) de dos bases de datos
Ambas bases de datos tienen que contener las mismas variables (columnas)
archivofinal <- rbind(archivo1,archivo2)
Estadística descriptiva
Necesitaremos estas 2 librerías:
library(summarytools)
library(compareGroups)Y cargar el archivo con el que se va a trabajar, “antro.csv”
antro <- read.csv2("~/Documents/TRABAJO/Curso R/antro.csv")Variables categóricas
freq(antro$sexo)Frequencies
antro$sexo
Type: Character
Freq % Valid % Valid Cum. % Total % Total Cum.
----------- ------- --------- -------------- --------- --------------
mujer 5081 48.16 48.16 48.16 48.16
varon 5470 51.84 100.00 51.84 100.00
<NA> 0 0.00 100.00
Total 10551 100.00 100.00 100.00 100.00Si queremos evitar los NA en caso de que existan
freq(antro$zona, report.nas = FALSE, headings = FALSE)
Freq % % Cum.
----------- ------- -------- --------
Zona1 1426 13.52 13.52
Zona2 1860 17.63 31.14
Zona3 1694 16.06 47.20
Zona4 5571 52.80 100.00
Total 10551 100.00 100.00Para todas las variables categóricas (para evitar que tabule las fechas, antes hay que convertirlas a formato fecha)
antro$fechanac <- as.Date(antro$fechanac, format="%d/%m/%y")
antro$fechavisit <- as.Date(antro$fechavisit, format="%d/%m/%y")freq(antro, report.nas = FALSE, headings = FALSE)antro$zona
Freq % % Cum.
----------- ------- -------- --------
Zona1 1426 13.52 13.52
Zona2 1860 17.63 31.14
Zona3 1694 16.06 47.20
Zona4 5571 52.80 100.00
Total 10551 100.00 100.00
antro$sexo
Freq % % Cum.
----------- ------- -------- --------
mujer 5081 48.16 48.16
varon 5470 51.84 100.00
Total 10551 100.00 100.00
antro$ip
Freq % % Cum.
----------- ------- -------- --------
-2.33 1694 16.06 16.06
-1.67 5571 52.80 68.86
-0.2 1426 13.52 82.37
1.13 1860 17.63 100.00
Total 10551 100.00 100.00
antro$rango
Freq % % Cum.
--------------- ------- -------- --------
Normopeso 7178 68.03 68.03
Obesidad 1368 12.97 81.00
Sobrepeso 2005 19.00 100.00
Total 10551 100.00 100.00
antro$ciap_ap
Freq % % Cum.
-------------------- ------- -------- --------
(Empty string) 11 0.10 0.10
NO 10214 96.81 96.91
SI 326 3.09 100.00
Total 10551 100.00 100.00
antro$ciap_imc
Freq % % Cum.
-------------------- ------- --------- ---------
(Empty string) 4 0.038 0.038
NO 7176 68.013 68.050
SI 3371 31.950 100.000
Total 10551 100.000 100.000
antro$niveledu
Freq % % Cum.
------------------- ------- -------- --------
ninguno 2434 23.07 23.07
primaria 5067 48.02 71.09
secundaria 2593 24.58 95.67
ticnico 230 2.18 97.85
universitario 227 2.15 100.00
Total 10551 100.00 100.00
antro$tabaco
Freq % % Cum.
---------------- ------- -------- --------
ex-fumador 779 7.38 7.38
fumador 980 9.29 16.67
no fumador 8792 83.33 100.00
Total 10551 100.00 100.00
antro$ets
Freq % % Cum.
----------- ------ -------- --------
no 6677 67.94 67.94
si 3151 32.06 100.00
Total 9828 100.00 100.00
antro$acos
Freq % % Cum.
----------- ------- -------- --------
no 7175 68.00 68.00
si 3376 32.00 100.00
Total 10551 100.00 100.00
antro$puntostest
Freq % % Cum.
----------- ------- --------- ---------
1 575 5.710 5.710
2 1254 12.453 18.163
3 1422 14.121 32.284
4 1393 13.833 46.117
5 1183 11.748 57.865
6 1047 10.397 68.262
7 870 8.640 76.902
8 662 6.574 83.476
9 466 4.628 88.103
10 373 3.704 91.807
11 266 2.642 94.449
12 223 2.214 96.663
13 110 1.092 97.756
14 79 0.785 98.540
15 77 0.765 99.305
16 26 0.258 99.563
17 22 0.218 99.782
18 10 0.099 99.881
20 6 0.060 99.940
22 3 0.030 99.970
25 3 0.030 100.000
Total 10070 100.000 100.000
antro$familia
Freq % % Cum.
-------------- ------ -------- --------
negativo 5222 59.23 59.23
positivo 3595 40.77 100.00
Total 8817 100.00 100.00Tabla contingencia
ctable(antro$zona, antro$sexo, prop="r")Cross-Tabulation, Row Proportions
zona * sexo
Data Frame: antro
------- ------ -------------- -------------- ----------------
sexo mujer varon Total
zona
Zona1 712 (49.9%) 714 (50.1%) 1426 (100.0%)
Zona2 878 (47.2%) 982 (52.8%) 1860 (100.0%)
Zona3 814 (48.1%) 880 (51.9%) 1694 (100.0%)
Zona4 2677 (48.1%) 2894 (51.9%) 5571 (100.0%)
Total 5081 (48.2%) 5470 (51.8%) 10551 (100.0%)
------- ------ -------------- -------------- ----------------El argumento prop nos sirve para indicar si queremos las proporciones por fila (‘r’) o columna (‘c’). Podemos eliminar la columna de missings indicando que el argumento useNA sea “no”
ctable(antro$zona, antro$sexo, useNA="no", prop="r")Cross-Tabulation, Row Proportions
zona * sexo
Data Frame: antro
------- ------ -------------- -------------- ----------------
sexo mujer varon Total
zona
Zona1 712 (49.9%) 714 (50.1%) 1426 (100.0%)
Zona2 878 (47.2%) 982 (52.8%) 1860 (100.0%)
Zona3 814 (48.1%) 880 (51.9%) 1694 (100.0%)
Zona4 2677 (48.1%) 2894 (51.9%) 5571 (100.0%)
Total 5081 (48.2%) 5470 (51.8%) 10551 (100.0%)
------- ------ -------------- -------------- ----------------Variables continuas
descr(antro$Hb)Descriptive Statistics
antro$Hb
N: 10551
Hb
----------------- ----------
Mean 12.48
Std.Dev 1.45
Min 10.00
Q1 11.20
Median 12.46
Q3 13.75
Max 15.00
MAD 1.88
IQR 2.55
CV 0.12
Skewness 0.02
SE.Skewness 0.02
Kurtosis -1.22
N.Valid 10543.00
Pct.Valid 99.92Para todas las variables continuas
descr(antro)Descriptive Statistics
antro
N: 10551
creat edad1emb edad1sex edaddiag edadmadre Hb ip
----------------- ---------- ---------- ---------- ---------- ----------- ---------- ----------
Mean 1.41 21.61 20.05 38.72 48.62 12.48 -1.08
Std.Dev 0.92 4.92 4.82 12.99 11.93 1.45 1.18
Min 0.06 10.00 5.00 14.00 20.00 10.00 -2.33
Q1 0.55 18.00 17.00 28.00 40.00 11.20 -1.67
Median 1.50 21.00 19.00 37.00 48.00 12.46 -1.67
Q3 2.24 24.00 23.00 47.00 57.00 13.75 -0.20
Max 3.00 52.00 52.00 80.00 84.00 15.00 1.13
MAD 1.29 4.45 4.45 13.34 13.34 1.88 0.00
IQR 1.69 6.00 6.00 19.00 17.00 2.55 1.47
CV 0.66 0.23 0.24 0.34 0.25 0.12 -1.09
Skewness 0.15 0.90 0.94 0.56 0.13 0.02 0.96
SE.Skewness 0.02 0.02 0.02 0.03 0.02 0.02 0.02
Kurtosis -1.43 0.94 1.63 -0.36 -0.65 -1.22 -0.57
N.Valid 10543.00 10070.00 10452.00 5307.00 10551.00 10543.00 10551.00
Pct.Valid 99.92 95.44 99.06 50.30 100.00 99.92 100.00
Table: Table continues below
peso puntostest talla
----------------- ---------- ------------ ----------
Mean 27.60 5.46 123.08
Std.Dev 12.66 3.24 19.06
Min 8.50 1.00 80.00
Q1 18.60 3.00 108.00
Median 23.10 5.00 119.00
Q3 34.00 7.00 137.00
Max 128.50 25.00 189.00
MAD 8.90 2.97 20.76
IQR 15.40 4.00 29.00
CV 0.46 0.59 0.15
Skewness 1.49 1.02 0.50
SE.Skewness 0.02 0.02 0.02
Kurtosis 2.69 1.22 -0.64
N.Valid 10539.00 10070.00 10540.00
Pct.Valid 99.89 95.44 99.90Descriptiva según una segunda variable
stby(antro, INDICES = antro$sexo,
FUN = descr, stats = "common", transpose = TRUE)Descriptive Statistics
antro
Group: sexo = mujer
N: 5081
Mean Std.Dev Min Median Max N.Valid Pct.Valid
---------------- -------- --------- ------- -------- -------- --------- -----------
creat 0.53 0.28 0.06 0.53 1.01 5076.00 99.90
edad1emb 21.50 4.89 10.00 21.00 43.00 4855.00 95.55
edad1sex 20.02 4.81 5.00 19.00 46.00 5034.00 99.07
edaddiag 38.26 12.56 15.00 36.00 80.00 2306.00 45.38
edadmadre 48.24 11.81 20.00 48.00 82.00 5081.00 100.00
Hb 12.48 1.46 10.00 12.48 15.00 5077.00 99.92
ip -1.09 1.18 -2.33 -1.67 1.13 5081.00 100.00
peso 27.57 12.55 8.50 23.10 94.80 5074.00 99.86
puntostest 5.53 3.26 1.00 5.00 25.00 4855.00 95.55
talla 123.03 19.39 80.00 118.50 183.00 5073.00 99.84
Group: sexo = varon
N: 5470
Mean Std.Dev Min Median Max N.Valid Pct.Valid
---------------- -------- --------- ------- -------- -------- --------- -----------
creat 2.22 0.45 1.44 2.21 3.00 5467.00 99.95
edad1emb 21.71 4.94 10.00 21.00 52.00 5215.00 95.34
edad1sex 20.08 4.83 5.00 19.00 52.00 5418.00 99.05
edaddiag 39.08 13.30 14.00 37.00 80.00 3001.00 54.86
edadmadre 48.97 12.03 20.00 49.00 84.00 5470.00 100.00
Hb 12.47 1.45 10.00 12.43 15.00 5466.00 99.93
ip -1.08 1.19 -2.33 -1.67 1.13 5470.00 100.00
peso 27.63 12.75 10.66 23.20 128.50 5465.00 99.91
puntostest 5.38 3.22 1.00 5.00 25.00 5215.00 95.34
talla 123.12 18.75 88.00 119.00 189.00 5467.00 99.95El argumento stats=“common” lo ponemos para que saque menos estadísticos (sólo los más comunes), pero si no ponemos nada los saca todos y el argumento transpose = TRUE sirve para transponer la tabla de resultados.
Y para todas las variables al mismo tiempo
dfSummary(antro)Data Frame Summary
antro
Dimensions: 10551 x 22
Duplicates: 0
---------------------------------------------------------------------------------------------------------------
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
---- ------------- -------------------------- ---------------------- --------------------- ---------- ---------
1 zona 1. Zona1 1426 (13.5%) II 10551 0
[character] 2. Zona2 1860 (17.6%) III (100.0%) (0.0%)
3. Zona3 1694 (16.1%) III
4. Zona4 5571 (52.8%) IIIIIIIIII
2 fechanac min : 2004-01-01 3959 distinct values . : . 10538 13
[Date] med : 2012-07-16 . . : : : : (99.9%) (0.1%)
max : 2017-12-24 : : : : : : :
range : 13y 11m 23d : : : : : : : :
. : : : : : : : : :
3 sexo 1. mujer 5081 (48.2%) IIIIIIIII 10551 0
[character] 2. varon 5470 (51.8%) IIIIIIIIII (100.0%) (0.0%)
4 peso Mean (sd) : 27.6 (12.7) 705 distinct values : : 10539 12
[numeric] min < med < max: : : (99.9%) (0.1%)
8.5 < 23.1 < 128.5 : :
IQR (CV) : 15.4 (0.5) : : :
: : : : .
5 talla Mean (sd) : 123.1 (19.1) 678 distinct values : 10540 11
[numeric] min < med < max: : : (99.9%) (0.1%)
80 < 119 < 189 . : : .
IQR (CV) : 29 (0.2) : : : : : :
: : : : : : :
6 fechavisit min : 2007-10-31 1390 distinct values : 10539 12
[Date] med : 2019-09-27 : (99.9%) (0.1%)
max : 2020-12-30 : :
range : 13y 1m 30d : :
. : : :
7 ip Mean (sd) : -1.1 (1.2) -2.33 : 1694 (16.1%) III 10551 0
[numeric] min < med < max: -1.67 : 5571 (52.8%) IIIIIIIIII (100.0%) (0.0%)
-2.3 < -1.7 < 1.1 -0.20 : 1426 (13.5%) II
IQR (CV) : 1.5 (-1.1) 1.13 : 1860 (17.6%) III
8 rango 1. Normopeso 7178 (68.0%) IIIIIIIIIIIII 10551 0
[character] 2. Obesidad 1368 (13.0%) II (100.0%) (0.0%)
3. Sobrepeso 2005 (19.0%) III
9 ciap_ap 1. (Empty string) 11 ( 0.1%) 10551 0
[character] 2. NO 10214 (96.8%) IIIIIIIIIIIIIIIIIII (100.0%) (0.0%)
3. SI 326 ( 3.1%)
10 ciap_imc 1. (Empty string) 4 ( 0.0%) 10551 0
[character] 2. NO 7176 (68.0%) IIIIIIIIIIIII (100.0%) (0.0%)
3. SI 3371 (31.9%) IIIIII
11 Hb Mean (sd) : 12.5 (1.5) 501 distinct values : : . . . . . . . . 10543 8
[numeric] min < med < max: : : : : : : : : : : (99.9%) (0.1%)
10 < 12.5 < 15 : : : : : : : : : :
IQR (CV) : 2.6 (0.1) : : : : : : : : : :
: : : : : : : : : :
12 creat Mean (sd) : 1.4 (0.9) 253 distinct values : : : 10543 8
[numeric] min < med < max: : : : . (99.9%) (0.1%)
0.1 < 1.5 < 3 : : : : : : : :
IQR (CV) : 1.7 (0.7) : : : : : : : :
: : : : : : : : : :
13 edad1sex Mean (sd) : 20.1 (4.8) 40 distinct values : 10452 99
[integer] min < med < max: : (99.1%) (0.9%)
5 < 19 < 52 : :
IQR (CV) : 6 (0.2) : : .
: : : : .
14 edad1emb Mean (sd) : 21.6 (4.9) 34 distinct values : 10070 481
[integer] min < med < max: : : (95.4%) (4.6%)
10 < 21 < 52 : : :
IQR (CV) : 6 (0.2) : : : .
. : : : : .
15 edaddiag Mean (sd) : 38.7 (13) 65 distinct values : 5307 5244
[integer] min < med < max: . . : (50.3%) (49.7%)
14 < 37 < 80 : : : :
IQR (CV) : 19 (0.3) : : : : : .
: : : : : : : : .
16 edadmadre Mean (sd) : 48.6 (11.9) 62 distinct values : : 10551 0
[integer] min < med < max: : : : : (100.0%) (0.0%)
20 < 48 < 84 : : : : :
IQR (CV) : 17 (0.2) . : : : : : :
. : : : : : : : .
17 niveledu 1. ninguno 2434 (23.1%) IIII 10551 0
[character] 2. primaria 5067 (48.0%) IIIIIIIII (100.0%) (0.0%)
3. secundaria 2593 (24.6%) IIII
4. ticnico 230 ( 2.2%)
5. universitario 227 ( 2.2%)
18 tabaco 1. ex-fumador 779 ( 7.4%) I 10551 0
[character] 2. fumador 980 ( 9.3%) I (100.0%) (0.0%)
3. no fumador 8792 (83.3%) IIIIIIIIIIIIIIII
19 ets 1. no 6677 (67.9%) IIIIIIIIIIIII 9828 723
[character] 2. si 3151 (32.1%) IIIIII (93.1%) (6.9%)
20 acos 1. no 7175 (68.0%) IIIIIIIIIIIII 10551 0
[character] 2. si 3376 (32.0%) IIIIII (100.0%) (0.0%)
21 puntostest Mean (sd) : 5.5 (3.2) 21 distinct values : 10070 481
[integer] min < med < max: : : : (95.4%) (4.6%)
1 < 5 < 25 : : :
IQR (CV) : 4 (0.6) : : : .
: : : : :
22 familia 1. negativo 5222 (59.2%) IIIIIIIIIII 8817 1734
[character] 2. positivo 3595 (40.8%) IIIIIIII (83.6%) (16.4%)
---------------------------------------------------------------------------------------------------------------Lo mismo puede hacerse con gráficos más elaborados (vista en
explorador)
view(dfSummary(antro))
Pruebas de hipótesis
Importante: no utilizar el test correcto va en contra del investigador ya que supone no encontrar diferencias cuando realmente existen. Siempre que sea posiblea, utilizar variables continuas, tienen más potencia.
Prueba de hipótesis para la media, μ, de una población normal
t.test(x, y = NULL, alternative = c(“two.sided”, “less”, “greater”), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, …)
Los argumentos a definir dentro de t.test para hacer la prueba son:
- x: vector numérico con los datos.
- alternative: tipo de hipótesis alterna. Los valores disponibles son “two.sided” cuando la hipótesis alterna es ≠, “less” para el caso < y “greater” para >
- mu: valor de referencia de la prueba.
- conf.level: nivel de confianza para reportar el intervalo de confianza asociado (opcional).
Test de normalidad
Para utilizar el t.test es preciso que la variable tenga una distribución normal (como se ha mencionado, podría utilizarse también en caso contrario pero el test pierde potencia por lo que es preferible utilizar otro si no es así)
Para ver si las variablers tienen una distribución normal, pueden utilizarse alguno de estos test incluidos en la librería nortest
library(nortest)Test de Kolmogorov-Smirnov
lillie.test(antro$Hb)
Lilliefors (Kolmogorov-Smirnov) normality test
data: antro$Hb
D = 0.061679, p-value < 2.2e-16Test de Anderson-Darling
ad.test(antro$Hb)
Anderson-Darling normality test
data: antro$Hb
A = 124.21, p-value < 2.2e-16También puede utilizarse el test de Shapiro si el tamaño de la muestra es menor de 5000
shapiro.test(antro$Hb)
Podemos hacer un histograma, en el que también se ve que no es una distribución normal
hist(antro$Hb)
Como no es normal, habría que utilizar un test no paramétrico. Si
utilizamos el t.test y queremos ver si la media de Hb es menor de 10 en
la población, sería así:
H0: Hb=10
H1: Hb>10
t.test(antro$Hb,
alternative="greater", mu=10)
One Sample t-test
data: antro$Hb
t = 175.08, df = 10542, p-value < 2.2e-16
alternative hypothesis: true mean is greater than 10
95 percent confidence interval:
12.45301 Inf
sample estimates:
mean of x
12.47627 En este caso, como p<0.05, rechazamos la hipótesis nula, es decir, el nivel de Hb de la población es mayor de 10.
Pruebas de hipótesis para la proporción, p, de una población
Prueba χ2 de Pearson
prop.test(x, n, p = NULL, alternative = c(“two.sided”, “less”, “greater”), conf.level = 0.95, correct = TRUE)
Los argumentos a definir dentro de prop.test para hacer la prueba son:
- x: número de éxitos en la muestra.
- n: número de observaciones en la muestra.
- alternative: tipo de hipótesis alterna. Los valores disponibles son “two.sided” cuando la alterna es ≠, “less” para el caso < y “greater” para >
- p: valor de referencia de la prueba.
- correct: valor lógico para indicar si se usa la corrección de Yates.
- conf.level: nivel de confianza para reportar el intervalo de confianza asociado (opcional).
Ejemplo, supongamos que queremos saber si en nuestro estudio la proporción de obesidad es del 12.5% ya que ese es el porcentaje que hay en España y se considera el de referencia
Entonces, H0: p=0.125, H1: p≠0.125
n sería el número de filas (pacientes) y x el número de obesos
nrow(antro)[1] 10551table(antro$rango)
Normopeso Obesidad Sobrepeso
7178 1368 2005 prop.test(x=1368, n=10551, p=0.125)
1-sample proportions test with continuity correction
data: 1368 out of 10551, null probability 0.125
X-squared = 2.0488, df = 1, p-value = 0.1523
alternative hypothesis: true p is not equal to 0.125
95 percent confidence interval:
0.1233343 0.1362491
sample estimates:
p
0.129656 Concluimos que no tenemos evidencias para afirmar que el porcentaje de obesos en nuestro estudio sea distinto al de España (0.125)
Prueba binomial exacta
Se hace cuando tenemos un tamaño muestral n, menor de 30
binom.test(x, n, p = 0.5, alternative = c(“two.sided”, “less”, “greater”), conf.level = 0.95)
Los argumentos a definir dentro de binom.test para hacer la prueba son:
- x: número de éxitos en la muestra.
- n: número de observaciones en la muestra.
- alternative: tipo de hipótesis alterna. Los valores disponibles son “two.sided” cuando la alterna es ≠, “less” para el caso < y “greater” para >
- p: valor de referencia de la prueba.
- conf.level: nivel de confianza para reportar el intervalo de confianza asociado (opcional).
binom.test(x=1371, n=10576, p=0.125)
Exact binomial test
data: 1371 and 10576
number of successes = 1371, number of trials = 10576, p-value = 0.1497
alternative hypothesis: true probability of success is not equal to 0.125
95 percent confidence interval:
0.1232867 0.1361842
sample estimates:
probability of success
0.1296331 Prueba de hipótesis para la razón de varianzas σ21/ σ22
Cuando se comparan dos medias, además de comprobar que la variable es normal, debemos comprobar que la varianza de esa variable en ambos grupos es la misma: homogeneidad de varianzas
var.test(creat ~ sexo, data=antro)
F test to compare two variances
data: creat by sexo
F = 0.37417, num df = 5075, denom df = 5466, p-value < 2.2e-16
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.3544995 0.3949593
sample estimates:
ratio of variances
0.3741727 Hay otros test para comprobar la homogeneidad de varianzas, el test de Barlett
bartlett.test(antro$creat~antro$sexo)
Bartlett test of homogeneity of variances
data: antro$creat by antro$sexo
Bartlett's K-squared = 1209.4, df = 1, p-value < 2.2e-16O el test de Levene (necesita la librería car)
library(car)
leveneTest(antro$creat~antro$sexo)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 1674.7 < 2.2e-16 ***
10541
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Y como H0 : σ21 / σ22 = 1, vemos que las varianzas no son homogéneas
Prueba de hipótesis para la igualdad de medias μA = μB
t.test(x, y = NULL, alternative = c(“two.sided”, “less”, “greater”), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, …)
Los argumentos a definir dentro de t.test para hacer la prueba son:
- x: vector numérico con la información de la muestra 1,
- y: vector numérico con la información de la muestra 2,
- alternative: tipo de hipótesis alterna. Los valores disponibles son “two.sided” cuando la alterna es ≠, “less” para el caso < y “greater” para >
- mu: valor de referencia de la prueba (opcinal, no necesario la mayoría de veces ya que siempre queremos comparar si son iguales lo que implica que la diferencia es 0 que es el valor por defecto).
- var.equal=TRUE: indica que las varianzas son desconocidas pero iguales. Si no lo son, basta con poner var.equal=FALSE.
- conf.level: nivel de confianza para reportar el intervalo de confianza asociado (opcional).
t.test(creat ~ sexo, data=antro,
var.equal=FALSE)
Welch Two Sample t-test
data: creat by sexo
t = -234.07, df = 9157.4, p-value < 2.2e-16
alternative hypothesis: true difference in means between group mujer and group varon is not equal to 0
95 percent confidence interval:
-1.700171 -1.671931
sample estimates:
mean in group mujer mean in group varon
0.5331954 2.2192464 El test asume que las varianzas son iguales, así que si no son iguales, es necesario poner var.equal=FALSE ya que por defecto ese argumento es TRUE
Prueba de hipótesis para la igualdad de dos o más medias (ANOVA)
H0: μA = μB = ⋯ =
μK
H1: algún par de medias es distinta
mod <- aov(ip ~ zona, data=antro)
summary(mod) Df Sum Sq Mean Sq F value Pr(>F)
zona 3 14775 4925 1.364e+29 <2e-16 ***
Residuals 10547 0 0
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A continuación, hay que ver entre qué grupos hay diferencias, y para ello se utilizan los post-hoc tests.
TukeyHSD(mod) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = ip ~ zona, data = antro)
$zona
diff lwr upr p adj
Zona2-Zona1 1.33 1.33 1.33 0
Zona3-Zona1 -2.13 -2.13 -2.13 0
Zona4-Zona1 -1.47 -1.47 -1.47 0
Zona3-Zona2 -3.46 -3.46 -3.46 0
Zona4-Zona2 -2.80 -2.80 -2.80 0
Zona4-Zona3 0.66 0.66 0.66 0Y vemos que hay diferencias entre todos los grupos.
Prueba de hipótesis para la igualdad de proporciones pA = pB
H0: pA = pB
H0: pA ≠ pB
El test de Chi-cuadrado que se calcula a partir de la tabla de contingencia, añadiendo chisq=TRUE
library(summarytools)ctable(antro$sexo, antro$rango,
useNA="no", prop="c", chisq=TRUE)Cross-Tabulation, Column Proportions
sexo * rango
Data Frame: antro
------- ------- --------------- --------------- --------------- ----------------
rango Normopeso Obesidad Sobrepeso Total
sexo
mujer 3534 ( 49.2%) 559 ( 40.9%) 988 ( 49.3%) 5081 ( 48.2%)
varon 3644 ( 50.8%) 809 ( 59.1%) 1017 ( 50.7%) 5470 ( 51.8%)
Total 7178 (100.0%) 1368 (100.0%) 2005 (100.0%) 10551 (100.0%)
------- ------- --------------- --------------- --------------- ----------------
----------------------------
Chi.squared df p.value
------------- ---- ---------
33.496 2 0
----------------------------O bien
chisq.test(antro$sexo, antro$rango)
Pearson's Chi-squared test
data: antro$sexo and antro$rango
X-squared = 33.496, df = 2, p-value = 5.327e-08Pruebas no paramétricas
Test de Mann-Whitney-Wilcoxon
Para la media de una población (cuando no es normal), H0: μ = 10
wilcox.test(antro$Hb, mu=10)
Wilcoxon signed rank test with continuity correction
data: antro$Hb
V = 55435185, p-value < 2.2e-16
alternative hypothesis: true location is not equal to 10Para la diferencia de medias, H0: μA = μB
wilcox.test(creat~sexo, data=antro)
Wilcoxon rank sum test with continuity correction
data: creat by sexo
W = 0, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0Test de Kruskall-Wallis
La alternativa a un test de ANOVA (este requiere normalidad), es utilizar el test de Kruskall-Wallis
kruskal.test(ip~zona, data=antro)
Kruskal-Wallis rank sum test
data: ip by zona
Kruskal-Wallis chi-squared = 10550, df = 3, p-value < 2.2e-16Al igual que con ANOVA, debemos realizar tests a posteriori para ver entre qué grupos hay differencias. Para el caso no paramétrico, lo podemos hacer con el test de Nemenyi de la librería DescTools
library(DescTools)Y al igual que en el caso del ANOVA, usaremos el método propuesto por Tukey. La variable grupal (argumento g) debe ser factor (también podemos convertir a factor esa variable antes de hacer el test)
pvals <- NemenyiTest(x = antro$ip,
g = as.factor(antro$zona), dist="tukey")
pvals
Nemenyi's test of multiple comparisons for independent samples (tukey)
mean.rank.diff pval
Zona2-Zona1 1643.0 <2e-16 ***
Zona3-Zona1 -7131.0 <2e-16 ***
Zona4-Zona1 -3498.5 <2e-16 ***
Zona3-Zona2 -8774.0 <2e-16 ***
Zona4-Zona2 -5141.5 <2e-16 ***
Zona4-Zona3 3632.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Importante: Los test paramétricos suelen funcionar bien con muchos casos. Los test no paramétricos muestran su mejor potencia cuando las bases de datos son pequeñas.
Test de Fisher
Es la versión no paramétrica de Chi2. Este test es necesario no sólo cuando hay poco tamaño muestral, sino también cuando hay pocos casos en algunas de las celdas de la tabla.
fisher.test(antro$sexo, antro$rango)
Fisher's Exact Test for Count Data
data: antro$sexo and antro$rango
p-value = 4.848e-08
alternative hypothesis: two.sidedModelos de regresión
La regresión lineal, es el método de regresión paramétrica que usamos cuando la variable de resultado o respuesta es continua. Cuando el resultado es binario, utilizamos la regresión logística.
Regresión lineal
Un modelo lineal, funcion lm, es paramétrico porque asumimos que la relación entre dos variables es lineal y puede ser definida por los parámetros de una recta.
En los modelos de regresión lineal (1 variable predictora) y de regresión lineal multivariante hay un intercept (β0), una o más variables predictoras (x1, x2…) y un error (e):
y= β0 + β1 x1 + β2 x2+ .. .+ e
En la base de datos mtcars, veremos si el consumo de combustible “mpg” está correlacionado con el peso del coche “wt”
data(mtcars)
simple_model <- lm(mpg~wt, data=mtcars)
summary(simple_model)
Call:
lm(formula = mpg ~ wt, data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-4.5432 -2.3647 -0.1252 1.4096 6.8727
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
wt -5.3445 0.5591 -9.559 1.29e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10Adjusted R-squared es 0,74, lo que indica que wt explica el 74% de la variabilidad de mpg
Añadimos otra variable predictora, una versión binaria de caballos de fuerza (hp), 0 para valores inferiores al promedio de hp y 1 para valores superiores al promedio
multivariable_model <- lm(mpg ~ wt + hp , data = mtcars)
summary(multivariable_model)
Call:
lm(formula = mpg ~ wt + hp, data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-3.941 -1.600 -0.182 1.050 5.854
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.22727 1.59879 23.285 < 2e-16 ***
wt -3.87783 0.63273 -6.129 1.12e-06 ***
hp -0.03177 0.00903 -3.519 0.00145 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.593 on 29 degrees of freedom
Multiple R-squared: 0.8268, Adjusted R-squared: 0.8148
F-statistic: 69.21 on 2 and 29 DF, p-value: 9.109e-12Vemos que ambas variables son significativas para explicar el consumo del coche.
Los resultados de la regresión solo son precisos si se dan un conjunto de supuestos (en orden de importancia):
1- Validez de los datos para responder a la pregunta de investigación. 2- Linealidad de la relación entre el resultado y las variables predictoras. 3- Independencia de los errores (en particular, sin correlación entre errores consecutivos como en el caso de los datos de series de tiempo). 4- Varianza igual de errores (homocedasticidad). 5- Normalidad de errores.
Mediante los siguientes gráficos podemos determinar si podemos usar nuestro modelo o no para realizar predicciones - Residuals vs fitted: en esta gráfica los residuals tienen que verse homogéneamente distribuidos a lo largo de la recta entre -2 y +2 y la línea tiene que ser recta, no curva. - Normal Q-Q: los puntos tienen que estar en la recta - Las otras 2 gráficas sirven para ver si hay outliers
par(mfrow=c(2,2))
plot(multivariable_model)
Y ya podemos ver que el modelo no es muy adecuado
Regresión logística
Los parámetros del modelo pueden estimarse con la función glm (modelo lineal generalizado) indicando que la familia que estamos modelando es la binomial. La variable resultado en este caso es binaria (0/1). Esto es importante, ya que si no indicamos nada, esta función glm realizará una estimación de los parámetros asumiendo que nuestra variable resultado es continua (regresión lineal).
La variable resultado (y) es el logaritmo de (y/1-y), es decir el log de la odds, que es la probabilidad de que ocurra (y) frente a que no ocurra (1-y).
Usaremos un ejemplo comparando el grupo de normopeso y exceso de peso (sobrepeso + obesidad)
table(antro$rango)
Normopeso Obesidad Sobrepeso
7178 1368 2005 Es necesario recodificar la variable y es aconsejable tener la información como 0/1 en lugar de “normopeso”, “exceso de peso”
antro$status <- ifelse(antro$rango=="Normopeso", 0, 1)
table(antro$status)
0 1
7178 3373 modelo_simple <- glm(status ~ puntostest, data=antro, family=binomial)
summary(modelo_simple)
Call:
glm(formula = status ~ puntostest, family = binomial, data = antro)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.3952 -0.9032 -0.8122 1.4032 1.6698
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.193688 0.042779 -27.90 <2e-16 ***
puntostest 0.084616 0.006484 13.05 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 12734 on 10069 degrees of freedom
Residual deviance: 12562 on 10068 degrees of freedom
(481 observations deleted due to missingness)
AIC: 12566
Number of Fisher Scoring iterations: 4Vemos que la puntuacion del test está asociada a tener exceso de peso. El coeficiente para el test es 0.085. La OR es el exp(0.085) = 0.915 por lo que la interpretación es que hay un 8.5% más de probabilidad (riesgo) de tener exceso de peso por cada punto que aumenta el test.
Introduzcamos en el modelo la variable antecedente familiar (variable “familiar”)
modelo_multivariable <- glm(status ~ puntostest + familia, data=antro, family=binomial)
summary(modelo_multivariable)
Call:
glm(formula = status ~ puntostest + familia, family = binomial,
data = antro)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9881 -0.2708 -0.2423 0.6888 2.7031
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.68359 0.10223 -36.033 < 2e-16 ***
puntostest 0.05657 0.01117 5.066 4.07e-07 ***
familiapositivo 4.54899 0.08887 51.188 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 10800.3 on 8461 degrees of freedom
Residual deviance: 5196.4 on 8459 degrees of freedom
(2089 observations deleted due to missingness)
AIC: 5202.4
Number of Fisher Scoring iterations: 6Ambas variables son estadísticamente significativas porque el p-valor es <0.05 en ambas.
¿Cuál de estos dos modelos es mejor? Para esta pregunta no usamos el criterio de información de Akaike (AIC). AIC menor indica mejor modelo.
AIC(modelo_simple)[1] 12566.32AIC(modelo_multivariable)[1] 5202.397Es mejor el modelo multivariable.
Creación de modelos
La función step() en R base automatiza la selección de variables paso a paso usando AIC. Primero tenemos que definir nuestro modelo completo, es decir, el modelo con las variables que queremos usar para la selección de variables. En nuestro caso supongamos que queremos ver qué variables son importantes entre Hb, creat, edad1sex, edad1emb, edaddiag, edadmadre, niveledu, tabaco, ets, acos, puntostest y familia.
La función step da error si hay datos con NA, para eso hay que seleccionar casos con datos completos:
antro_complete<- antro[complete.cases(antro),]El modelo sería entonces:
mod <- glm(status ~ Hb + creat + edad1sex + edad1emb + edaddiag + edadmadre + niveledu + tabaco + ets + acos + puntostest + familia,
data = antro_complete, family="binomial")
modF <- step(mod, trace = F, direction = "backward")
summary(modF)
Call:
glm(formula = status ~ Hb + creat + edad1sex + edaddiag + edadmadre +
niveledu + tabaco + ets + familia, family = "binomial", data = antro_complete)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.5653 -0.3003 -0.1721 0.4360 3.2250
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.292402 0.668733 -6.419 1.37e-10 ***
Hb 0.089222 0.042833 2.083 0.037250 *
creat 0.319190 0.066400 4.807 1.53e-06 ***
edad1sex -0.053895 0.014091 -3.825 0.000131 ***
edaddiag 0.080645 0.006559 12.294 < 2e-16 ***
edadmadre -0.052017 0.007636 -6.812 9.62e-12 ***
niveleduprimaria -0.085594 0.159127 -0.538 0.590649
niveledusecundaria -0.517196 0.196838 -2.628 0.008601 **
niveleduticnico -1.086820 0.369938 -2.938 0.003305 **
niveleduuniversitario -1.572987 0.452104 -3.479 0.000503 ***
tabacofumador 0.585288 0.249226 2.348 0.018853 *
tabacono fumador -0.092012 0.192289 -0.479 0.632286
etssi 0.352714 0.120391 2.930 0.003393 **
familiapositivo 4.347322 0.132130 32.902 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 5071.8 on 3944 degrees of freedom
Residual deviance: 2182.5 on 3931 degrees of freedom
AIC: 2210.5
Number of Fisher Scoring iterations: 6Para evitar escribir todas las variables se ppne . y puede quitarse alguna variable del modelo con “-”
mod <- glm(status ~ . -puntostest,
data = antro_complete, family="binomial")
summary(mod)
Call:
glm(formula = status ~ . - puntostest, family = "binomial", data = antro_complete)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.409e-06 -2.409e-06 -2.409e-06 2.409e-06 2.409e-06
Coefficients: (1 not defined because of singularities)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.657e+01 1.100e+06 0 1
zonaZona2 -6.430e-10 2.095e+04 0 1
zonaZona3 -3.007e-10 2.176e+04 0 1
zonaZona4 1.027e-09 1.764e+04 0 1
fechanac -2.025e-11 4.405e+01 0 1
sexovaron 1.911e-08 1.138e+05 0 1
peso -5.673e-10 1.584e+03 0 1
talla 3.274e-10 1.324e+03 0 1
fechavisit -5.617e-12 2.029e+01 0 1
ip NA NA NA NA
rangoObesidad 5.313e+01 4.380e+05 0 1
rangoSobrepeso 5.313e+01 4.377e+05 0 1
ciap_apNO -3.718e-09 2.062e+05 0 1
ciap_apSI 2.557e-07 2.097e+05 0 1
ciap_imcNO -2.732e-09 3.567e+05 0 1
ciap_imcSI 3.951e-08 2.526e+05 0 1
Hb -1.500e-08 3.760e+04 0 1
creat -1.137e-08 6.680e+04 0 1
edad1sex 5.710e-12 2.389e+03 0 1
edad1emb -5.668e-11 2.152e+03 0 1
edaddiag 8.747e-11 6.882e+02 0 1
edadmadre -2.206e-11 7.994e+02 0 1
niveleduprimaria 5.313e-09 1.665e+04 0 1
niveledusecundaria 4.480e-09 2.028e+04 0 1
niveleduticnico 2.774e-09 3.664e+04 0 1
niveleduuniversitario 5.470e-09 4.186e+04 0 1
tabacofumador 6.184e-10 2.570e+04 0 1
tabacono fumador -2.244e-09 2.023e+04 0 1
etssi 9.578e-10 1.255e+04 0 1
acossi -5.327e-10 1.304e+04 0 1
familiapositivo -1.968e-09 1.792e+04 0 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 5.0718e+03 on 3944 degrees of freedom
Residual deviance: 2.2887e-08 on 3915 degrees of freedom
AIC: 60
Number of Fisher Scoring iterations: 25modF <- step(mod, trace = F, direction = "backward")
summary(modF)
Call:
glm(formula = status ~ rango, family = "binomial", data = antro_complete)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.409e-06 -2.409e-06 -2.409e-06 2.409e-06 2.409e-06
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -26.57 6993.59 -0.004 0.997
rangoObesidad 53.13 16558.25 0.003 0.997
rangoSobrepeso 53.13 14479.33 0.004 0.997
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 5.0718e+03 on 3944 degrees of freedom
Residual deviance: 2.2887e-08 on 3942 degrees of freedom
AIC: 6
Number of Fisher Scoring iterations: 25trace=F surve para que no se muestren los resultados. La dirfección por defecto es “both”
También podemos filtrar una base de datos para hacer unos análisis específicos en un subgrupo de individuos. Para ello se utiliza la función filter() o subset. EJ:
familia.positivo <- filter(antro, status=="1" & familia=="positivo")
dim(familia.positivo)[1] 2726 23Nomogramas
Ver https://www.r-bloggers.com/2017/04/generating-dynamic-nomograms-using-dynnom/
Para hacer nomogramas las variables deben ser factor o numéricas y esto puede hacerse al importar los datos (read.delim (stringsAsfactors=T)) o bien modificando el tipo de variable tras la importación
antro <- read.csv2("~/Documents/TRABAJO/Curso R/antro.csv", stringsAsFactors=TRUE)
str(antro)'data.frame': 10551 obs. of 22 variables:
$ zona : Factor w/ 4 levels "Zona1","Zona2",..: 4 4 4 4 4 4 4 4 4 4 ...
$ fechanac : Factor w/ 3960 levels "","1/1/04","1/1/07",..: 50 444 2652 335 335 3053 80 1896 2955 3928 ...
$ sexo : Factor w/ 2 levels "mujer","varon": 1 1 1 1 1 1 1 1 1 1 ...
$ peso : num 50.5 33.9 27.5 33.4 33.9 19.5 45.4 17.7 51.9 50.1 ...
$ talla : num 156 150 132 148 149 ...
$ fechavisit: Factor w/ 1391 levels "","1/10/13","1/10/14",..: 346 345 1277 30 30 1288 980 852 1251 807 ...
$ ip : num -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 -1.67 ...
$ rango : Factor w/ 3 levels "Normopeso","Obesidad",..: 1 1 1 1 1 1 1 1 1 1 ...
$ ciap_ap : Factor w/ 3 levels "","NO","SI": 3 2 2 2 2 2 2 2 2 2 ...
$ ciap_imc : Factor w/ 3 levels "","NO","SI": 2 2 2 2 2 2 2 2 2 2 ...
$ Hb : num 15 15 15 15 15 ...
$ creat : num 1 1 1 1 1 1 1 1 1 1 ...
$ edad1sex : int 16 14 23 29 22 25 33 36 18 17 ...
$ edad1emb : int 16 14 24 29 22 28 33 36 18 17 ...
$ edaddiag : int 40 NA 30 45 42 NA NA NA 23 45 ...
$ edadmadre : int 64 48 41 51 58 57 64 56 41 52 ...
$ niveledu : Factor w/ 5 levels "ninguno","primaria",..: 2 1 2 2 2 2 1 2 2 1 ...
$ tabaco : Factor w/ 3 levels "ex-fumador","fumador",..: 1 2 3 3 3 1 1 3 2 1 ...
$ ets : Factor w/ 2 levels "no","si": 2 1 1 1 1 1 1 1 1 2 ...
$ acos : Factor w/ 2 levels "no","si": 1 1 1 2 1 2 1 1 1 1 ...
$ puntostest: int 9 11 3 1 8 5 2 6 3 13 ...
$ familia : Factor w/ 2 levels "negativo","positivo": 1 2 1 1 1 1 1 1 1 1 ...Hay que cargar estas librerías
library(magrittr)
library(DynNom)Crear el nomograma
DynNom(modF, antro_complete)
Creación de tablas descriptivas
Se realiza con la librería compareGroups
library(compareGroups)Ver documentación: https://cran.r-project.org/web/packages/compareGroups/vignettes/compareGroups_vignette.html
La librería compareGroups tiene 3 funciones clave:
compareGroups(): genera los cálculos
createTable(): crea la tabla descriptiva creada con compareGroups (). Se puede customizar excluyendo categorías de las variables, cambiando el número de decimales, etc.
export2…(): exporta las tablas a EXCEL, Word, LaTeX, Rmarkdown, etc.
Para los ejemplos se utiliza el dataframe incluido “predimed” (grupo control, grupo con nueces, grupo con aceite de oliva)
data(predimed)Comparar los grupos y generar los cálculos
Se hace con la función compareGroups()
Las variables descritas se colocan en al lado derecho
de ~ separadas por el signo +,
mientras que la variable que indica los grupos se
coloca en el lado izquierdo de la fórmula. Para seleccionar todas las
variables podemos usar ‘.’ y para eliminar variables,
‘-’.
No se permiten transformaciones en el entorno de fórmulas. Si es necesario, deben realizarse antes de llamar a la función compareGroups()
Cada variable es descrita en función del tipo (numérica o categórica). El argumento method fuerza a que los análisis de una variable se haga asumiendo que la variable es:
method = 1 normalmente distribuida;
method = 2 continua no normal;
method = 3 categórica;
method = NA, lleva a cabo un test de normalidad y decide si la variable
es normal o no (valor por defecto).
Así, si quisiéramos describir todas (.) las variables de nuestros datos menos la variable toevent (- toevent) y que el programa decidiera qué variables son normales o no (method = NA) excepto para la variable wtn y p14 que queremos que se reporten con la mediana y los cuartiles en vez de la media y desviación estandard, ejecutaríamos
res <- compareGroups(group ~ . - toevent, data = predimed, method = c(wtn = 2, p14 = 2))
res
-------- Summary of results by groups of 'Intervention group'---------
var N p.value method
1 Sex 6324 <0.001** categorical
2 Age 6324 0.003** continuous normal
3 Smoking 6324 0.444 categorical
4 Body mass index 6324 <0.001** continuous normal
5 Waist circumference 6324 0.045** continuous normal
6 Waist-to-height ratio 6324 <0.001** continuous normal
7 Hypertension 6324 0.249 categorical
8 Type-2 diabetes 6324 0.017** categorical
9 Dyslipidemia 6324 0.423 categorical
10 Family history of premature CHD 6324 0.581 categorical
11 Hormone-replacement therapy 5661 0.850 categorical
12 MeDiet Adherence score 6324 <0.001** continuous non-normal
13 AMI, stroke, or CV Death 6324 0.064* categorical
selection
1 ALL
2 ALL
3 ALL
4 ALL
5 ALL
6 ALL
7 ALL
8 ALL
9 ALL
10 ALL
11 ALL
12 ALL
13 ALL
-----
Signif. codes: 0 '**' 0.05 '*' 0.1 ' ' 1 En esta tabla vemos el número de individuos sin valores faltantes (missings) para cada variable, el p-valor correspondiente al test de normalidad y el tipo de variable que considera para el análisis.
Si queremos ver las descriptivas, basta con usar la función summary ().
Los símbolos [ ] se utilizan para seleccionar variables. Podemos seleccionarlas con la posición 1:2 o con el nombre c(“sex”, “age”)
summary(res, c(“sex”, “age”))
summary(res[1:3])
--- Descriptives of each row-variable by groups of 'Intervention group' ---
-------------------
row-variable: Sex
Male Female Male% Female% p.overall p.trend
[ALL] 2679 3645 42.36243 57.63757
Control 812 1230 39.76494 60.23506 8.1e-05 0.388386
MedDiet + Nuts 968 1132 46.09524 53.90476
MedDiet + VOO 899 1283 41.20073 58.79927
p.Control vs MedDiet + Nuts p.Control vs MedDiet + VOO
[ALL]
Control 0.000133 0.358324
MedDiet + Nuts
MedDiet + VOO
p.MedDiet + Nuts vs MedDiet + VOO
[ALL]
Control 0.002076
MedDiet + Nuts
MedDiet + VOO
-------------------
row-variable: Age
N mean sd lower upper p.overall p.trend
[ALL] 6324 67.0117 6.17499 66.85948 67.16392
Control 2042 67.34231 6.27992 67.06977 67.61485 0.002666 0.101163
MedDiet + Nuts 2100 66.6819 6.016395 66.42444 66.93937
MedDiet + VOO 2182 67.01971 6.212578 66.75889 67.28052
p.Control vs MedDiet + Nuts p.Control vs MedDiet + VOO
[ALL]
Control 0.001672 0.20596
MedDiet + Nuts
MedDiet + VOO
p.MedDiet + Nuts vs MedDiet + VOO
[ALL]
Control 0.172672
MedDiet + Nuts
MedDiet + VOO
-------------------
row-variable: Smoking
Never Current Former Never% Current% Former% p.overall
[ALL] 3892 858 1574 61.54333 13.56736 24.88931
Control 1282 270 490 62.78159 13.22233 23.99608 0.444354
MedDiet + Nuts 1259 296 545 59.95238 14.09524 25.95238
MedDiet + VOO 1351 292 539 61.91567 13.38222 24.70211
p.trend p.Control vs MedDiet + Nuts p.Control vs MedDiet + VOO
[ALL]
Control 0.572674 0.517748 0.834021
MedDiet + Nuts
MedDiet + VOO
p.MedDiet + Nuts vs MedDiet + VOO
[ALL]
Control 0.630418
MedDiet + Nuts
MedDiet + VOO Crear tabla descriptiva
Podemos crear la tabla descriptiva aplicando la función createTable() al objeto que hemos obtenido tras usar la función compareGroups() (res).
Ver documentación en
https://cran.r-project.org/web/packages/compareGroups/vignettes/compareGroups_vignette.html
Argumentos de createTable():
- type: personaliza cómo se muestran las variables categóricas. 1 porcentajes, 3 frecuencias absolutas y 2 o NA ambas (opción por defecto)
- hide.no = “no”: oculta el nivel “no” para las variables binarias.
- hide: para una varible categórica, indica qué categoría se oculta. También se aplica a variables categóricas con más de dos categorías.
- digits: especifica el número de dígitos decimales que se mostrarán
- show.n = TRUE: muestra la N total de cada variable
- show.descr = FALSE solo se muestran los p-values
- show.all = TRUE muestra una columna con descriptivos de todos los datos
- show.p.overall = FALSE omite p-values
- show.p.trend = TRUE: si la variable respuesta tiene más de 2 categorías se muestra un p-value para la tendencia
restab <- createTable(res, digits = c(p14 = 0, hormo=1), type = 2,
hide = c(sex = "Male"), hide.no = "no", show.n = TRUE)
restab
--------Summary descriptives table by 'Intervention group'---------
________________________________________________________________________________________
Control MedDiet + Nuts MedDiet + VOO p.overall N
N=2042 N=2100 N=2182
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Sex: Female 1230 (60.2%) 1132 (53.9%) 1283 (58.8%) <0.001 6324
Age 67.3 (6.28) 66.7 (6.02) 67.0 (6.21) 0.003 6324
Smoking: 0.444 6324
Never 1282 (62.8%) 1259 (60.0%) 1351 (61.9%)
Current 270 (13.2%) 296 (14.1%) 292 (13.4%)
Former 490 (24.0%) 545 (26.0%) 539 (24.7%)
Body mass index 30.3 (3.96) 29.7 (3.77) 29.9 (3.71) <0.001 6324
Waist circumference 101 (10.8) 100 (10.6) 100 (10.4) 0.045 6324
Waist-to-height ratio 0.63 (0.07) 0.62 (0.06) 0.63 (0.06) <0.001 6324
Hypertension 1711 (83.8%) 1738 (82.8%) 1786 (81.9%) 0.249 6324
Type-2 diabetes 970 (47.5%) 950 (45.2%) 1082 (49.6%) 0.017 6324
Dyslipidemia 1479 (72.4%) 1539 (73.3%) 1560 (71.5%) 0.423 6324
Family history of premature CHD 462 (22.6%) 460 (21.9%) 507 (23.2%) 0.581 6324
Hormone-replacement therapy 31 (1.7%) 30 (1.6%) 36 (1.8%) 0.850 5661
MeDiet Adherence score 8 [7;10] 9 [8;10] 9 [8;10] <0.001 6324
AMI, stroke, or CV Death 97 (4.75%) 70 (3.33%) 85 (3.90%) 0.064 6324
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ header.labels: para cambiar cabeceras
Por ejemplo, para cambiar “p.overall” por “p-value”
print(restab, header.labels = c(p.overall = "p-value"))
--------Summary descriptives table by 'Intervention group'---------
______________________________________________________________________________________
Control MedDiet + Nuts MedDiet + VOO p-value N
N=2042 N=2100 N=2182
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Sex: Female 1230 (60.2%) 1132 (53.9%) 1283 (58.8%) <0.001 6324
Age 67.3 (6.28) 66.7 (6.02) 67.0 (6.21) 0.003 6324
Smoking: 0.444 6324
Never 1282 (62.8%) 1259 (60.0%) 1351 (61.9%)
Current 270 (13.2%) 296 (14.1%) 292 (13.4%)
Former 490 (24.0%) 545 (26.0%) 539 (24.7%)
Body mass index 30.3 (3.96) 29.7 (3.77) 29.9 (3.71) <0.001 6324
Waist circumference 101 (10.8) 100 (10.6) 100 (10.4) 0.045 6324
Waist-to-height ratio 0.63 (0.07) 0.62 (0.06) 0.63 (0.06) <0.001 6324
Hypertension 1711 (83.8%) 1738 (82.8%) 1786 (81.9%) 0.249 6324
Type-2 diabetes 970 (47.5%) 950 (45.2%) 1082 (49.6%) 0.017 6324
Dyslipidemia 1479 (72.4%) 1539 (73.3%) 1560 (71.5%) 0.423 6324
Family history of premature CHD 462 (22.6%) 460 (21.9%) 507 (23.2%) 0.581 6324
Hormone-replacement therapy 31 (1.7%) 30 (1.6%) 36 (1.8%) 0.850 5661
MeDiet Adherence score 8 [7;10] 9 [8;10] 9 [8;10] <0.001 6324
AMI, stroke, or CV Death 97 (4.75%) 70 (3.33%) 85 (3.90%) 0.064 6324
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ También podemos crear tablas descriptivas que no estén separadas por una variable grupal. Basta con dejar el lado izquierdo de ~ vacío.
resNoGroups <- compareGroups(~ . , predimed)
restabNoGroups <- createTable(resNoGroups, hide.no = "no")
print(restabNoGroups, header.labels = c("all" = "Entire cohort"))
--------Summary descriptives table ---------
__________________________________________________
Entire cohort N
N=6324
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Intervention group: 6324
Control 2042 (32.3%)
MedDiet + Nuts 2100 (33.2%)
MedDiet + VOO 2182 (34.5%)
Sex: 6324
Male 2679 (42.4%)
Female 3645 (57.6%)
Age 67.0 (6.17) 6324
Smoking: 6324
Never 3892 (61.5%)
Current 858 (13.6%)
Former 1574 (24.9%)
Body mass index 30.0 (3.82) 6324
Waist circumference 100 (10.6) 6324
Waist-to-height ratio 0.63 (0.07) 6324
Hypertension 5235 (82.8%) 6324
Type-2 diabetes 3002 (47.5%) 6324
Dyslipidemia 4578 (72.4%) 6324
Family history of premature CHD 1429 (22.6%) 6324
Hormone-replacement therapy 97 (1.71%) 5661
MeDiet Adherence score 8.68 (1.94) 6324
follow-up to main event (years) 4.36 (1.69) 6324
AMI, stroke, or CV Death 252 (3.98%) 6324
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ Visualización
Finalmente, con compareGroups también podemos resumir los datos de forma gráfica mediante la función genérica plot (). Por ejemplo, podemos visualizar las dos primeras variables con
plot(res[1:2])
También podemos ver la distribución según la variable grupal
plot(res[1:2], bivar = TRUE)
Modelos estadísticos (OR y HR)
Odds ratio
show.ratio =TRUE: si la variable respuesta es dicotómica o se ha definido como class survival, pueden mostrarse en los resultados las Odds Ratios y Hazard Ratios
Para los estudios de casos y controles, puede ser interesante calcular la odds ratio (OR) para cada variable entre casos y controles. Aunque el estudio PREDIMED no es un caso-control, usaremos variable event como si fuese una variable caso-control. La forma más sencilla de obtener una tabla similar a la que hay en los artículos científicos es:
resOR <- compareGroups(event ~ . - toevent, predimed)
restabOR <- createTable(resOR, show.ratio = TRUE,
show.p.overall = FALSE)
restabOR
--------Summary descriptives table by 'AMI, stroke, or CV Death'---------
__________________________________________________________________________________
No Yes OR p.ratio
N=6072 N=252
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Intervention group:
Control 1945 (32.0%) 97 (38.5%) Ref. Ref.
MedDiet + Nuts 2030 (33.4%) 70 (27.8%) 0.69 [0.50;0.95] 0.021
MedDiet + VOO 2097 (34.5%) 85 (33.7%) 0.81 [0.60;1.09] 0.173
Sex:
Male 2528 (41.6%) 151 (59.9%) Ref. Ref.
Female 3544 (58.4%) 101 (40.1%) 0.48 [0.37;0.62] <0.001
Age 66.9 (6.14) 69.4 (6.65) 1.07 [1.04;1.09] <0.001
Smoking:
Never 3778 (62.2%) 114 (45.2%) Ref. Ref.
Current 809 (13.3%) 49 (19.4%) 2.01 [1.41;2.82] <0.001
Former 1485 (24.5%) 89 (35.3%) 1.99 [1.49;2.64] <0.001
Body mass index 30.0 (3.81) 29.8 (3.92) 0.98 [0.95;1.02] 0.365
Waist circumference 100 (10.6) 102 (10.6) 1.01 [1.00;1.03] 0.016
Waist-to-height ratio 0.63 (0.07) 0.63 (0.07) 3.64 [0.55;23.9] 0.178
Hypertension:
No 1047 (17.2%) 42 (16.7%) Ref. Ref.
Yes 5025 (82.8%) 210 (83.3%) 1.04 [0.75;1.48] 0.826
Type-2 diabetes:
No 3231 (53.2%) 91 (36.1%) Ref. Ref.
Yes 2841 (46.8%) 161 (63.9%) 2.01 [1.55;2.62] <0.001
Dyslipidemia:
No 1645 (27.1%) 101 (40.1%) Ref. Ref.
Yes 4427 (72.9%) 151 (59.9%) 0.56 [0.43;0.72] <0.001
Family history of premature CHD:
No 4694 (77.3%) 201 (79.8%) Ref. Ref.
Yes 1378 (22.7%) 51 (20.2%) 0.87 [0.63;1.18] 0.363
Hormone-replacement therapy:
No 5341 (98.2%) 223 (99.6%) Ref. Ref.
Yes 96 (1.77%) 1 (0.45%) 0.29 [0.01;1.27] 0.117
MeDiet Adherence score 8.70 (1.94) 8.24 (1.94) 0.89 [0.84;0.95] <0.001
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ Esta función se puede personalizar. Por ejemplo el argumento fact.ratio nos permite determinar la interpretación de la OR según el incremento de la variable independiente (fact.ratio = c(waist=10)). También podemos cambiar las cabeceras de la tabla como hemos explicado anteriormente. Para más posibilidades ver los argumentos en ?compareGroups.
resOR <- compareGroups(event ~ . - toevent, predimed,
fact.ratio = c(waist=10))
restabOR <- createTable(resOR, show.ratio = TRUE,
show.p.overall = FALSE,
hide.no = "no",
hide = c(sex = "Male"), type=1)
print(restabOR, header.labels = c(p.ratio = "p-value"))
--------Summary descriptives table by 'AMI, stroke, or CV Death'---------
________________________________________________________________________________
No Yes OR p-value
N=6072 N=252
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Intervention group:
Control 32.0% 38.5% Ref. Ref.
MedDiet + Nuts 33.4% 27.8% 0.69 [0.50;0.95] 0.021
MedDiet + VOO 34.5% 33.7% 0.81 [0.60;1.09] 0.173
Sex: Female 58.4% 40.1% 0.48 [0.37;0.62] <0.001
Age 66.9 (6.14) 69.4 (6.65) 1.07 [1.04;1.09] <0.001
Smoking:
Never 62.2% 45.2% Ref. Ref.
Current 13.3% 19.4% 2.01 [1.41;2.82] <0.001
Former 24.5% 35.3% 1.99 [1.49;2.64] <0.001
Body mass index 30.0 (3.81) 29.8 (3.92) 0.98 [0.95;1.02] 0.365
Waist circumference 100 (10.6) 102 (10.6) 1.15 [1.03;1.30] 0.016
Waist-to-height ratio 0.63 (0.07) 0.63 (0.07) 3.64 [0.55;23.9] 0.178
Hypertension 82.8% 83.3% 1.04 [0.75;1.48] 0.826
Type-2 diabetes 46.8% 63.9% 2.01 [1.55;2.62] <0.001
Dyslipidemia 72.9% 59.9% 0.56 [0.43;0.72] <0.001
Family history of premature CHD 22.7% 20.2% 0.87 [0.63;1.18] 0.363
Hormone-replacement therapy 1.77% 0.45% 0.29 [0.01;1.27] 0.117
MeDiet Adherence score 8.70 (1.94) 8.24 (1.94) 0.89 [0.84;0.95] <0.001
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ Hazard ratio
Cuando analizamos estudios de cohortes, solemos estar interesados en calcular los hazard ratio (HR) en vez de OR ya que de esta manera tenemos en cuenta el tiempo transcurrido hasta el evento y los posibles valores censurados (análisis de supervivencia).
Para llevar a cabo estos análisis con compareGroups primero debemos calcular la variable tiempo hasta el evento como una variable de supervivencia mediante la función Surv de la librería survival (que no hace falta instalar porque ya está en R).
library(tidyverse)
library(survival)
predimed <- mutate(predimed,
eventSurv = Surv(toevent, event == "Yes"))Tras esto, podemos usar esta variable a lado izquierdo de ~ y usar una sintaxis similar a la que usamos para calcualr OR simplemente cambiando event por eventSurv.
resHZ <- compareGroups(eventSurv ~ . - toevent, predimed)
restabHZ <- createTable(resOR, show.ratio = TRUE,
show.p.overall = FALSE)
restabHZ
--------Summary descriptives table by 'AMI, stroke, or CV Death'---------
__________________________________________________________________________________
No Yes OR p.ratio
N=6072 N=252
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Intervention group:
Control 1945 (32.0%) 97 (38.5%) Ref. Ref.
MedDiet + Nuts 2030 (33.4%) 70 (27.8%) 0.69 [0.50;0.95] 0.021
MedDiet + VOO 2097 (34.5%) 85 (33.7%) 0.81 [0.60;1.09] 0.173
Sex:
Male 2528 (41.6%) 151 (59.9%) Ref. Ref.
Female 3544 (58.4%) 101 (40.1%) 0.48 [0.37;0.62] <0.001
Age 66.9 (6.14) 69.4 (6.65) 1.07 [1.04;1.09] <0.001
Smoking:
Never 3778 (62.2%) 114 (45.2%) Ref. Ref.
Current 809 (13.3%) 49 (19.4%) 2.01 [1.41;2.82] <0.001
Former 1485 (24.5%) 89 (35.3%) 1.99 [1.49;2.64] <0.001
Body mass index 30.0 (3.81) 29.8 (3.92) 0.98 [0.95;1.02] 0.365
Waist circumference 100 (10.6) 102 (10.6) 1.15 [1.03;1.30] 0.016
Waist-to-height ratio 0.63 (0.07) 0.63 (0.07) 3.64 [0.55;23.9] 0.178
Hypertension:
No 1047 (17.2%) 42 (16.7%) Ref. Ref.
Yes 5025 (82.8%) 210 (83.3%) 1.04 [0.75;1.48] 0.826
Type-2 diabetes:
No 3231 (53.2%) 91 (36.1%) Ref. Ref.
Yes 2841 (46.8%) 161 (63.9%) 2.01 [1.55;2.62] <0.001
Dyslipidemia:
No 1645 (27.1%) 101 (40.1%) Ref. Ref.
Yes 4427 (72.9%) 151 (59.9%) 0.56 [0.43;0.72] <0.001
Family history of premature CHD:
No 4694 (77.3%) 201 (79.8%) Ref. Ref.
Yes 1378 (22.7%) 51 (20.2%) 0.87 [0.63;1.18] 0.363
Hormone-replacement therapy:
No 5341 (98.2%) 223 (99.6%) Ref. Ref.
Yes 96 (1.77%) 1 (0.45%) 0.29 [0.01;1.27] 0.117
MeDiet Adherence score 8.70 (1.94) 8.24 (1.94) 0.89 [0.84;0.95] <0.001
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ Exportar las tablas
Las tablas creadas se pueden exportar a Word (export2word ()), CSV (export2csv ()) y otros formatos (pdf, LaTeX, Mardown). Si usamos export2md () lo exportaremos a Markdown y si ponemos esto en un documento de R Mardown la visualización de la tabla será mucho mejor. Por ejemplo:
por defecto
export2md(restabOR)| No | Yes | OR | p.ratio | |
|---|---|---|---|---|
| N=6072 | N=252 | |||
| Intervention group: | ||||
| Control | 32.0% | 38.5% | Ref. | Ref. |
| MedDiet + Nuts | 33.4% | 27.8% | 0.69 [0.50;0.95] | 0.021 |
| MedDiet + VOO | 34.5% | 33.7% | 0.81 [0.60;1.09] | 0.173 |
| Sex: Female | 58.4% | 40.1% | 0.48 [0.37;0.62] | <0.001 |
| Age | 66.9 (6.14) | 69.4 (6.65) | 1.07 [1.04;1.09] | <0.001 |
| Smoking: | ||||
| Never | 62.2% | 45.2% | Ref. | Ref. |
| Current | 13.3% | 19.4% | 2.01 [1.41;2.82] | <0.001 |
| Former | 24.5% | 35.3% | 1.99 [1.49;2.64] | <0.001 |
| Body mass index | 30.0 (3.81) | 29.8 (3.92) | 0.98 [0.95;1.02] | 0.365 |
| Waist circumference | 100 (10.6) | 102 (10.6) | 1.15 [1.03;1.30] | 0.016 |
| Waist-to-height ratio | 0.63 (0.07) | 0.63 (0.07) | 3.64 [0.55;23.9] | 0.178 |
| Hypertension | 82.8% | 83.3% | 1.04 [0.75;1.48] | 0.826 |
| Type-2 diabetes | 46.8% | 63.9% | 2.01 [1.55;2.62] | <0.001 |
| Dyslipidemia | 72.9% | 59.9% | 0.56 [0.43;0.72] | <0.001 |
| Family history of premature CHD | 22.7% | 20.2% | 0.87 [0.63;1.18] | 0.363 |
| Hormone-replacement therapy | 1.77% | 0.45% | 0.29 [0.01;1.27] | 0.117 |
| MeDiet Adherence score | 8.70 (1.94) | 8.24 (1.94) | 0.89 [0.84;0.95] | <0.001 |
o bien coloreada por filas de variables
export2md(restabOR, strip = TRUE, first.strip = TRUE)| No | Yes | OR | p.ratio | |
|---|---|---|---|---|
| N=6072 | N=252 | |||
| Intervention group: | ||||
| Control | 32.0% | 38.5% | Ref. | Ref. |
| MedDiet + Nuts | 33.4% | 27.8% | 0.69 [0.50;0.95] | 0.021 |
| MedDiet + VOO | 34.5% | 33.7% | 0.81 [0.60;1.09] | 0.173 |
| Sex: Female | 58.4% | 40.1% | 0.48 [0.37;0.62] | <0.001 |
| Age | 66.9 (6.14) | 69.4 (6.65) | 1.07 [1.04;1.09] | <0.001 |
| Smoking: | ||||
| Never | 62.2% | 45.2% | Ref. | Ref. |
| Current | 13.3% | 19.4% | 2.01 [1.41;2.82] | <0.001 |
| Former | 24.5% | 35.3% | 1.99 [1.49;2.64] | <0.001 |
| Body mass index | 30.0 (3.81) | 29.8 (3.92) | 0.98 [0.95;1.02] | 0.365 |
| Waist circumference | 100 (10.6) | 102 (10.6) | 1.15 [1.03;1.30] | 0.016 |
| Waist-to-height ratio | 0.63 (0.07) | 0.63 (0.07) | 3.64 [0.55;23.9] | 0.178 |
| Hypertension | 82.8% | 83.3% | 1.04 [0.75;1.48] | 0.826 |
| Type-2 diabetes | 46.8% | 63.9% | 2.01 [1.55;2.62] | <0.001 |
| Dyslipidemia | 72.9% | 59.9% | 0.56 [0.43;0.72] | <0.001 |
| Family history of premature CHD | 22.7% | 20.2% | 0.87 [0.63;1.18] | 0.363 |
| Hormone-replacement therapy | 1.77% | 0.45% | 0.29 [0.01;1.27] | 0.117 |
| MeDiet Adherence score | 8.70 (1.94) | 8.24 (1.94) | 0.89 [0.84;0.95] | <0.001 |
o bien cambiar el tamaño
export2md(restabOR, size=6)| No | Yes | OR | p.ratio | |
|---|---|---|---|---|
| N=6072 | N=252 | |||
| Intervention group: | ||||
| Control | 32.0% | 38.5% | Ref. | Ref. |
| MedDiet + Nuts | 33.4% | 27.8% | 0.69 [0.50;0.95] | 0.021 |
| MedDiet + VOO | 34.5% | 33.7% | 0.81 [0.60;1.09] | 0.173 |
| Sex: Female | 58.4% | 40.1% | 0.48 [0.37;0.62] | <0.001 |
| Age | 66.9 (6.14) | 69.4 (6.65) | 1.07 [1.04;1.09] | <0.001 |
| Smoking: | ||||
| Never | 62.2% | 45.2% | Ref. | Ref. |
| Current | 13.3% | 19.4% | 2.01 [1.41;2.82] | <0.001 |
| Former | 24.5% | 35.3% | 1.99 [1.49;2.64] | <0.001 |
| Body mass index | 30.0 (3.81) | 29.8 (3.92) | 0.98 [0.95;1.02] | 0.365 |
| Waist circumference | 100 (10.6) | 102 (10.6) | 1.15 [1.03;1.30] | 0.016 |
| Waist-to-height ratio | 0.63 (0.07) | 0.63 (0.07) | 3.64 [0.55;23.9] | 0.178 |
| Hypertension | 82.8% | 83.3% | 1.04 [0.75;1.48] | 0.826 |
| Type-2 diabetes | 46.8% | 63.9% | 2.01 [1.55;2.62] | <0.001 |
| Dyslipidemia | 72.9% | 59.9% | 0.56 [0.43;0.72] | <0.001 |
| Family history of premature CHD | 22.7% | 20.2% | 0.87 [0.63;1.18] | 0.363 |
| Hormone-replacement therapy | 1.77% | 0.45% | 0.29 [0.01;1.27] | 0.117 |
| MeDiet Adherence score | 8.70 (1.94) | 8.24 (1.94) | 0.89 [0.84;0.95] | <0.001 |
o en formato word
export2word(restabOR, file="table.docx")o a pdf, (requiere TEX)
export2pdf(restabOR, file="table.pdf")
Gráficos en R
Gráficos para una variable cuantitativa
Gráfico de tallo y hoja
antro <- read.csv2("~/Documents/TRABAJO/Curso R/antro.csv")
stem(antro$talla)
The decimal point is 1 digit(s) to the right of the |
8 | 0466677888888888999999999999
9 | 00000000000000000001111111111111111111111111111122222222222222222222+744
10 | 00000000000000000000000000000000000000000000000000000000000000000000+1872
11 | 00000000000000000000000000000000000000000000000000000000000000000000+2522
12 | 00000000000000000000000000000000000000000000000000000000000000000000+1421
13 | 00000000000000000000000000000000000000000000000000000000000000000000+1208
14 | 00000000000000000000000000000000000000000000000000000000000000000000+934
15 | 00000000000000000000000000000000000000000000000000000000000000000000+850
16 | 00000000000000000000000000000000000000000000000000000001111111111111+270
17 | 0000000000111111111222233333333444455566679999
18 | 00389Boxplot
La estructura de la función boxplot con los argumentos más comunes de uso se muestran a continuación:
boxplot(x, …)
- x: vector numérico con los datos para crear el
boxplot.
- formula: fórmula con la estructura x ~ g para
indicar que las observaciones en el vector x van a ser agrupadas de
acuerdo a los niveles del factor g.
- data: base de datos (data frame) con las
variables.
- range: valor numérico que indica la extensión de
los bigotes. Si es positivo, los bigotes se extenderán hasta el punto
más extremo de tal manera que el bigote no supere veces el rango
intercuatílico (IQ). Un valor de cero hace que los bigotes se extiendan
hasta los datos extremos.
- col: vector con los colores a usar en el cuerpo de
las cajas.
- log: para indicar si las coordenadas x o y o serán graficadas en escala logarítmica.
En R usamos de forma indistinta “ o
’ para las etiquetas.
En R, si queremos ver dos gráficas en una misma figura, podemos usar la
función par() con el argumento mfrow.
par(mfrow=c(1,2))
boxplot(x=antro$peso, ylab='Peso', col = 'blue')
boxplot(x=antro$talla, xlab='Talla', horizontal=TRUE, col = 'red')
Para volver a la visualización de un gráfico
par(mfrow=c(1,1))Los colores pueden verse con la función colors(). Los 30 primeros son:
head(colors(), n=30) [1] "white" "aliceblue" "antiquewhite" "antiquewhite1"
[5] "antiquewhite2" "antiquewhite3" "antiquewhite4" "aquamarine"
[9] "aquamarine1" "aquamarine2" "aquamarine3" "aquamarine4"
[13] "azure" "azure1" "azure2" "azure3"
[17] "azure4" "beige" "bisque" "bisque1"
[21] "bisque2" "bisque3" "bisque4" "black"
[25] "blanchedalmond" "blue" "blue1" "blue2"
[29] "blue3" "blue4" Es posible crear boxplots para comparar dos o varios grupos definidos por 1 o 2 variables cualitativas
par(mfrow=c(1, 2))
boxplot(talla ~ rango, data=antro,
col=c('lightblue', 'pink'),
xlab='Rango', main='A',
ylab='Talla')
boxplot(talla ~ rango*zona, data=antro,
col=c('lightblue', 'pink'),
xlab='Rango y Zona', main='B',
ylab='Talla')
Histograma
par(mfrow=c(1, 1))
hist(antro$talla, col='blue', main='',
ylab='Talla')
density(x, …)
- x: es el vector con los datos para los cuales se
quiere la densidad.
- bw: ancho de banda (se usa para suavizar más o
menos la gráfica).
- kernel: núcleo de suavización a usar, los posibles
valores son gaussian, rectangular, triangular, epanechnikov, biweight,
cosine o optcosine, el valor por defecto es gaussian. Esto es muy
avanzado y no suele cambiarse.
- na.rm: valor lógico, si es TRUE se eliminan los valores con NA para construir la densidad, el valor por defecto es FALSE.
plot(density(antro$peso, na.rm=TRUE))
El siguiente código muestra cómo afecta el ancho de banda en la
estimación.
Estos valores de suavizado dependen de la escala que tengamos en nuestra
variable de interés
par(mfrow=c(2,2))
plot(density(antro$peso, na.rm=TRUE), main='')
plot(density(antro$peso, bw=0.01, na.rm=TRUE), main='')
plot(density(antro$peso, bw=5, na.rm=TRUE), main='')
plot(density(antro$peso, bw=20, na.rm=TRUE), main='')
También podemos obtener la misma distribución separada para casos y controles o grupos (aumentamos el limite del eje y con ylim=c(0, 0.60))
datos.casos <- subset(antro, zona=='Zona3')
datos.control <- subset(antro, zona=='Zona4')
den.casos <- density(datos.casos$peso, na.rm=TRUE)
den.control <- density(datos.control$peso, na.rm=TRUE)
par(mfrow=c(1,1))
plot(den.casos,
main='Edad peso por zona', ylab='Densidad',
xlab='Zona (caso/control)', lwd=4, col='blue', ylim=c(0, 0.08))
lines(den.control, lwd=4, col='red')
legend('topright', legend=c('Caso', 'Control'), bty='n',
lwd=3, col=c('blue', 'red'))
Gráficos para dos variables cuantitativas
plot(x, y, main, sub, xlab, ylab, …)
- x: vector numérico con las coordenadas del eje
horizontal.
- y: vector numérico con las coordenadas del eje
vertical.
- type: tipo de gráfico a crear Las opciones
son:
- ‘p’ para obtener puntos, esta es la opción por
defecto.
- ‘l’ para obtener líneas.
- ‘b’ para obtener los puntos y líneas que unen los
puntos.
- ‘c’ para obtener sólo las líneas y dejando los
espacios donde estaban los puntos obtenidos con la opción ‘b’.
- ‘o’ para obtener los puntos y líneas
superpuestas.
- ‘h’ para obtener líneas verticales desde el origen
hasta el valor yi de cada punto, similar a un histograma.
- ‘s’ para obtener escalones.
- ‘S’ similar al anterior.
- ‘n’ para que no dibuje (útil para hacer gráficos ad
hoc).
- …: otros parámetros gráficos (ver ?plot).
- main: título del gráfico
- sub: subtítulo del gráfico
- xlab: etiqueta eje X
- ylab: etiqueta eje Y
library("MASS")Esta librería contiene una base de datos de coches
plot(Cars93$Weight, Cars93$EngineSize,
ylab="Engine Size",
xlab="Weight",
main="My plot")
Podemos colorear (col) y cambiar la forma de los puntos (pch) en función de otras variables y añadir líneas. ‘lty’ define el típo de línea, ‘lwd’ el grosor y ‘col’ el color
plot(Cars93$Weight, Cars93$EngineSize,
col=Cars93$Type,
pch=as.numeric(Cars93$Type))
lines(x=c(min(Cars93$Weight), max(Cars93$Weight)),
y=c(min(Cars93$EngineSize),
max(Cars93$EngineSize)), lwd=4, lty=3, col="blue")
abline(h=3, lty=1, lwd= 3, col='green')
abline(v=2000, lty=6, lwd=4, col='red')
También podemos cambiar un solo punto:
plot(Cars93$Weight, Cars93$EngineSize, ylab="Engine Size",
xlab="Weight", main="My plot")
points(x=min(Cars93$Weight), y=min(Cars93$EngineSize),
pch=16, col="red")
points(x=max(Cars93$Weight), y=max(Cars93$EngineSize),
pch=16, col="blue")
Y podemos añadir texto con las funciones:
text() (se añade en el gráfico)
mtext() (se añade en el margen del gráfico). side
indica 1: abajo, 2: izquierda, 3: arriba y 4: derecha).
legend() añade una leyenda dentro del gráfico en la
posición que queramos según una de estas opciones: c(“bottomright”,
“bottom”, “bottomleft”, “left”, “topleft”, “top”, “topright”, “right”,
“center”).
plot(Cars93$Weight, Cars93$EngineSize, ylab="Engine Size",
xlab="Weight", main="My plot")
text(x=2000, y=5, "texto con coordenadas")
mtext(side=3, "subtítulo", line=0.45)
legend("bottomright", legend=c("Data Points"), pch="o")
Podemos añadir la recta de regresión mediante:
plot(Cars93$Weight, Cars93$EngineSize, ylab="Engine Size",
xlab="Weight", main="My plot")
mod <- lm(EngineSize ~ Weight, data=Cars93)
abline(mod, lwd=3, col="blue")
Si queremos visualizar relaciones no lineales, podemos usar una estimación suavizada mediante regresión local o loess:
plot(Cars93$EngineSize ~ Cars93$Weight)
lines(loess.smooth(Cars93$Weight, Cars93$EngineSize), col='blue')
La función pairs() crea un panel con todas los gráficos de dispersión entre todas las variables de la base de datos que le pasemos. Por ejemplo, para las primeras 7 variables:
pairs(Cars93[,1:7])
Podemos eliminar los diagramas de dispersión bajo la diagonal (lower.panel) porque son figuras redundantes, editar los nombres de las variables (labels), darles mayor tamaño a la letra de equiquetas (cex.labels), usando puntos (pch=1) o cruces (pch=3), cambiarles el color (col), el tamaño (cex), cambiar la escala a logaritmica (log) y poner marcas horizontales en el eje vertical (las=1) o vertical (las=2). El símbolo barra invertida n es para introducir un salto de línea en la etiqueta.
pairs(Cars93[,1:7], lower.panel=NULL,
cex.labels=1.4, log='xy',
main='Matriz de dispersión', las=1,
labels=c('Fabricante', 'Modelo', 'Tipo', 'Precio \n mínimo', 'Precio' , 'Precio \n máximo', 'Ciudad'),
pch=3, cex=0.4, col='dodgerblue2')
Gráficos para variables cualitativas
Gráficos de barras
barplot() se usa para obtener un gráfico de barras. La estructura de la función barplot con los argumentos más comunes de uso se muestra a continuación:
- height: vector o matriz con la información de las
frecuencias absolutas o relativas.
- beside: valor lógico para indicar si las barras
deben estar pegadas (útil cuando la información es una matriz)
- horiz: valor lógico para indicar si el diagrama de
barras debe ser horizontal, por defecto es FALSE.
- las: para rotar 90º las etiquetas de las categorías y mejorar la visualización en caso de ser muy largas.
Para representar cuántos individuos hay en cada zona, necesitamos tener esta información en modo de tabla.
table(antro$zona)
Zona1 Zona2 Zona3 Zona4
1426 1860 1694 5571 info <- prop.table(table(antro$zona))
barplot(info, las=2)
Añadimos ylim=c(0, 0.60) para aumentar el eje Y
Luego se usa la función text() para incluir un texto en
las coordenadas x=xx y y=info, el parámetro pos=3 coloca el texto en la
parte superior de las barras y el parámetro label sirve para indicar lo
que se desea escribir en las coordenadas indicadas, en este caso son las
frecuencias relativas, redondeadas a 2 decimales, almacenadas en
info.
xx <- barplot(info, col='darkblue', ylim=c(0, 0.60),
xlab='', main='Zonas de estudio ',
ylab='Frecuencia relativa', las=2)
text(x=xx, y=info, pos=3, cex=0.8, col="red",
label=round(info, 2))
Supongamos que ahora queremos obtener la misma información por sexos
table(antro$sexo, antro$zona)
Zona1 Zona2 Zona3 Zona4
mujer 712 878 814 2677
varon 714 982 880 2894info2 <- table(antro$sexo, antro$zona)
barplot(info2, las=2) 
o bien,
barplot(info2, beside=TRUE, las=2)
Lo mejoramos
barplot(info2, beside = TRUE, las=1,
xlab ='', ylab ='Frecuencia',
ylim = c(0,3000),
main = 'Individuos por zona y sexo',
col = c("mistyrose","lightblue"),)
legend('topleft', legend = rownames(info2),
bty='n', fill=c("mistyrose","lightblue"))
Gráficos circulares
pie(x, labels)
- x: vector con elementos no negativos que
representan las frecuencias de los niveles de la variable
cualitativa.
- labels: vector con los nombres a visualizar en cada parte del pastel, por defecto se usan los nombres del vector x.
pie(info, main="Individuos por zona")
Gráficos de puntos
Los gráficos de puntos son muy útiles para representar tablas de frecuencias (de 1 o 2 vías) o tablas de resumen en relación a una o dos variables. La estructura de la función dotchart () se muestra a continuación.
dotchart(x, labels=NULL, groups=NULL, gdata=NULL, pt.cex, pch, color, lcolor, …)
- x: vector o matriz con la información de las
frecuencias o medida de resumen a representar. Si x es una matriz las
columnas representarán agrupaciones.
- labels: vector con los nombres a usar para los
puntos, por defecto toma los nombres de las filas de la matriz x.
- groups: vector con los nombres a usar para los
grupos, por defecto toma los nombres de las columnas de la matriz
x.
- cex: tamaño de los nombres a mostrar en los
ejes.
- pt.cex: tamaño del punto.
- pch: tipo de punto a usar (igual que en
points).
- color: tipo de color usar para los puntos.
- lcolor: color para la línea asociada a cada
punto.
- …: otros parámetros gráficos que pueden ser pasados como argumentos.
Supongamos que queremos graficar el número de individuos estudiados en cada zona
dotchart(x = info,
cex=1.2, xlab="Individuos del estudio por zona")
Librería LATTICE
Esta librería es muy potente ya que permite hacer muchos de los
gráficos que hemos visto hasta ahora, de forma sencilla.
También es muy útil para crear gráficos condicionales. Es decir,
observar cómo cambia y con x en los niveles de z. Estos tipos de
gráficos son útiles para analizar datos multidimensionales.
library(lattice)Panel de figuras que nos describa la relación entre la variable MPG.city y Price separada por cada tipo de coche (Type).
xyplot(MPG.city ~ EngineSize | Type, data=Cars93)
Mismo gráfico separado por tipo de coche y origen,
xyplot(MPG.city ~ EngineSize | Type * Origin, data = Cars93)
Si queremos obtener el mismo gráfico con la información de cada tipo superpuesta en cada panel, usamos panel = panel.superpose y el argumento groups para separar por origen
xyplot(MPG.city ~ EngineSize|Type,
panel=panel.superpose, groups=Origin,
data = Cars93)
Podemos añadir una leyenda con key
pars <- trellis.par.get("superpose.symbol")
xyplot(MPG.city ~ EngineSize | Type,
panel = panel.superpose, groups=Origin,
key = list(columns = 2,
text = list(levels(Cars93$Origin)),
points = Rows(pars,1:2)),
data=Cars93)
Existen otras funciones en la librería lattice que hacen funciones similares para otros tipos de gráfico:
- splom( ~ data.frame) # Scatterplot matrix
- bwplot(factor ~ numeric, …) # Boxplot
- qqmath(factor ~ numeric, …) # Q-Q plot
- dotplot(factor ~ numeric, …) # 1-D display
- stripplot(factor ~ numeric, …)
- barchar(character ~ numeric, …)
- histogram( ~ numeric, …)
- densityplot( ~ numeric, …) # Smoothed version of histogram
Por ejemplo, podríamos hacer un boxplot para el peso según la zona y rango:
bwplot(peso ~ rango | zona,
data = antro)
ggplot2
Con ggplot2 pueden hacerse gráficos avanzados mucho más versátiles y vistosos, pero requiere más aprendizaje.
En esta web con ejemplos, puede buscarse el gráfico deseado y copiar el código para modificarlo y adaptarlo a nuestras preferencias.
La estructura básica del código para hacer un gráfico en ggplot es:
plot = data + aesthetics + geometry
y para facilitar la visualización suele ponerse en líneas con un + al final:
ggplot(data + # nuestra base de datos
aes(x= , y= , …) + # estética de los datos
geom_tipo() # define el tipo de gráfico
La estética (aes) de los datos se refiere a: - posición (en los ejes) - color exterior (color) y de relleno (fill) - forma de puntos (shape) - tipo de línea (linetype) - tamaño (size)
Y la geometría (geom) al tipo de gráfico (puntos,
líneas, etc.):
- geom_point (puntos)
- geom_lines (lineas)
- geom_histogram (histograma)
- geom_boxplot (boxplot)
- geom_bar (barras)
- geom_smooth (líneas suavizadas) - geom_polygons (polígonos en un mapa)
- …
Con el comando help.search(“geom_,” package = “ggplot2”) se puede ver el listado de objetos geométricos.
Por ejemplo para hacer una gráfica de puntos de la base de datos iris que represente la relación de longitud de sépalos y pétalos:
library(ggplot2)
ggplot(iris) +
aes(x = Sepal.Length, y=Petal.Length) +
geom_point()
También se puede escribir así:
ggplot(iris) +
geom_point(mapping = aes(x = Sepal.Length, y = Petal.Length))
Se puede añadir un tema (theme) que da un formato estético determinado - theme_dark() - theme_minimal() - theme_classic() - theme_void() - theme_test()
ggplot(iris) +
aes(x = Sepal.Length, y=Petal.Length) +
geom_point() +
theme_minimal()
Para variables categóricas
ggplot(data = antro) +
geom_bar(mapping = aes(x = zona))
Para variables continuas
ggplot(data = iris) +
geom_histogram(mapping = aes(x = Sepal.Width), binwidth = 0.1)
Variable continua según una variable categórica
ggplot(data = iris, mapping = aes(x = Sepal.Width)) +
geom_freqpoly(mapping = aes(colour = Species), binwidth = 0.1)
Mejor con un boxplot y la variable continua en el eje y:
ggplot(data = iris, mapping = aes(y = Sepal.Width, x=Species)) +
geom_boxplot(mapping = aes(colour = Species))
Podemos mejorarlo:
ggplot(iris, aes(x=Species, y=Sepal.Width) ) +
geom_boxplot(alpha=0.8, outlier.colour = "darkblue") +
geom_point(stat= "summary", fun.y=median,
shape=15, size=2.5, color="red")
Si todo este código lo guardamos en un objeto, podremos añadir más código de manera más sencilla:
box <- ggplot(iris, aes(x=Species, y=Sepal.Width) ) +
geom_boxplot(alpha=0.8, outlier.colour = "darkblue") +
geom_point(stat= "summary", fun.y=median,
shape=15, size=2.5, color="red")Añadimos código para girar las coordenadas:
box +
coord_flip()
O añadirle más geometrías:
box +
geom_jitter(width = 0.1, alpha = 0.2)
Para dos variables continuas
ggplot(antro) +
geom_point(mapping = aes(x = peso, y = talla)) +
theme_minimal()
Cuando hay muchos puntos y se superponen, podemos usar la estética alpha:
ggplot(antro) +
geom_point(mapping = aes(x = peso, y = talla),
alpha = 1 / 100)
ggplot requiere más aprendizaje pero pueden hacerse gráficos muy interesantes, como el de este otro ejemplo, al que se le han añadido la media y la mediana con ‘stat_summary’:
library(ggplot2)
ggplot(antro, aes(x=zona, y=peso, fill=zona, color=zona)) +
geom_violin(width=1, size=1) +
geom_boxplot(width=0.1, color="black", alpha=0.2) +
theme_minimal() +
theme(legend.position="none") +
xlab("Zona de salud") +
ylab("Peso (kg)")+
stat_summary(fun = "median",
geom = "crossbar",
width=0.1,
colour = "blue")+
stat_summary(fun = "mean",
geom = "point",
colour = "red")+
scale_y_continuous(breaks = seq(0,130,by = 10))
Patchwork es una librería para colocar varios
gráficos juntos de forma sencilla
https://patchwork.data-imaginist.com/index.html
library(patchwork)Con esta librería, las gráficas pueden unirse sencillamente, por ejemplo:
G1 + G2 + G3 G1 + G2 / G3 G1 | G2 + G3
R Markdown
Este documento que estás leyendo es un documento de R Markdown. Markdown es una sintaxis de formato simple para crear documentos HTML, PDF y MS Word entre otros.
Para más información sobre el uso de R Markdown ver http://rmarkdown.rstudio.com en donde está esta guía completa https://bookdown.org/yihui/rmarkdown/
Crear un fichero Markdown
En la barra de menú, hacer clic en File -> New File -> R Markdown… o simplemente hacer click en el signo más verde en la esquina superior izquierda de RStudio.
A continuación hay que elegir el documento que quiere crearse, html, doc o pdf…
Una vez creado el documento, al hacer clic en el botón Knit, se generará el documento que incluirá tanto su texto como la salida de cualquier fragmento de código R incrustado dentro del documento.
Cabecera YAML
YAML es una estructura de lista anidada que incluye metadatos del documento. Está escrito entre dos líneas de tres guiones — y RStudio lo escribe automáticamente cuando se crea el documento con la barra de menus.
title: “RMarkdown”
author: “Juan José Lasarte Velillas”
date: “06 April, 2023”
output: html_document
Se puede agregar una tabla de contenidos al documento creado usando la opción toc en el encabezado YAML. Las opciones son:
- toc: true (por defecto false)
- toc_float: true para una tabla de contenidod flotante
- toc_depth: Cuántos niveles de secciones debe incluir la tabla de contenidos
- doc_depth: 3 incluirá secciones (cabeceras) con ###.
la palabra toc debe comenzar en la 5ª columna
Sintaxis
En estos archivos rmd pueden hacerse secciones con una o más almohadillas (#):
ejemplo con 4 almohadillas
ejemplo con 5 almohadillas
ejemplo con 6 almohadillas
y también escribirse:
itálica (entre 1 asterisco a cada lado)
negrilla (entre 2 asteriscos a cada lado)
subíndice (idem con ~)
superíndice (idem con ^)
Para ver más consultar:
https://raw.githubusercontent.com/rstudio/cheatsheets/main/translations/spanish/rmarkdown_es.pdf
Para ecuaciones se utiliza código LaTeX
Código de R en línea
El código en línea se crea usando una tilde inversa (`) y la letra r seguida de otra tilde inversa.
Por ejemplo, si hemos hecho un t.test para ver si hay diferencias significativas en el consumo de combustible (mpg) según las transmisiones sean manuales o automáticas (variable am) (dataframe mtcars). Esto nos dará un p.value
res <- t.test(mpg ~ am, data=mtcars)
res$p.value[1] 0.001373638Y este p.value podemos incluirlo en el texto poniendo, por ejemplo: “el p valor es (aquí incluir ‘`’r res$p.value’´’)”
Bloques (chunks) de código R
En estos bloques se puede insertar código R y se ejecutará automáticamente, evitando así copiar y pegar los resultados. También permite mostrar el código de R utilizado en los cálculos.
Estos chunks pueden insertarse con el botón Insert a new code chunk (verde en la parte superior) o desde el menú Code Insert Chunk o, lo que es más práctico, con las teclas Ctrl+Alt+i (Windows) o Option+Cmd+i (Mac)
Opciones de los bloques de código
En estos chunks pueden añadirse algunas opciones (=TRUE o FALSE) separadas por comas:
- echo, para mostrar o no el código
- warning, para mostrar o no las advertencias
- message, para mostrar o no los mensajes
- error, para mostrar o no los errores
Insertar chunk con el código (echo=TRUE) Especificar echo=T porque en el setup es FALSE
res <- t.test(mpg ~ am, data=mtcars)
res$p.value[1] 0.001373638Insertar chunk sin código (echo=FALSE)
[1] 3.162278Al inicio del documento puede insertarse un chunk con las opciones básicas para que afecte a todo el documento y que podemos nombrar como “setup”. Luego, si es necesario pueden cambiarse para un chunk en particular:
{r setup, include=FALSE, echo=FALSE}
knitr::opts_chunk$set(comment= ““, warning=FALSE, message = FALSE, echo
= FALSE)
También conviene poner nombre a los chunks a continuación de la “r” para que cuando se produzcan errores sepamos dónde están.
Insertar figuras
También puede incrustar gráficos, por ejemplo:
En este chunk se añadió el parámetro echo = FALSE para
evitar la impresión del código R que generó el gráfico.
Opciones del bloque relacionadas con las figuras:
fig.width and fig.height. Por
defecto: fig.width = 7, fig.height = 7
fig.align: Cómo alinear la figura, “left”, “right”, and
“center”
fig.path: directorio donde se quiere guardar la figura
creada por el bloque de código de R. Por defecto, ‘figure/’
dat <- rnorm(100)
boxplot(dat, col='red', xlab='Datos simulados')
Insertar tablas
- kable: Está dentro de la librería knitr, no tiene muchos argumentos, pero produce tablas muy visuales de forma sencilla.
- pander: Está dentro de la librería pander, tiene muchas más opciones y se puede personalizar. Es útil para poner en negrita ciertos valores (por ejemplo, valores por debajo de un umbral).
- xtable: útil para aquellos que trabajen con LaTeX
Es necesario instalar las librerías knitr y pander
Cargamos las librerías:
library(knitr)
library(pander)Cargamos los datos
antro <- read.csv2("~/Documents/TRABAJO/Curso R/antro.csv")Veamos la información que hay en las 10 primeras variables de nuestra base de datos
Con kable
kable(head(antro[,1:10]))| zona | fechanac | sexo | peso | talla | fechavisit | ip | rango | ciap_ap | ciap_imc |
|---|---|---|---|---|---|---|---|---|---|
| Zona4 | 1/2/04 | mujer | 50.5 | 155.5 | 16/2/17 | -1.67 | Normopeso | SI | NO |
| Zona4 | 12/2/04 | mujer | 33.9 | 149.5 | 16/2/16 | -1.67 | Normopeso | NO | NO |
| Zona4 | 28/2/04 | mujer | 27.5 | 132.0 | 7/2/14 | -1.67 | Normopeso | NO | NO |
| Zona4 | 11/4/04 | mujer | 33.4 | 148.0 | 1/7/16 | -1.67 | Normopeso | NO | NO |
| Zona4 | 11/4/04 | mujer | 33.9 | 149.0 | 1/7/16 | -1.67 | Normopeso | NO | NO |
| Zona4 | 30/4/04 | mujer | 19.5 | 110.0 | 7/5/10 | -1.67 | Normopeso | NO | NO |
Con pander
pander (head (antro[ , 1:10]))| zona | fechanac | sexo | peso | talla | fechavisit | ip | rango |
|---|---|---|---|---|---|---|---|
| Zona4 | 1/2/04 | mujer | 50.5 | 155.5 | 16/2/17 | -1.67 | Normopeso |
| Zona4 | 12/2/04 | mujer | 33.9 | 149.5 | 16/2/16 | -1.67 | Normopeso |
| Zona4 | 28/2/04 | mujer | 27.5 | 132 | 7/2/14 | -1.67 | Normopeso |
| Zona4 | 11/4/04 | mujer | 33.4 | 148 | 1/7/16 | -1.67 | Normopeso |
| Zona4 | 11/4/04 | mujer | 33.9 | 149 | 1/7/16 | -1.67 | Normopeso |
| Zona4 | 30/4/04 | mujer | 19.5 | 110 | 7/5/10 | -1.67 | Normopeso |
| ciap_ap | ciap_imc |
|---|---|
| SI | NO |
| NO | NO |
| NO | NO |
| NO | NO |
| NO | NO |
| NO | NO |
library(summarytools)
library(compareGroups)Ahora podemos ver una tabla descriptiva de las variables continuas con pander
pander(stby(antro, INDICES = antro$sex,
FUN = descr, stats = "common", transpose = TRUE))mujer:
Mean Std.Dev Min Median Max N.Valid Pct.Valid creat 0.5332 0.2751 0.06 0.53 1.01 5076 99.9 edad1emb 21.5 4.893 10 21 43 4855 95.55 edad1sex 20.02 4.809 5 19 46 5034 99.07 edaddiag 38.26 12.56 15 36 80 2306 45.38 edadmadre 48.24 11.81 20 48 82 5081 100 Hb 12.48 1.457 10 12.48 15 5077 99.92 ip -1.086 1.177 -2.33 -1.67 1.13 5081 100 peso 27.57 12.55 8.5 23.1 94.8 5074 99.86 puntostest 5.535 3.261 1 5 25 4855 95.55 talla 123 19.39 80 118.5 183 5073 99.84 varon:
Mean Std.Dev Min Median Max N.Valid Pct.Valid creat 2.219 0.4497 1.44 2.21 3 5467 99.95 edad1emb 21.71 4.935 10 21 52 5215 95.34 edad1sex 20.08 4.825 5 19 52 5418 99.05 edaddiag 39.08 13.3 14 37 80 3001 54.86 edadmadre 48.97 12.03 20 49 84 5470 100 Hb 12.47 1.448 10 12.43 15 5466 99.93 ip -1.082 1.189 -2.33 -1.67 1.13 5470 100 peso 27.63 12.75 10.66 23.2 128.5 5465 99.91 puntostest 5.385 3.222 1 5 25 5215 95.34 talla 123.1 18.75 88 119 189 5467 99.95