#install.packages("conflicted")
#install.packages("dplyr")
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(gtsummary)
library(gt)
library(readxl)
library(conflicted)
library(dplyr)
setwd(dir = "/Users/lorandacalderonzamora/Downloads/")
data <- read_excel("JuanLuis_DB.xlsx")
class(data)
## [1] "tbl_df"     "tbl"        "data.frame"
names(data)
## [1] "px"         "grupo"      "sexo"       "tiempo"     "OD"        
## [6] "caradental" "ICDAS"      "edad"
str(data)
## tibble [3,600 × 8] (S3: tbl_df/tbl/data.frame)
##  $ px        : num [1:3600] 1 1 1 1 1 1 1 1 1 1 ...
##  $ grupo     : num [1:3600] 1 1 1 1 1 1 1 1 1 1 ...
##  $ sexo      : num [1:3600] 0 0 0 0 0 0 0 0 0 0 ...
##  $ tiempo    : num [1:3600] 0 0 0 0 0 0 0 0 0 0 ...
##  $ OD        : num [1:3600] 0 0 0 1 1 1 2 2 2 3 ...
##  $ caradental: num [1:3600] 0 1 2 0 1 2 0 1 2 0 ...
##  $ ICDAS     : num [1:3600] 0 0 0 0 0 0 0 0 2 2 ...
##  $ edad      : num [1:3600] 2 2 2 2 2 2 2 2 2 2 ...
data <- as.data.frame(data)

# Recodificar variables, asegurándonos de incluir todas las categorías
data <- data %>%
  mutate(
    Aplicación_de_Fluor = dplyr::recode(grupo, `0` = "Control", `1` = "Tratamiento"),
    Sexo = dplyr::recode(sexo, `0` = "Femenino", `1` = "Masculino"),
    Duracion_del_tratamiento = dplyr::recode(tiempo, 
                                             `0` = "Antes de los frenos", 
                                             `1` = "2M después de los frenos", 
                                             `2` = "4M después de los frenos"),
    Cara_dental = dplyr::recode(caradental, 
                                `0` = "oclusal", 
                                `1` = "palatina", 
                                `2` = "vestibular"),
    Edad = dplyr::recode(edad, 
                         `0` = "15-25 años", 
                         `1` = "26-35 años", 
                         `2` = "36-50 años"),
    ICDAS = as.numeric(ICDAS)  # Asegurarse de que ICDAS sea numérico para análisis
  )

# Resumir los datos por paciente, tomando el promedio de las variables ICDAS y edad
summary_data <- data %>%
  group_by(px) %>%
  summarise(
    Aplicación_de_Fluor = first(Aplicación_de_Fluor),  # Mantener la categoría de tratamiento/control
    Sexo = first(Sexo),                                # Mantener el sexo
    Edad = first(Edad),                                # Mantener la edad
    ICDAS = mean(ICDAS, na.rm = TRUE)                  # Calcular el promedio del ICDAS
  )

# Crear la tabla de resumen por paciente
summary_table <- summary_data %>%
  tbl_summary(
    by = Aplicación_de_Fluor,  # Comparar por grupo de tratamiento
    missing = "no",
    type = list(
      Edad ~ "categorical",
      Sexo ~ "categorical",
      ICDAS ~ "continuous"  # ICDAS es continua
    )
  ) %>%
  add_overall() %>%
  add_p() %>%
  modify_header(label = "**Variable**") %>%
  modify_caption("**Resumen de pacientes según la aplicación de fluor**") %>%
  as_gt()
## The following warnings were returned during `as_gt()`:
## ! For variable `ICDAS` (`Aplicación_de_Fluor`) and "estimate", "statistic",
##   "p.value", "conf.low", and "conf.high" statistics: cannot compute exact
##   p-value with ties
## ! For variable `ICDAS` (`Aplicación_de_Fluor`) and "estimate", "statistic",
##   "p.value", "conf.low", and "conf.high" statistics: cannot compute exact
##   confidence intervals with ties
# Mostrar la tabla
summary_table
Resumen de pacientes según la aplicación de fluor
Variable Overall
N = 201		 Control
N = 101		 Tratamiento
N = 101		 p-value2
px 10.5 (5.5, 15.5) 15.5 (13.0, 18.0) 5.5 (3.0, 8.0) <0.001
Sexo


>0.9
    Femenino 10 (50%) 5 (50%) 5 (50%)
    Masculino 10 (50%) 5 (50%) 5 (50%)
Edad


0.3
    15-25 años 7 (35%) 5 (50%) 2 (20%)
    26-35 años 1 (5.0%) 0 (0%) 1 (10%)
    36-50 años 12 (60%) 5 (50%) 7 (70%)
ICDAS 0.83 (0.64, 0.98) 0.87 (0.64, 0.92) 0.76 (0.57, 1.04) >0.9
1 Median (Q1, Q3); n (%)
2 Wilcoxon rank sum exact test; Pearson’s Chi-squared test; Fisher’s exact test; Wilcoxon rank sum test 
library(ggplot2)
library(ggpubr)

# Resumir los datos por paciente, tomando el promedio de ICDAS por individuo
summary_data <- data %>%
  group_by(px) %>%
  summarise(
    Aplicación_de_Fluor = first(Aplicación_de_Fluor),  # Tratamiento/control
    ICDAS = mean(ICDAS, na.rm = TRUE)                 # Promedio del ICDAS por paciente
  )
ggplot(summary_data, aes(x = Aplicación_de_Fluor, y = ICDAS, fill = Aplicación_de_Fluor)) +
  geom_boxplot() +
  labs(title = "Comparación de ICDAS entre Tratamiento y Control",
       x = "Grupo de Tratamiento",
       y = "Promedio de ICDAS") +
  theme_minimal() +
  scale_fill_manual(values = c("Control" = "blue", "Tratamiento" = "orange")) + 
  stat_compare_means(method = "t.test", label = "p.format", label.x = 1.5)  

data$tiempo <- factor(data$tiempo, levels = c(0, 1, 2), labels = c("Antes", "2 meses", "4 meses"))

# Crear el boxplot y ajustar la escala del eje y
ggplot(data, aes(x = tiempo, y = ICDAS, fill = tiempo)) +
  geom_boxplot() +
  labs(title = "Comparación de ICDAS en diferentes tiempos",
       x = "Tiempo",
       y = "ICDAS") +
  theme_minimal() +
  scale_fill_manual(values = c("Antes" = "blue", "2 meses" = "orange", "4 meses" = "darkgreen")) +
  scale_y_continuous(limits = c(0, 6), breaks = seq(0, 6, by = 1)) +  # Ajustar la escala del eje y de 0 a 6
  stat_compare_means(method = "anova", label = "p.format")  # Mostrar valor p con ANOVA
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_compare_means()`).

data$Edad <- factor(data$edad, levels = c(0, 1, 2), labels = c("15-25", "26-35", "36-50"))
# Realizar ANOVA para comparar los valores de ICDAS entre los grupos de edad
anova_model <- aov(ICDAS ~ Edad, data = data)
summary(anova_model)
##               Df Sum Sq Mean Sq F value  Pr(>F)   
## Edad           2     11   5.356   4.867 0.00775 **
## Residuals   3597   3958   1.100                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Visualizar los resultados del ANOVA
ggplot(data, aes(x = Edad, y = ICDAS, fill = Edad)) +
  geom_boxplot() +
  labs(title = "Comparación de ICDAS entre los Grupos de Edad",
       x = "Grupo de Edad",
       y = "ICDAS") +
  theme_minimal() +
  scale_fill_manual(values = c("15-25" = "blue", "26-35" = "orange", "36-50" = "darkgreen")) +
  scale_y_continuous(limits = c(0, 6), breaks = seq(0, 6, by = 1)) +
  stat_compare_means(method = "anova", label = "p.format")  # Mostrar valor p en el gráfico
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_compare_means()`).

# Tomar solo la primera observación por cada individuo y cara dental para el análisis
data_ajustada <- data %>%
  group_by(px, Cara_dental) %>%
  dplyr::filter(ICDAS > 0) %>%
  summarise(ICDAS = first(ICDAS))  # Tomar la primera observación del ICDAS por individuo y cara dental
## `summarise()` has grouped output by 'px'. You can override using the `.groups`
## argument.
# Realizar ANOVA para comparar ICDAS entre las caras dentales
anova_cara_dental <- aov(ICDAS ~ Cara_dental, data = data_ajustada)
summary(anova_cara_dental)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Cara_dental  2  20.83  10.417   40.95 9.46e-12 ***
## Residuals   57  14.50   0.254                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(data_ajustada, aes(x = Cara_dental, y = ICDAS, fill = Cara_dental)) +
  geom_boxplot() +
  labs(title = "Comparación de ICDAS entre las caras dentales (por individuo)",
       x = "Cara dental",
       y = "ICDAS") +
  theme_minimal() +
  scale_fill_manual(values = c("oclusal" = "blue", "palatina" = "orange", "vestibular" = "darkgreen")) +  # Escala correcta para caras dentales
  scale_y_continuous(limits = c(0, 6), breaks = seq(0, 6, by = 1)) +  
  stat_compare_means(method = "anova", label = "p.format")  

t_test_genero <- t.test(ICDAS ~ Sexo, data = data)
print(t_test_genero)
## 
##  Welch Two Sample t-test
## 
## data:  ICDAS by Sexo
## t = -1.7462, df = 3587.3, p-value = 0.08085
## alternative hypothesis: true difference in means between group Femenino and group Masculino is not equal to 0
## 95 percent confidence interval:
##  -0.129724579  0.007502357
## sample estimates:
##  mean in group Femenino mean in group Masculino 
##               0.7716667               0.8327778
ggplot(data, aes(x = Sexo, y = ICDAS, fill = Sexo)) +
  geom_boxplot() +
  labs(title = "Comparación de ICDAS entre hombres y mujeres",
       x = "Género",
       y = "ICDAS") +
  theme_minimal() +
  scale_fill_manual(values = c("Femenino" = "pink", "Masculino" = "blue")) +
  scale_y_continuous(limits = c(0, 6), breaks = seq(0, 6, by = 1)) +
  stat_compare_means(method = "t.test", label = "p.format") 
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_compare_means()`).