Bioestadística aplicada a la investigación

install.packages('tidyverse')

Installing package into ‘/cloud/lib/x86_64-pc-linux-gnu-library/4.3’
(as ‘lib’ is unspecified)
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/tidyverse_2.0.0.tar.gz'
Content type 'application/x-gzip' length 425237 bytes (415 KB)
==================================================
downloaded 415 KB

* installing *binary* package ‘tidyverse’ ...
* DONE (tidyverse)

The downloaded source packages are in
    ‘/tmp/Rtmp3bw7Hk/downloaded_packages’

install.packages("rstatix")

Installing package into ‘/cloud/lib/x86_64-pc-linux-gnu-library/4.3’
(as ‘lib’ is unspecified)
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/rstatix_0.7.2.tar.gz'
Content type 'application/x-gzip' length 604777 bytes (590 KB)
==================================================
downloaded 590 KB

* installing *binary* package ‘rstatix’ ...
* DONE (rstatix)

The downloaded source packages are in
    ‘/tmp/Rtmp3bw7Hk/downloaded_packages’

library(tidyverse)

── Attaching core tidyverse packages ────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.4.4     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.0
✔ purrr     1.0.2     ── Conflicts ──────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(rstatix)

Attaching package: ‘rstatix’

The following object is masked from ‘package:stats’:

    filter

Introduccion a R

1. Vectores

Primero agregamos valores de edades

edad = c(23,26,27,28,32,34,28,29,31,25)

Luego agregamos nombres

nombres = c("Maria", "Gustavo", "Ramon", "Silvia", "Alma", "Eduardo", "Enrique", "Erick", "Hugo", "Jesus")

Para revisarlos

nombres

 [1] "Maria"   "Gustavo" "Ramon"   "Silvia"  "Alma"    "Eduardo" "Enrique" "Erick"   "Hugo"   
[10] "Jesus"  

Generación de vectores aleatorios. Genera 10 datos aleatorios en un rango entre 8,000 a 30,000.

salario = sample(8000:30000,size = 10, replace = FALSE)
salario

 [1] 17998  9212 11821 22614 29972 14395 22702 14514 14754 21666

Nota: Si replace = FALSE, entonces no hay reemplazo.

Ahora agregamos 10 datos aleatorios referentes a sexo

sexo = sample(c("H", "M"), size = 10, replace = TRUE)
sexo

 [1] "M" "H" "M" "H" "H" "H" "H" "M" "H" "M"

Tenemos diferentes variables - Cualitativas - Nominales - Ordinales - Cuantitativas - Discretas - Continuas

2. Bases de datos (data frames)

Tenemos bases de datos - Externas - Internas

Este sería un ejemplo de bases de datos internas de R

ejemplo_1 <- datasets::state.x77
head(ejemplo_1)

           Population Income Illiteracy Life Exp Murder HS Grad Frost   Area
Alabama          3615   3624        2.1    69.05   15.1    41.3    20  50708
Alaska            365   6315        1.5    69.31   11.3    66.7   152 566432
Arizona          2212   4530        1.8    70.55    7.8    58.1    15 113417
Arkansas         2110   3378        1.9    70.66   10.1    39.9    65  51945
California      21198   5114        1.1    71.71   10.3    62.6    20 156361
Colorado         2541   4884        0.7    72.06    6.8    63.9   166 103766

La forma de leer una base de datos externa, archivo separado por comas .csv, es con la función read.csv(directorio.csv)

datos1 <- read.csv("/cloud/project/Bioestadistica_R/datos_1.csv")
head(datos1) #Con la función head() se observan los primeros 6 valores del data frame

str(datos1)

'data.frame':   200 obs. of  11 variables:
 $ Clave             : chr  "A001" "A002" "A003" "A004" ...
 $ Sexo              : chr  "M" "H" "M" "H" ...
 $ Edad              : int  21 35 38 18 32 21 39 23 37 35 ...
 $ Estatus           : chr  "Control" "Caso" "Control" "Control" ...
 $ Exposicion        : chr  "Sin exposicion" "Humo de biomasa" "Fumador" "Fumador" ...
 $ Enfermedad        : chr  "Sin alteraciones" "Capacidad pulmonar alterada" "Sin alteraciones" "Sin alteraciones" ...
 $ Cotinina          : num  1 8 183.2 225.5 2.3 ...
 $ Anos.de.exposicion: int  0 5 20 4 0 3 15 5 20 32 ...
 $ Talla             : num  1.65 1.56 1.66 1.77 1.5 1.78 1.8 1.55 1.68 1.7 ...
 $ Peso              : int  88 55 68 104 56 115 98 74 99 65 ...
 $ Padres.fumadores  : chr  "No" "No" "Si" "No" ...

Para cargar la segunda base de datos

datos2 <- read_csv("/cloud/project/Bioestadistica_R/datos_2.csv")     #Aquí usamos la función read_csv del paquete tidyverse

Rows: 200 Columns: 6── Column specification ──────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Clave
dbl (5): Cigarros al dia, Cartuchos al dia, TNFA, CRP, PvCO2
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(datos2)

Creamos una base de datos a partir de los vectores realizados previamente

datos3 <- data.frame(nombres, edad, sexo, salario) 
as_tibble(datos3)

Para cambiar el nombre de una variable utilizamos rename()

rename(datos3, nombre = nombres)   # primero ponemos el (dataframe, nombre nuevo = nombre viejo)

Recordar que ejercicio2 = datos1 y ejercicio3=datos2

Revisar algunas funciones de tidyverse Función select

ejercicio5 <- datos1 %>% 
  select(Clave, Edad, Cotinina, Anos.de.exposicion, Talla, Peso)
ejercicio5

Función filter

ejercicio6 <- datos1 %>% 
  filter(Edad>=25)
ejercicio6

Función filter con variable cualitativa

ejercicio6 <- datos1 %>% 
  filter(Exposicion=="Fumador"|Exposicion=="Vapeador"|Exposicion=="Sin exposicion")
ejercicio6

Ver la función mutate para crear variables nuevas a partir de variables existentes

ejercicio2 <- datos1 %>% 
  mutate(IMC = Peso / Talla^2)
ejercicio2

Podemos cambiar el tipo de variable de cuantitativa a cualitativa con la función mutate

ejercicio2 <- ejercicio2 %>% 
  mutate(IMC_escala = case_when(IMC<18 ~ "Desnutricion",
                                IMC<25 ~ "Normopeso",
                                IMC<30 ~ "Sobrepeso",
                                IMC>=30 ~ "Obesidad"))
ejercicio2

Ahora veremos las rutas de análisis

• Depuración de bases de datos
  • Datos perdidos(na.omit)
  • Filtrar datos (filter)
  • Seleccionar las variables de uso (select)
  • Conversiones / Correcciones (mutate)
  • Analizar la estructura de los datos (str)

str(ejercicio2)    

'data.frame':   200 obs. of  13 variables:
 $ Clave             : chr  "A001" "A002" "A003" "A004" ...
 $ Sexo              : chr  "M" "H" "M" "H" ...
 $ Edad              : int  21 35 38 18 32 21 39 23 37 35 ...
 $ Estatus           : chr  "Control" "Caso" "Control" "Control" ...
 $ Exposicion        : chr  "Sin exposicion" "Humo de biomasa" "Fumador" "Fumador" ...
 $ Enfermedad        : chr  "Sin alteraciones" "Capacidad pulmonar alterada" "Sin alteraciones" "Sin alteraciones" ...
 $ Cotinina          : num  1 8 183.2 225.5 2.3 ...
 $ Anos.de.exposicion: int  0 5 20 4 0 3 15 5 20 32 ...
 $ Talla             : num  1.65 1.56 1.66 1.77 1.5 1.78 1.8 1.55 1.68 1.7 ...
 $ Peso              : int  88 55 68 104 56 115 98 74 99 65 ...
 $ Padres.fumadores  : chr  "No" "No" "Si" "No" ...
 $ IMC               : num  32.3 22.6 24.7 33.2 24.9 ...
 $ IMC_escala        : chr  "Obesidad" "Normopeso" "Normopeso" "Obesidad" ...

La datos que salen es el análisis de estructura

También podemos cambiar el tipo de variable. Por ejemplo, la variable Sexo está categorizada como chr (character), pero la cambiamos a factor. Conviene cambiarlas para realizar adecuadamente el análisis estadístico, ya que las variables tipo chr no son útiles para ese análisis.

ejercicio2$Sexo <- as.factor(ejercicio2$Sexo)
ejercicio2$Estatus <- as.factor(ejercicio2$Estatus)
ejercicio2$Exposicion <- as.factor(ejercicio2$Exposicion)
ejercicio2$Enfermedad <- as.factor(ejercicio2$Enfermedad)
ejercicio2$Padres.fumadores <- as.factor(ejercicio2$Padres.fumadores)
ejercicio2$IMC_escala <- as.factor(ejercicio2$IMC_escala)
str(ejercicio2) 

'data.frame':   200 obs. of  13 variables:
 $ Clave             : chr  "A001" "A002" "A003" "A004" ...
 $ Sexo              : Factor w/ 2 levels "H","M": 2 1 2 1 2 1 1 2 1 2 ...
 $ Edad              : int  21 35 38 18 32 21 39 23 37 35 ...
 $ Estatus           : Factor w/ 2 levels "Caso","Control": 2 1 2 2 1 2 2 1 1 2 ...
 $ Exposicion        : Factor w/ 4 levels "Fumador","Humo de biomasa",..: 3 2 1 1 3 1 1 4 1 2 ...
 $ Enfermedad        : Factor w/ 5 levels "Asma","Bronquitis",..: 5 3 5 5 2 5 5 2 2 5 ...
 $ Cotinina          : num  1 8 183.2 225.5 2.3 ...
 $ Anos.de.exposicion: int  0 5 20 4 0 3 15 5 20 32 ...
 $ Talla             : num  1.65 1.56 1.66 1.77 1.5 1.78 1.8 1.55 1.68 1.7 ...
 $ Peso              : int  88 55 68 104 56 115 98 74 99 65 ...
 $ Padres.fumadores  : Factor w/ 2 levels "No","Si": 1 1 2 1 2 2 2 1 1 2 ...
 $ IMC               : num  32.3 22.6 24.7 33.2 24.9 ...
 $ IMC_escala        : Factor w/ 4 levels "Desnutricion",..: 3 2 2 3 2 3 3 3 3 2 ...

Procedemos a hacer un análisis estadístico

resultados <- ejercicio2 %>% 
  get_summary_stats()    #Solo es con el paquete rstatix
resultados

Para hacer un análisis estratificado se requiere una variable cualitativa como factor

resultados <- ejercicio2 %>% 
  group_by(Enfermedad) %>% 
  get_summary_stats()    #Solo es con el paquete rstatix
resultados

Para realizar la tarea abrimos el archivo de datos de framingham, luego seleccionamos las variables de las cuales queremos obtener las medidas de tendencia central y las medidas de dispersión utilizando el paquete rstatix

datos_framingham <- read_csv("/cloud/project/Bioestadistica_R/framingham.csv")

Rows: 4238 Columns: 16── Column specification ──────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (16): male, age, education, currentSmoker, cigsPerDay, BPMeds, prevalentStroke, preval...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

View(datos_framingham)

str(datos_framingham)

spc_tbl_ [4,238 × 16] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ male           : num [1:4238] 1 0 1 0 0 0 0 0 1 1 ...
 $ age            : num [1:4238] 39 46 48 61 46 43 63 45 52 43 ...
 $ education      : num [1:4238] 4 2 1 3 3 2 1 2 1 1 ...
 $ currentSmoker  : num [1:4238] 0 0 1 1 1 0 0 1 0 1 ...
 $ cigsPerDay     : num [1:4238] 0 0 20 30 23 0 0 20 0 30 ...
 $ BPMeds         : num [1:4238] 0 0 0 0 0 0 0 0 0 0 ...
 $ prevalentStroke: num [1:4238] 0 0 0 0 0 0 0 0 0 0 ...
 $ prevalentHyp   : num [1:4238] 0 0 0 1 0 1 0 0 1 1 ...
 $ diabetes       : num [1:4238] 0 0 0 0 0 0 0 0 0 0 ...
 $ totChol        : num [1:4238] 195 250 245 225 285 228 205 313 260 225 ...
 $ sysBP          : num [1:4238] 106 121 128 150 130 ...
 $ diaBP          : num [1:4238] 70 81 80 95 84 110 71 71 89 107 ...
 $ BMI            : num [1:4238] 27 28.7 25.3 28.6 23.1 ...
 $ heartRate      : num [1:4238] 80 95 75 65 85 77 60 79 76 93 ...
 $ glucose        : num [1:4238] 77 76 70 103 85 99 85 78 79 88 ...
 $ TenYearCHD     : num [1:4238] 0 0 0 1 0 0 1 0 0 0 ...
 - attr(*, "spec")=
  .. cols(
  ..   male = col_double(),
  ..   age = col_double(),
  ..   education = col_double(),
  ..   currentSmoker = col_double(),
  ..   cigsPerDay = col_double(),
  ..   BPMeds = col_double(),
  ..   prevalentStroke = col_double(),
  ..   prevalentHyp = col_double(),
  ..   diabetes = col_double(),
  ..   totChol = col_double(),
  ..   sysBP = col_double(),
  ..   diaBP = col_double(),
  ..   BMI = col_double(),
  ..   heartRate = col_double(),
  ..   glucose = col_double(),
  ..   TenYearCHD = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 

datos_framingham <- rename(datos_framingham, sexos = male)
datos_framingham

resultados_framingham <- datos_framingham %>% 
  select(age, sexos, cigsPerDay, totChol, sysBP, diaBP, BMI, heartRate, glucose) %>% 
  filter(!complete.cases(.)) %>%
  group_by(sexos) %>% 
  get_summary_stats()

View(resultados_framingham)

Ahora obtenemos el rango de manera manual calculando la diferencia entre el máximo y el mínimo.

rango_age <- max(datos_framingham$age, na.rm = TRUE) - min(datos_framingham$age, na.rm = TRUE)
rango_cigs <- max(datos_framingham$cigsPerDay, na.rm = TRUE) - min(datos_framingham$cigsPerDay, na.rm = TRUE)
rango_chol <- max(datos_framingham$totChol, na.rm = TRUE) - min(datos_framingham$totChol, na.rm = TRUE)
rango_sbp <- max(datos_framingham$sysBP, na.rm = TRUE) - min(datos_framingham$sysBP, na.rm = TRUE)
rango_dbp <- max(datos_framingham$diaBP, na.rm = TRUE) - min(datos_framingham$diaBP, na.rm = TRUE)
rango_BMI <- max(datos_framingham$BMI, na.rm = TRUE) - min(datos_framingham$BMI, na.rm = TRUE)
rango_hr <- max(datos_framingham$heartRate, na.rm = TRUE) - min(datos_framingham$heartRate, na.rm = TRUE)
rango_gluc <- max(datos_framingham$glucose, na.rm = TRUE) - min(datos_framingham$glucose, na.rm = TRUE)

rango_age

[1] 38

rango_cigs

[1] 70

rango_chol

[1] 589

rango_sbp

[1] 211.5

3. Anális descriptivo y normalidad (Sesión 3), con enfoque en el análisis de normalidad

Tener precaución para utilizar read.csvo read_csvporque el manejo de las tablas es diferente

library(tidyverse)

tarea_1 <- read.csv("/cloud/project/Bioestadistica_R/framingham.csv")
str(tarea_1)

'data.frame':   4238 obs. of  16 variables:
 $ male           : int  1 0 1 0 0 0 0 0 1 1 ...
 $ age            : int  39 46 48 61 46 43 63 45 52 43 ...
 $ education      : int  4 2 1 3 3 2 1 2 1 1 ...
 $ currentSmoker  : int  0 0 1 1 1 0 0 1 0 1 ...
 $ cigsPerDay     : int  0 0 20 30 23 0 0 20 0 30 ...
 $ BPMeds         : int  0 0 0 0 0 0 0 0 0 0 ...
 $ prevalentStroke: int  0 0 0 0 0 0 0 0 0 0 ...
 $ prevalentHyp   : int  0 0 0 1 0 1 0 0 1 1 ...
 $ diabetes       : int  0 0 0 0 0 0 0 0 0 0 ...
 $ totChol        : int  195 250 245 225 285 228 205 313 260 225 ...
 $ sysBP          : num  106 121 128 150 130 ...
 $ diaBP          : num  70 81 80 95 84 110 71 71 89 107 ...
 $ BMI            : num  27 28.7 25.3 28.6 23.1 ...
 $ heartRate      : int  80 95 75 65 85 77 60 79 76 93 ...
 $ glucose        : int  77 76 70 103 85 99 85 78 79 88 ...
 $ TenYearCHD     : int  0 0 0 1 0 0 1 0 0 0 ...

Cambiamos el tipo de variable malecomo factor

tarea_1$male <- as.factor(tarea_1$male)
str(tarea_1)

'data.frame':   4238 obs. of  16 variables:
 $ male           : Factor w/ 2 levels "0","1": 2 1 2 1 1 1 1 1 2 2 ...
 $ age            : int  39 46 48 61 46 43 63 45 52 43 ...
 $ education      : int  4 2 1 3 3 2 1 2 1 1 ...
 $ currentSmoker  : int  0 0 1 1 1 0 0 1 0 1 ...
 $ cigsPerDay     : int  0 0 20 30 23 0 0 20 0 30 ...
 $ BPMeds         : int  0 0 0 0 0 0 0 0 0 0 ...
 $ prevalentStroke: int  0 0 0 0 0 0 0 0 0 0 ...
 $ prevalentHyp   : int  0 0 0 1 0 1 0 0 1 1 ...
 $ diabetes       : int  0 0 0 0 0 0 0 0 0 0 ...
 $ totChol        : int  195 250 245 225 285 228 205 313 260 225 ...
 $ sysBP          : num  106 121 128 150 130 ...
 $ diaBP          : num  70 81 80 95 84 110 71 71 89 107 ...
 $ BMI            : num  27 28.7 25.3 28.6 23.1 ...
 $ heartRate      : int  80 95 75 65 85 77 60 79 76 93 ...
 $ glucose        : int  77 76 70 103 85 99 85 78 79 88 ...
 $ TenYearCHD     : int  0 0 0 1 0 0 1 0 0 0 ...

tarea_1 <- tarea_1 %>% 
  mutate(male = case_when(male == "0" ~ "M",
                          male == "1" ~ "H"))

str(tarea_1)

'data.frame':   4238 obs. of  16 variables:
 $ male           : chr  "H" "M" "H" "M" ...
 $ age            : int  39 46 48 61 46 43 63 45 52 43 ...
 $ education      : int  4 2 1 3 3 2 1 2 1 1 ...
 $ currentSmoker  : int  0 0 1 1 1 0 0 1 0 1 ...
 $ cigsPerDay     : int  0 0 20 30 23 0 0 20 0 30 ...
 $ BPMeds         : int  0 0 0 0 0 0 0 0 0 0 ...
 $ prevalentStroke: int  0 0 0 0 0 0 0 0 0 0 ...
 $ prevalentHyp   : int  0 0 0 1 0 1 0 0 1 1 ...
 $ diabetes       : int  0 0 0 0 0 0 0 0 0 0 ...
 $ totChol        : int  195 250 245 225 285 228 205 313 260 225 ...
 $ sysBP          : num  106 121 128 150 130 ...
 $ diaBP          : num  70 81 80 95 84 110 71 71 89 107 ...
 $ BMI            : num  27 28.7 25.3 28.6 23.1 ...
 $ heartRate      : int  80 95 75 65 85 77 60 79 76 93 ...
 $ glucose        : int  77 76 70 103 85 99 85 78 79 88 ...
 $ TenYearCHD     : int  0 0 0 1 0 0 1 0 0 0 ...

View(tarea_1)

Realizamos el análisis descriptivo

resultados_tarea_1 <- tarea_1 %>% 
  group_by(male) %>% 
  get_summary_stats(age, cigsPerDay, totChol, sysBP, diaBP, BMI, heartRate, glucose)

resultados_tarea_1 <- resultados_tarea_1 %>% 
  mutate(Rango_1 = max - min)
  
View(resultados_tarea_1)

Para crear un archivo .csv

write.csv(resultados_tarea_1, file = "resultados_tarea1")

Temas de estadística

Definir población y muestra. La muestra es el conjunto de elementos representativos de la población y se debe seleccionar aleatoriamente.

Escalas de medición: 4 niveles de escala

• Nominal: Se utilzan nombres, símbolos o numerales, que son empleados como etiquetas. Sirve para separar por grupos. No se pueden aplicar operaciones matemáticas.

• Ordinal: Los números o etiquetas utilizadas indican un orden. Sirve para organizar.

• Intervalos: Establece un orden determinado. El 0 significa ausencia de valor. Se puede utilizar para operaciones matemáticas.

• Proporción/cociente/razón: Son secuencias numéricas. Se pueden realizar operaciones matemáticas. Es la escala que permite comparar de mejor manera los datos.

Distribución de datos

• Distribución normal
• Distribución de Poisson
• Distribución binomial

Asimetría de datos

• Negativa
• Simétrica
• Positiva

Curtosis

• Leptocúrtica
• Mesocúrtica
• Platicúrtica

4. La normalidad…

Puede ser evaluada mediante 3 modelos

• Modelos gráficos: histogramas y qqplot
• Media = Mediana = Moda
• Pruebas de normalidad:
  • Shapiro-Wild: para grupos con menos de 50 observaciones
  • Kolmogorov-Smirnov/Lillifors: para grupos con más de 50 observaciones

Representación de los datos a partir de las medidas de tendencia central

Utilizamos el valor de la edad

ejercicio_1 <- tarea_1 %>% 
  get_summary_stats(age)

conteo <- tarea_1 %>% 
  count(age)

#mean: 49.58
#median: 49
#moda: 40

Solo con estos parámetros se podría decir que la muestra no cumple con la normalidad. Veamos su historgrama con la función hist(), del cual podemos modificar su título, las etiquetas y el color.

hist(tarea_1$age,
     main = "Histograma de edad",
     ylab = "Frecuencia",
     xlab = "Edad (años)",
     col = "salmon",
     breaks = 50)

Podemos realizar las gráficas con el paquete ggplot que permite crear gráficos más estéticos

library(ggplot)

Error in library(ggplot) : there is no package called ‘ggplot’

Para graficar el histograma

Construyamos los gráficos de densidad

Revisar la asimetría y la curtosis

library(car)

Loading required package: carData

Attaching package: ‘car’

The following object is masked from ‘package:psych’:

    logit

La curtosis ideal debería ser cercano a 0 para decir que es mesocúrtica. Es decir,

• k = 0 -> mesocúrtica
• k > 0 -> leptocúrtica
• k < 0 -> platicúrtica

Para calcular la asimetría, decimos que

• a = 0 (simétrico)
• a > 0 (asimetría positiva)
• a < 0 (asimetría negativa)

Para calcular la curtosis utilizamos el paquete psych

kurtosi(tarea_1$age)

[1] -0.991

Cuando el valor de curtosis es < 2 o > 2, entonces no tiene un comportamiento normal

Calculamos la curtosis y el sesgo para hombres y para mujeres

kurtosi(tarea_1$age[tarea_1$male == "H"])

[1] -0.963

kurtosi(tarea_1$age[tarea_1$male == "M"])

[1] -1.01

#Para calcular el sesgo (simetría o asimetría)
skew(tarea_1$age)

[1] 0.228

skew(tarea_1$age[tarea_1$male == "H"])

[1] 0.272

skew(tarea_1$age[tarea_1$male == "M"])

[1] 0.195

Veamos los QQ-plot

qqPlot(tarea_1$age,
       main = "QQ plot edad",
       ylab = "Edad (años)",
       xlab = "Cuantiles")

[1] 1625 3138

Combinar análisis de asimetría con histograma, y curtosis con QQ-plot. Los QQplot permiten observar el peso de las colas. Si es una distribución normal, entonces se debería de comportar como una línea recta.

Podemos realizar el QQplot estratificado por hombres y mujeres.

qqPlot(tarea_1$age[tarea_1$male == "H"],
       main = "QQ plot edad",
       ylab = "Edad (años)",
       xlab = "Cuantiles",
       col = "blue",
       col.lines = "red")

Gráficos de cajas y bigotes son importantes para el análisis

Utilizando ggplot

Pruebas de normalidad

• Shapiro-Wilks: Se utiliza cuando n < 50 (incluso menor a 40). Cuando hay menos de 10 datos ya no se recomienda usar.
• Kolmogorov-Smirnov: Se utiliza cuando n > 50

Para realizar la de Kolmogorov-Smirnov…Con los requisitos mínimo Tiene como \(H_0\): La distribución real y la distribución teórica son iguales o similares (normalidad) \(H_a\): La distribución real y la distribución teórica son diferentes (no normalidad)

Cuando p > 0.05 no se rechaza la \(H_0\)

OJO: ya no se usa así, se usa una mejorada

ks.test(tarea_1$age, "pnorm",
        mean = 49.6, sd = 8.57,
        alternative = "two.sided")

Warning: ties should not be present for the Kolmogorov-Smirnov test

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  tarea_1$age
D = 0.08, p-value <2e-16
alternative hypothesis: two-sided

Como el el valor de p obtenido es < 0.05 (p-value < 2e-16), se rechaza la hipótesis nula y concluimos que la muestra no proviene de una distribución normal.

La prueba de Kolmogorov-Smirnov-Lilliefors es la corregida

# Se utiliza la librería de nortest
lillie.test(tarea_1$age[tarea_1$male == "H"])

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  tarea_1$age[tarea_1$male == "H"]
D = 0.09, p-value <2e-16

Para hacer la de Shapiro-Wilk

shapiro.test(tarea_1$age[tarea_1$male == "H"])

    Shapiro-Wilk normality test

data:  tarea_1$age[tarea_1$male == "H"]
W = 1, p-value <2e-16

shapiro.test(tarea_1$age[tarea_1$male == "M"])

    Shapiro-Wilk normality test

data:  tarea_1$age[tarea_1$male == "M"]
W = 1, p-value <2e-16

4. Pruebas de contraste estadístico

Pruebas paramétricas (t de student y ANOVA)

Requieren cumplir varios supuestos

• Normalidad
• Igualdad de varianzas (Prueba de Levene): Requiere el planteamiento de las siguientes hipótesis
  • \(H_0\): la varianza de la variable es igual o semejante en casos y controles (\(p>0.05\))
  • \(H_a\): las varianzas son diferentes entre grupos (\(p<0.05\))

Las hipotesis planteadas en las pruebas de comparación de medias son:

• \(H_0: \mu_1 = \mu_2\)
• \(H_a: \mu_1 \neq \mu_2\)

Pruebas no paramétricas

• Signos de Wilcoxon
• U de Mann-Whitney
• Friedman
• Kruskal-Wallis

Empecemos con las pruebas paramétricas

Cargamos nuestra base de datos

sesion_4 <- read.csv("/cloud/project/Bioestadistica_R/datos_1.csv")
View(sesion_4)

Pasos a seguir para realizar la prueba:

1. Realizar el análisis de normalidad
2. Análisis de igualdad varianza (Prueba de Levene)
3. Prueba de t de student
4. Cálculo de tamaño de efecto
5. Conclusión

Comparemos la talla entre casos y controles. ¿Existe diferencia en esta variable entre ambos grupos?

1. Análisis de normalidad

summary(sesion_4$Estatus) #Debe de estar en tipo factor

   Caso Control 
    100     100 

Dado que tenemos > 50 observaciones para cada grupo, utilizamos la prueba de Lilliefors. El análisis debe de ser realizado para cada grupo por separado, recordando que

\(H_0\): La distribución real y teórica son iguales ($p > 0.05) \(H_a\): La distribución real y teórica son diferentes ($p < 0.05)

lillie.test(sesion_4$Talla[sesion_4$Estatus == "Control"])

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  sesion_4$Talla[sesion_4$Estatus == "Control"]
D = 0.1, p-value = 0.02

lillie.test(sesion_4$Talla[sesion_4$Estatus == "Caso"])

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  sesion_4$Talla[sesion_4$Estatus == "Caso"]
D = 0.1, p-value = 0.006

Dado que \(p-value < 0.05\), se recomienda análisis no paramétrico para comparar grupos.

2. Ahora verificamos el análisis de igualdad de varianza (Prueba de Levene) \(H_0\): Las varianzas entre grupos son iguales o similares ($p > 0.05) \(H_a\): Las varianzas entre grupos son diferentes ($p < 0.05)

leveneTest(sesion_4$Talla ~ sesion_4$Estatus,
            center = "median")

Levene's Test for Homogeneity of Variance (center = "median")
       Df F value Pr(>F)
group   1       0   0.96
      198               

leveneTest(sesion_4$Talla ~ sesion_4$Estatus,
            center = "mean")

Levene's Test for Homogeneity of Variance (center = "mean")
       Df F value Pr(>F)
group   1       0   0.99
      198               

El p-value (Pr(>F)) > 0.05, entonces se acepta la \(H_0\). Las varianzas entre grupos son iguales o similares.

3. Ahora realizamos la prueba de T

Las hipótesis son

\(H_0: \mu_1 = \mu_2\) \(H_a: \mu_1 \neq \mu_2\)

#Si var.eq = TRUE, entonces hace la t de student normal
t.test(sesion_4$Talla ~ sesion_4$Estatus,
       alternative = "two.sided",
       conf.level = 0.95,
       var.equal = TRUE,
       paired = FALSE)

    Two Sample t-test

data:  sesion_4$Talla by sesion_4$Estatus
t = -0.6, df = 198, p-value = 0.5
alternative hypothesis: true difference in means between group Caso and group Control is not equal to 0
95 percent confidence interval:
 -0.0324  0.0164
sample estimates:
   mean in group Caso mean in group Control 
                 1.67                  1.67 

#Si var.eq = FALSE, entonces hace la corrección de Welch porque las varianzas entre grupos son diferentes
#Esta prueba es para varianzas NO iguales
t.test(sesion_4$Talla ~ sesion_4$Estatus,
       alternative = "two.sided",
       conf.level = 0.95,
       var.equal = FALSE,
       paired = FALSE)

    Welch Two Sample t-test

data:  sesion_4$Talla by sesion_4$Estatus
t = -0.6, df = 198, p-value = 0.5
alternative hypothesis: true difference in means between group Caso and group Control is not equal to 0
95 percent confidence interval:
 -0.0324  0.0164
sample estimates:
   mean in group Caso mean in group Control 
                 1.67                  1.67 

#Pacticamente esta es igual a la anterior
t.test(sesion_4$Talla ~ sesion_4$Estatus)

    Welch Two Sample t-test

data:  sesion_4$Talla by sesion_4$Estatus
t = -0.6, df = 198, p-value = 0.5
alternative hypothesis: true difference in means between group Caso and group Control is not equal to 0
95 percent confidence interval:
 -0.0324  0.0164
sample estimates:
   mean in group Caso mean in group Control 
                 1.67                  1.67 

El intervalo de confianza (IC) son -0.0324 y 0.0164

Dado que el p-valor (0.5) > 0.05, entonces no se rechaza la \(H_0\) y las medias entre ambos grupos son iguales. El cálculo del tamaño del efecto NO se aplica cuando p > 0.05, pero vamos a verificarlo

cohen.d(sesion_4$Talla, sesion_4$Estatus)

Call: cohen.d(x = sesion_4$Talla, group = sesion_4$Estatus)
Cohen d statistic of difference between two means
     lower effect upper
[1,] -0.19   0.09  0.37

Multivariate (Mahalanobis) distance between groups
[1] 0.092
r equivalent of difference between two means
data 
0.05 

Otra opción de anális usando rstatix

library(rstatix)
library(tidyverse)
var_igualdad <- sesion_4 %>% 
  levene_test(Talla ~ Estatus, center = "median")

efecto <- sesion_4 %>% 
  cohens_d(Talla ~ Estatus)

5. ANOVA (Sesión 5)

Cargamos la base de datos

Preparación de la base de datos

str(pulmonar)

'data.frame':   130 obs. of  10 variables:
 $ Clave      : chr  "P001" "P002" "P003" "P004" ...
 $ Estatus    : Factor w/ 2 levels "Case","Control": 2 2 1 2 1 2 2 2 1 1 ...
 $ Sexo       : Factor w/ 2 levels "H","M": 1 2 1 2 2 2 2 1 2 1 ...
 $ Dosis      : Factor w/ 3 levels "Clásico","Nuevo",..: 3 1 1 2 3 3 2 3 1 3 ...
 $ Edad       : int  62 55 79 69 70 63 58 50 47 88 ...
 $ Peso       : num  65.9 68.5 72.1 92.9 84.9 ...
 $ Talla      : num  2.28 2.1 1.52 2.1 1.86 1.55 1.73 2.71 1.79 1.93 ...
 $ Fun_pulm   : num  86.9 89.3 63.6 86.3 107.8 ...
 $ PvCO2      : num  13.5 20.2 28.8 19.8 42.1 ...
 $ Neutrofilos: num  8.66 7.35 23.05 7.44 21.79 ...

head(pulmonar)

Empecemos con el proceso de análisis.

1. Primero contemos el número de datos, n.

Dado que para cada uno de los grupos su n es < 50, utilizamos la prueba de Shapiro-Wilks

p <- pulmonar %>% 
  ggplot(pulmonar)+
  geom_bar(aes(x = Dosis, y = PvCO2))

Error in `ggplot()`:
! `mapping` should be created with `aes()`.
✖ You've supplied a <data.frame> object
Backtrace:
 1. pulmonar %>% ggplot(pulmonar)
 3. ggplot2:::ggplot.default(., pulmonar)

Las hipótesis son

• \(H_0\): La distribución real \(=\) la distribución teórica (normalidad, p-value > 0.05)
• \(H_a\): La distribución real \(\neq\) la distribución teórica son (no normalidad, p-value < 0.05)

Dado que el análisis de normalidad indicó que el valor p de la variable PvCO2 para cada uno de los grupos fue menor que 0.05, entonces se rechaza la hipótesis nula y aceptamos la \(H_a\), que es que no provienen de una distribución normal. Ahorita omitimos que sean no normales. Para realizar la ANOVA, todos los grupos deben cumplir el supuesto de normalidad

EL siguiente paso es comprobar la igualdad de varianzas

2. Igualdad de varianzas. Para este caso se usa la paquetería car y carData

leveneTest(pulmonar$PvCO2 ~ pulmonar$Dosis, 
           center = "median")

Levene's Test for Homogeneity of Variance (center = "median")
       Df F value Pr(>F)
group   2    0.41   0.67
      127               

Las hipótesis son - \(H_0\): Las varianzas de PvCO2 \(=\) en todos los grupos (p-value > 0.05) - \(H_a\): Las varianzas de PvCO2 \(\neq\) en al menos un grupo (no normalidad, p-value < 0.05)

Usando el paquete rstatix

pulmonar %>% 
  levene_test(PvCO2 ~ Dosis, center = "median")

Dado que el p-value es > 0.05 (p = 0.667), entonces la conclusión del analisis de igualdad de varianzas es que las varianzas son similares o iguales. Procedemos a realizas la prueba de ANOVA

3. ANOVA Las hipótesis son

• \(H_0: \mu_1 = \mu_2 = \mu_3\)
• \(H_a: \mu_1 \neq \mu_2 \neq \mu_3\):

La \(H_a\) dice que hay diferencias en la media de PvCO2 entre alguno de los grupos

anova

ANOVA Table (type II tests)

  Effect DFn DFd   F    p p<.05   ges
1  Dosis   2 127 1.4 0.25       0.022

Dado que p = 0.25, que es > 0.05, entonces aceptamos la \(H_0\) y diríamos que no hay diferencia entre grupos. Haremos las pruebas post-hoc

Finalmente, realizamos la conclusión

5. Los tratamientos no provocan una diferencia en la PvCO2 (p-value = 0.25)

ANOVA de 2 vías

El análisis consiste

1. Normalidad

pulmonar %>% 
  group_by(Dosis) %>% 
  shapiro_test(Neutrofilos)


summary(pulmonar$Estatus)

   Case Control 
     40      90 

#Para los casos n=40, usamos la de Shapiro-Wilks
shapiro_test(pulmonar$Neutrofilos[pulmonar$Estatus == "Case"])

#Para controles n=90 usamos la de lillieforce
library(nortest)
lillie.test(pulmonar$Neutrofilos[pulmonar$Estatus == "Control"])

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  pulmonar$Neutrofilos[pulmonar$Estatus == "Control"]
D = 0.06, p-value = 0.6

Para el % de neutrofilos estratificado por estatus cumplen la normalidad. SIn embargo, al estratificar por el tratamiento no cumple.

normalidad_condicionada <- pulmonar %>% 
  group_by(Dosis, Estatus) %>% 
  summary() %>% 
  shapiro_test(Neutrofilos)

Error in UseMethod("select") : 
  no applicable method for 'select' applied to an object of class "table"

Cuando se condiciona el análisis de normalidad de dos factores, el % de neutrófilos cumple con el supuesto de normalidad

2. Igualdad de varianzas

La varianza de % de neutrófilos es diferente entre los grupos (Estatus). La varianza de neutrófilos es igual entre grupos (Dosis)

3. ANOVA

result_2vias <- pulmonar %>% 
  anova_test(Neutrofilos ~ Dosis*Estatus, 
             effect.size = "ges") #Eta cuadrada generalizada

Con el análisis post-hoc

ppst_2vias <-pulmonar %>% 
  tukey_hsd(Neutrofilos ~ Dosis*Estatus)

4. Conclusión Analizando los subgrupos dosis:casos se encontró que al comparar los múltiples tratamientos entre casos y controles, no hay diferencias significativas. Al realizar un análisis individual de cada factor, se detectó una gran diferencia en el factor Estatus (p < 0.01, ges = 0.66), por lo que se sospecha que las diferencias que puedan encontrarse son atribuibles al factor Estatus.

ANOVA de 2 vías de medidas repetidas

Cargamos la base de datos y convertimos las variables caracter a factor.

str(seguimiento)

'data.frame':   84 obs. of  10 variables:
 $ Clave              : chr  "C-001" "C-002" "C-004" "C-005" ...
 $ Tiempo             : Factor w/ 3 levels "Antes","Despues",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Sexo               : Factor w/ 2 levels "H","M": 2 1 2 1 2 2 1 2 1 1 ...
 $ Fumadores          : Factor w/ 2 levels "No","Si": 2 2 2 2 2 2 2 2 2 2 ...
 $ Edad.de.inicio     : int  24 21 23 19 18 15 14 21 18 23 ...
 $ Cigarros.al.día    : int  43 11 45 39 38 30 44 27 49 28 ...
 $ Padres.fumadores   : Factor w/ 2 levels "No","Si": 2 2 1 2 2 1 2 2 2 2 ...
 $ Problemas.cardiacos: Factor w/ 2 levels "No","Si": 2 1 2 1 1 2 2 1 2 2 ...
 $ CHRNA5             : Factor w/ 2 levels "No","Si": 2 2 1 2 1 2 2 1 2 1 ...
 $ rs578776           : Factor w/ 3 levels "CC","CT","TT": 2 2 3 3 1 2 2 1 3 1 ...

Hay que asegurarnos que están todos los datos completos de la base de datos. Es decir, que la n sea la misma para todos los tiempos .

Analicemos la normalidad

Los valores de Mantenimiento no provienen de una distribución normal. En este caso, no podríamos continuar con la ANOVA y tendríamos que hacer otra prueba no paramétrica. Sin embargo, vamos a continuar con la ANOVA.

Ahora realicemos la prueba de igualdad de varianzas. Prueba de Levene

Dado el valor de p > 0.05 (p-value = 0.143), podemos decir que las varianzas son iguales/similares en todos los grupos.

Planteemos la prueba de ANOVA

seguimiento %>% 
  anova_test(dv = Cigarros.al.día, 
             wid = Clave,
             within = Tiempo)

ANOVA Table (type III tests)

$ANOVA
  Effect DFn DFd    F        p p<.05   ges
1 Tiempo   2  54 40.2 2.06e-11     * 0.362

$`Mauchly's Test for Sphericity`
  Effect     W     p p<.05
1 Tiempo 0.901 0.256      

$`Sphericity Corrections`
  Effect  GGe      DF[GG]   p[GG] p[GG]<.05   HFe      DF[HF]    p[HF] p[HF]<.05
1 Tiempo 0.91 1.82, 49.12 1.4e-10         * 0.972 1.94, 52.47 3.75e-11         *

Ya que el valor de p de ANOVA es < 0.05 (p-value = 2.06e-11), podemos decir que hay diferencias entre grupos.

Para el test de esfericidad

• \(H_0\): los datos no varian en cada nivel de tiempo (p >0.05)
• \(H_a\): los datos varían en cada nivel de tiempo (p < 0.05)

Si no se cumple el supuesto de esfericidad por el test de Mauchly, habrá que ver los datos de la corrección de esfericidad. Habrá que considerar el valor de p de p[GG] del test de corrección de esfericidad como el que se tiene que reportar de la prueba de ANOVA.

Hasta aquí concluimos parcialmente que existe diferencia en el número de cigarros por día después de cada etapa de tratamiento.

##6. Pruebas no paramétricas

1. Son pruebas para comparar datos con distribuciones no normales
2. No son sensibles a datos extremos atípicos, por lo que su presencia no afectará el análisis en conjunto
3. Fundament: ordena los datos de menor a mayor, cada dato se le asigna posición y un valor numérico (rango de 1 a n)

• Comparaciones cuantitativas: prueba de Wilcoxon, U de Mann-Whitney, Kruskal-Wallis, Friedman
• Comparaciones cualitativas: chi cuadrada, Fisher

str(sesion_6)

'data.frame':   200 obs. of  11 variables:
 $ Clave             : chr  "A001" "A002" "A003" "A004" ...
 $ Sexo              : chr  "M" "H" "M" "H" ...
 $ Edad              : int  21 35 38 18 32 21 39 23 37 35 ...
 $ Estatus           : chr  "Control" "Caso" "Control" "Control" ...
 $ Exposicion        : Factor w/ 4 levels "Fumador","Humo de biomasa",..: 3 2 1 1 3 1 1 4 1 2 ...
 $ Enfermedad        : chr  "Sin alteraciones" "Capacidad pulmonar alterada" "Sin alteraciones" "Sin alteraciones" ...
 $ Cotinina          : num  1 8 183.2 225.5 2.3 ...
 $ Anos.de.exposicion: int  0 5 20 4 0 3 15 5 20 32 ...
 $ Talla             : num  1.65 1.56 1.66 1.77 1.5 1.78 1.8 1.55 1.68 1.7 ...
 $ Peso              : int  88 55 68 104 56 115 98 74 99 65 ...
 $ Padres.fumadores  : chr  "No" "No" "Si" "No" ...

Hay que probar la asimetría (skew)

skew(sesion_6$Cotinina[sesion_6$Exposicion== "Sin exposicion"])

[1] 0.908

skew(sesion_6$Cotinina[sesion_6$Exposicion== "Fumador"])

[1] -0.115

skew(sesion_6$Cotinina[sesion_6$Exposicion== "Vapeador"])

[1] 1.99

skew(sesion_6$Cotinina[sesion_6$Exposicion== "Humo de biomasa"])

[1] 0.734

Asimetría < -2 o > 2 ya era una simetría muy marcada. Se espera tener una asimetría = 0 (simétricos), o cuando \(a \neq 0\) (asimétricos )

Realizamos el análisis de igualdad de varianzas

Realizamos la prueba de U de Mann-Whitney para datos no relacionados

Existe diferencia entre el grupo sin exposición vs fumadores (se report el p.adj.signif), y sin exposicion vs vapeadores.

La prueba post-hoc de Kruskal-Wallis es la siguiente

Para responder a la pregunta de qué tan grande es la diferencia, entonces calculamos el tamaño del efecto

Los niveles de cotinina varían entre el grupo sin exposición vs fumadores (p = 3.19e-16) y sin exposición vs vapeadores. Las diferencias encontradas son grandes según la prueba del tamaño del efecto.

Revisemos como se comparan entre todos los grupos en los niveles de Cotinina. Realizamos la prueba de Kruskal-Wallis

No sabemos entre que grupos hay diferencia, pero si existe porque el valor de p arrojado por la prueba de Kruskal-Wallis (p = 2.27e-32). Para saber entre que grupos hacemos la prueba post-hoc de dunn

O hacemos una prueba post-hoc

Para ver el tamaño del efecto

Pruebas cualitativas

Realicemos una tabla de contingencia (principalmente para variables cualitativas)

estatus_sexo

    Estatus
Sexo Caso Control
   H   58      58
   M   42      42

Una tabla de nx2 (Primero filas, luego columnas)

table(sesion_6$Sexo, sesion_6$Exposicion,
      dnn = c("Sexo", "Exposicion"))

    Exposicion
Sexo Fumador Humo de biomasa Sin exposicion Vapeador
   H      55              17             18       26
   M      35              20             14       15

Para realizar la prueba

chisq.test(estatus_exposicion)

    Pearson's Chi-squared test

data:  estatus_exposicion
X-squared = 3, df = 3, p-value = 0.4

cuando p > 0.05 la freucuencia de datos es igual, cuando p < 0.05 la frecuenca de datos es diferente