Análisis descriptivo

library(readxl)
library(dplyr)
BASE <- read_excel("BD_articuloaceptabilidad_autotomaVPH.xlsx") #cargue base datos
BASE1 <- BASE %>% select(edad, tipo_afiliacion, vacuna_vph, meses_lastscreen, n_partos, resultado_vph) #seleccionar 
BASE1 <- as.data.frame(BASE1) #manejo como dataframe
#str(BASE1)


BASE1$tipo_afiliacion <- factor(BASE1$tipo_afiliacion, levels = c(1,2,3), labels = c("Contributivo", "Subsidiado", "Otro"))
BASE1$vacuna_vph <- factor(BASE1$vacuna_vph, levels = c(1,2), labels = c("Si", "No"))
BASE1$resultado_vph <- replace(BASE1$resultado_vph, BASE1$resultado_vph == 3, NA)
BASE1$resultado_vph <- factor(BASE1$resultado_vph, levels = c(1,2), labels = c("Detectado", "No detectado"))

library(gtsummary)
tabla1 <- BASE1 %>%
  select(edad, tipo_afiliacion, vacuna_vph, meses_lastscreen, n_partos, resultado_vph) %>%
  tbl_summary(
    by = resultado_vph,
    type = list(
      all_continuous()                    ~ "continuous2",
      all_categorical()                   ~ "categorical"
    ),
    label = list(
      edad ~ "Edad", 
      tipo_afiliacion ~ "Afiliacion", 
      vacuna_vph ~"Vacunacion vph", 
      meses_lastscreen ~ "Tiempo en meses último tamizaje", 
      n_partos ~ "Paridad"
          ),
    statistic = list(
      all_continuous2()  ~ c(
        "{mean} ({sd})",
        "{median} ({p25}, {p75})",
        "{min} - {max}"
      ),
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = list(
      all_continuous()  ~ 2,
      all_categorical() ~ 2
    ),
    missing = "ifany",
    missing_text = "Sin dato"
  ) %>%
  add_p(
    test = list(
      all_continuous()  ~ "kruskal.test",  
      all_categorical() ~ "chisq.test"
    ),
    pvalue_fun = ~ style_pvalue(.x, digits = 3)
  ) %>%
  add_overall() %>%
  modify_header(label ~ "**Variable**") %>%
  bold_labels() %>%
  bold_p(t = 0.05) %>%
  modify_caption("**Tabla 1**")

tabla1
Tabla 1
Variable Overall
N = 869
1
Detectado
N = 226
1
No detectado
N = 643
1
p-value2
Edad


0.263
    Mean (SD) 43.91 (9.87) 43.42 (10.24) 44.09 (9.73)
    Median (Q1, Q3) 42.00 (36.00, 51.00) 41.00 (35.00, 50.00) 43.00 (36.00, 51.00)
    Min - Max 24.00 - 77.00 25.00 - 77.00 24.00 - 69.00
    Sin dato 3 0 3
Afiliacion


0.030
    Contributivo 129.00 (15.58%) 28.00 (13.02%) 101.00 (16.48%)
    Subsidiado 697.00 (84.18%) 185.00 (86.05%) 512.00 (83.52%)
    Otro 2.00 (0.24%) 2.00 (0.93%) 0.00 (0.00%)
    Sin dato 41 11 30
Vacunacion vph


0.731
    Si 98.00 (13.80%) 27.00 (14.84%) 71.00 (13.45%)
    No 612.00 (86.20%) 155.00 (85.16%) 457.00 (86.55%)
    Sin dato 159 44 115
Tiempo en meses último tamizaje


0.168
    Mean (SD) 49.45 (435.50) 32.26 (28.32) 55.37 (504.69)
    Median (Q1, Q3) 25.82 (16.89, 39.39) 25.07 (13.96, 37.06) 25.89 (17.74, 40.90)
    Min - Max 0.36 - 12,015.80 0.56 - 234.09 0.36 - 12,015.80
    Sin dato 108 31 77
Paridad


0.214
    Mean (SD) 2.39 (1.47) 2.31 (1.52) 2.41 (1.45)
    Median (Q1, Q3) 2.00 (2.00, 3.00) 2.00 (1.00, 3.00) 2.00 (2.00, 3.00)
    Min - Max 0.00 - 10.00 0.00 - 8.00 0.00 - 10.00
    Sin dato 4 1 3
1 n (%)
2 Kruskal-Wallis rank sum test; Pearson’s Chi-squared test

Variables aceptabilidad

BASE2 <- BASE %>% select(encuesta_post_toma,    como_se_sintio, razones_no_satisfaccion___1,    razones_no_satisfaccion___2,    razones_no_satisfaccion___3,    razones_no_satisfaccion___4,    razones_no_satisfaccion___5,    otra_no_satisfaccion,   confianza_prueba,   toma_prueba)
#lapply(BASE2, table)

BASE2$emocion <- ifelse(
  BASE2$razones_no_satisfaccion___1 == 1 | 
    BASE2$razones_no_satisfaccion___2 == 1 |
    grepl("inseguridad", BASE2$otra_no_satisfaccion, ignore.case = TRUE),
    2, 1)

BASE2$tecnico <- ifelse(
  BASE2$razones_no_satisfaccion___3 == 1 | 
  BASE2$razones_no_satisfaccion___4 == 1 |
  grepl("instrumento", BASE2$otra_no_satisfaccion, ignore.case = TRUE),
  2, 1
)
BASE2$confianza <- ifelse(BASE2$confianza_prueba == 1, 1, 2)

BASE2$acept_autotoma <- ifelse(BASE2$toma_prueba == 1, 1,2)

Análisis descriptivo

library(dplyr)
library(ggplot2)

BASE3 <- as.data.frame(BASE[, c("resultado_vph", "como_se_sintio", 
                  "razones_no_satisfaccion___1", "razones_no_satisfaccion___2",
                  "razones_no_satisfaccion___3", "razones_no_satisfaccion___4",
                  "razones_no_satisfaccion___5", "otra_no_satisfaccion",
                  "confianza_prueba", "toma_prueba")])

BASE3$como_se_sintio <- factor(BASE3$como_se_sintio, levels = c(1,2), labels = c("Bien", "No Bien"))

BASE3$resultado_vph <- replace(BASE3$resultado_vph, BASE3$resultado_vph == 3, NA)
BASE3$resultado_vph <- factor(BASE3$resultado_vph, levels = c(1,2), labels = c("Detectado", "No detectado"))

BASE3$barrera_emo <- factor(ifelse(
  BASE3$razones_no_satisfaccion___1 == 1 | 
    BASE3$razones_no_satisfaccion___2 == 1 |
    grepl("inseguridad", BASE3$otra_no_satisfaccion, ignore.case = TRUE),
    "Si", "No"))

BASE3$barrera_tec <- factor(ifelse(
  BASE3$razones_no_satisfaccion___3 == 1 | 
  BASE3$razones_no_satisfaccion___4 == 1 |
  grepl("instrumento", BASE3$otra_no_satisfaccion, ignore.case = TRUE),
  "Si", "No"
))

BASE3$confianza_prueba <- ifelse(BASE2$confianza_prueba == 1, 1, 2)
BASE3$confianza_prueba <- factor(BASE3$confianza_prueba, levels = c(1,2), labels = c("Si", "No"))
BASE3$toma_prueba <- factor(BASE3$toma_prueba, levels = c(1,2,3), labels = c("Misma mujer", "Profesional salud", "Sin preferencia"))

BASE3 <- BASE3 %>%
 filter(!is.na(toma_prueba))  #870 --> 856

BASE3 <- BASE3 %>%
 filter(!is.na(resultado_vph))  #856 --> 855


#BASE2$consumo_diario_std_3cat[is.na(BASE1$consumo_diario_std_3cat)] <- "Consumo Ligero"

datos_sankey <- as.data.frame(BASE3[, c("resultado_vph", "como_se_sintio", "barrera_emo", "barrera_tec", "confianza_prueba", "toma_prueba")])
#dput(names(BASE3))


library(forcats)

datos_sankey <- datos_sankey %>%
  rename(VPH = resultado_vph, Sintio_toma = como_se_sintio, Barrera_emocional = barrera_emo, Barrera_tecnica = barrera_tec, Confianza = confianza_prueba, Preferencia_prueba = toma_prueba)

BASE_transformado <- datos_sankey %>% group_by(VPH, Sintio_toma, Barrera_emocional, Barrera_tecnica, Confianza, Preferencia_prueba) %>% summarise(Freq = n()) %>% ungroup() 

library(easyalluvial) 

alluvial_wide( data = datos_sankey, 
               max_variables = 8, fill_by = 'last_variable',
               col_vector_flow = c("skyblue", "salmon", "#52b788"))

library(gt)
BASE_transformado %>%
  arrange(desc(Freq)) %>%       
  gt() %>%
  tab_header(title = "Conteos") %>%
  cols_label_with(fn = toupper)
Conteos
VPH SINTIO_TOMA BARRERA_EMOCIONAL BARRERA_TECNICA CONFIANZA PREFERENCIA_PRUEBA FREQ
No detectado Bien No No Si Misma mujer 297
Detectado Bien No No Si Misma mujer 87
No detectado Bien No No Si Profesional salud 87
No detectado Bien No No Si Sin preferencia 81
Detectado Bien No No Si Profesional salud 45
No detectado No Bien Si No Si Profesional salud 34
No detectado Bien No No No Profesional salud 32
Detectado Bien No No Si Sin preferencia 25
No detectado Bien Si No Si Profesional salud 17
No detectado No Bien Si No No Profesional salud 17
Detectado Bien No No No Profesional salud 11
No detectado Bien No No No Misma mujer 11
No detectado Bien Si No Si Misma mujer 11
No detectado No Bien Si No Si Misma mujer 9
Detectado No Bien Si No No Profesional salud 8
Detectado No Bien Si No Si Profesional salud 7
No detectado No Bien No Si Si Profesional salud 6
No detectado No Bien Si No Si Sin preferencia 5
No detectado No Bien Si Si No Profesional salud 5
Detectado Bien Si No Si Misma mujer 4
Detectado Bien Si No No Profesional salud 4
Detectado No Bien Si No Si Misma mujer 4
No detectado Bien Si No Si Sin preferencia 4
No detectado No Bien Si No No Misma mujer 4
No detectado Bien No No No Sin preferencia 3
Detectado Bien No No No Misma mujer 2
Detectado Bien Si No Si Profesional salud 2
Detectado No Bien Si No Si Sin preferencia 2
Detectado No Bien Si No No Misma mujer 2
Detectado No Bien Si Si No Profesional salud 2
No detectado Bien No Si Si Misma mujer 2
No detectado Bien No Si No Profesional salud 2
No detectado Bien Si No No Misma mujer 2
No detectado Bien Si No No Profesional salud 2
No detectado No Bien No No No Profesional salud 2
No detectado No Bien No Si No Profesional salud 2
No detectado No Bien Si Si Si Misma mujer 2
Detectado Bien No No No Sin preferencia 1
Detectado Bien Si No Si Sin preferencia 1
Detectado No Bien No No Si Profesional salud 1
Detectado No Bien No No No Profesional salud 1
Detectado No Bien No Si No Profesional salud 1
Detectado No Bien Si Si Si Misma mujer 1
Detectado No Bien Si Si Si Profesional salud 1
No detectado Bien Si No No Sin preferencia 1
No detectado Bien Si Si Si Misma mujer 1
No detectado No Bien No No Si Sin preferencia 1
No detectado No Bien No Si Si Misma mujer 1
No detectado No Bien Si No No Sin preferencia 1
No detectado No Bien Si Si Si Profesional salud 1
tabla_acept <- BASE3 %>%
  select(como_se_sintio, barrera_emo, barrera_tec, 
         confianza_prueba, toma_prueba, resultado_vph) %>%
  tbl_summary(
    by = resultado_vph,
    label = list(
      como_se_sintio ~ "¿Cómo se sintió?",
      barrera_emo ~ "Barrera emocional",
      barrera_tec ~ "Barrera técnica",
      confianza_prueba ~ "Confianza en la prueba",
      toma_prueba ~ "Preferencia de toma"
    ),
    statistic = all_categorical() ~ "{n} ({p}%)",
    missing = "ifany",
    missing_text = "Sin dato"
  ) %>%
  add_p(pvalue_fun = ~ style_pvalue(.x, digits = 3)) %>%
  add_overall() %>%
  bold_labels() %>%
  bold_p(t = 0.05) %>%
  modify_caption("**Aceptabilidad de autotoma según resultado VPH**")
tabla_acept
Aceptabilidad de autotoma según resultado VPH
Characteristic Overall
N = 855
1
Detectado
N = 212
1
No detectado
N = 643
1
p-value2
¿Cómo se sintió?


0.955
    Bien 735 (86%) 182 (86%) 553 (86%)
    No Bien 120 (14%) 30 (14%) 90 (14%)
Barrera emocional


0.970
    No 701 (82%) 174 (82%) 527 (82%)
    Si 154 (18%) 38 (18%) 116 (18%)
Barrera técnica


0.443
    No 828 (97%) 207 (98%) 621 (97%)
    Si 27 (3.2%) 5 (2.4%) 22 (3.4%)
Confianza en la prueba


0.454
    Si 739 (86%) 180 (85%) 559 (87%)
    No 116 (14%) 32 (15%) 84 (13%)
Preferencia de toma


0.178
    Misma mujer 440 (51%) 100 (47%) 340 (53%)
    Profesional salud 290 (34%) 83 (39%) 207 (32%)
    Sin preferencia 125 (15%) 29 (14%) 96 (15%)
1 n (%)
2 Pearson’s Chi-squared test

Analisis clase latente

library(poLCA)
set.seed(2026)

Class <- cbind(como_se_sintio, emocion, confianza, acept_autotoma) ~ 1

lca_2 <- poLCA(Class, data = BASE2, nclass = 2, nrep = 10, maxiter = 1000, graphs = FALSE)
## Model 1: llik = -1440.129 ... best llik = -1440.129
## Model 2: llik = -1440.129 ... best llik = -1440.129
## Model 3: llik = -1440.129 ... best llik = -1440.129
## Model 4: llik = -1440.129 ... best llik = -1440.129
## Model 5: llik = -1440.129 ... best llik = -1440.129
## Model 6: llik = -1440.129 ... best llik = -1440.129
## Model 7: llik = -1440.129 ... best llik = -1440.129
## Model 8: llik = -1440.129 ... best llik = -1440.129
## Model 9: llik = -1440.129 ... best llik = -1440.129
## Model 10: llik = -1440.129 ... best llik = -1440.129
## Conditional item response (column) probabilities,
##  by outcome variable, for each class (row) 
##  
## $como_se_sintio
##           Pr(1) Pr(2)
## class 1:  0.191 0.809
## class 2:  1.000 0.000
## 
## $emocion
##            Pr(1)  Pr(2)
## class 1:  0.1309 0.8691
## class 2:  0.9646 0.0354
## 
## $confianza
##            Pr(1)  Pr(2)
## class 1:  0.6320 0.3680
## class 2:  0.9132 0.0868
## 
## $acept_autotoma
##            Pr(1)  Pr(2)
## class 1:  0.1856 0.8144
## class 2:  0.5843 0.4157
## 
## Estimated class population shares 
##  0.1733 0.8267 
##  
## Predicted class memberships (by modal posterior prob.) 
##  0.1787 0.8213 
##  
## ========================================================= 
## Fit for 2 latent classes: 
## ========================================================= 
## number of observations: 856 
## number of estimated parameters: 9 
## residual degrees of freedom: 6 
## maximum log-likelihood: -1440.129 
##  
## AIC(2): 2898.258
## BIC(2): 2941.028
## G^2(2): 43.33248 (Likelihood ratio/deviance statistic) 
## X^2(2): 42.29904 (Chi-square goodness of fit) 
## 
lca_3 <- poLCA(Class, data = BASE2, nclass = 3, nrep = 10, maxiter = 1000, graphs = FALSE)
## Model 1: llik = -1420.693 ... best llik = -1420.693
## Model 2: llik = -1420.692 ... best llik = -1420.692
## Model 3: llik = -1420.694 ... best llik = -1420.692
## Model 4: llik = -1419.552 ... best llik = -1419.552
## Model 5: llik = -1420.694 ... best llik = -1419.552
## Model 6: llik = -1420.694 ... best llik = -1419.552
## Model 7: llik = -1420.694 ... best llik = -1419.552
## Model 8: llik = -1420.693 ... best llik = -1419.552
## Model 9: llik = -1419.551 ... best llik = -1419.551
## Model 10: llik = -1419.552 ... best llik = -1419.551
## Conditional item response (column) probabilities,
##  by outcome variable, for each class (row) 
##  
## $como_se_sintio
##            Pr(1)  Pr(2)
## class 1:  0.2111 0.7889
## class 2:  0.8729 0.1271
## class 3:  0.9975 0.0025
## 
## $emocion
##            Pr(1)  Pr(2)
## class 1:  0.0000 1.0000
## class 2:  0.9943 0.0057
## class 3:  0.9668 0.0332
## 
## $confianza
##            Pr(1)  Pr(2)
## class 1:  0.6448 0.3552
## class 2:  0.5419 0.4581
## class 3:  0.9677 0.0323
## 
## $acept_autotoma
##            Pr(1)  Pr(2)
## class 1:  0.2037 0.7963
## class 2:  0.0000 1.0000
## class 3:  0.6716 0.3284
## 
## Estimated class population shares 
##  0.1553 0.1247 0.72 
##  
## Predicted class memberships (by modal posterior prob.) 
##  0.1612 0.0736 0.7652 
##  
## ========================================================= 
## Fit for 3 latent classes: 
## ========================================================= 
## number of observations: 856 
## number of estimated parameters: 14 
## residual degrees of freedom: 1 
## maximum log-likelihood: -1419.551 
##  
## AIC(3): 2867.102
## BIC(3): 2933.634
## G^2(3): 2.17683 (Likelihood ratio/deviance statistic) 
## X^2(3): 2.134992 (Chi-square goodness of fit) 
##  
## ALERT: iterations finished, MAXIMUM LIKELIHOOD NOT FOUND 
## 
library(tidyverse)

fit_table <- tibble(
  Clases = c(2, 3),
  LogLik = c(lca_2$llik, lca_3$llik),
  AIC    = c(lca_2$aic,  lca_3$aic),
  BIC    = c(lca_2$bic,  lca_3$bic),
  Chi2   = c(lca_2$Chisq, lca_3$Chisq),
  df     = c(lca_2$resid.df, lca_3$resid.df)
)

library(gt)

fit_table %>%
  gt() %>%
  tab_header(title = "Comparación de modelos LCA") %>%
  fmt_number(columns = c(LogLik, AIC, BIC, Chi2), decimals = 2) %>%
  tab_footnote("Menor AIC y BIC indican mejor ajuste. Entropía > 0.80 indica buena clasificación.") %>%
  tab_style(
    style = cell_fill(color = "#d4edda"),
    locations = cells_body(rows = BIC == min(BIC))
  )
Comparación de modelos LCA
Clases LogLik AIC BIC Chi2 df
2 −1,440.13 2,898.26 2,941.03 42.30 6
3 −1,419.55 2,867.10 2,933.63 2.13 1
Menor AIC y BIC indican mejor ajuste. Entropía > 0.80 indica buena clasificación.