Carga de base

Carga de librerías

library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.4.3
## ── 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   4.0.0     ✔ 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(knitr)
library(tableone)
library(lme4)
## Cargando paquete requerido: Matrix
## 
## Adjuntando el paquete: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library(sjPlot)
## #refugeeswelcome
## 
## Adjuntando el paquete: 'sjPlot'
## 
## The following object is masked from 'package:ggplot2':
## 
##     set_theme
library(performance)
## Warning: package 'performance' was built under R version 4.4.3
library(lmerTest)
## Warning: package 'lmerTest' was built under R version 4.4.3
## 
## Adjuntando el paquete: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step
#library(rworldmap) # no es necesario 
#library(countrycode) # no es necesario
library(dplyr)
library(performance)
library(DHARMa)
## Warning: package 'DHARMa' was built under R version 4.4.3
## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(ggeffects)
library(naniar)
## Warning: package 'naniar' was built under R version 4.4.3
library(marginaleffects)
## Warning: package 'marginaleffects' was built under R version 4.4.3
library(data.table)
## Warning: package 'data.table' was built under R version 4.4.3
## 
## Adjuntando el paquete: 'data.table'
## 
## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## 
## The following object is masked from 'package:purrr':
## 
##     transpose
library(haven)
library(lme4)
library(geepack)
## Warning: package 'geepack' was built under R version 4.4.3
library(performance)
library(psych)
## 
## Adjuntando el paquete: 'psych'
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(sjstats)
## 
## Adjuntando el paquete: 'sjstats'
## 
## The following object is masked from 'package:psych':
## 
##     phi
## 
## The following objects are masked from 'package:performance':
## 
##     icc, r2
options(scipen=999, digits=2)

Datos faltantes

vis_miss(TES)

Eliminar datos faltantes de las columnas de caract del pie

TES <- TES %>%
  filter(!if_any(16:25, is.na))

Pasar variables nuéricas a factor

TES[sapply(TES, is.character)] <- lapply(TES[sapply(TES, is.character)], as.factor)
TES$Cicatrización <- as.factor(TES$Cicatrización)


TES <- TES %>%
  mutate(across(16:25, as.factor))
TES <- TES %>%
  mutate(across(5:13, as.factor))
TES <- TES %>%
  mutate(across(OM, as.factor))

Evaluar cuántos pacientes cicatrizaron

table1<-table(TES$Cicatrización)
table1
## 
##   0   1 
## 124 175
prop.table(table1)
## 
##    0    1 
## 0.41 0.59

Distribución de variables según outcomes

TES$cicat=TES$Cicatrización
#Isquemia
TES%>%
  dplyr::group_by(isq) %>%
  dplyr::summarise(tar = mean (as.numeric(cicat)-1)) %>%
  ggplot(aes(x = isq, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", col = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## `geom_smooth()` using formula = 'y ~ x'

TES %>%
  dplyr::group_by(isq) %>%
  dplyr::summarise(tar = mean(as.numeric(cicat)-1)) %>%
  ggplot(aes(x = isq, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "isquemia", y = "Proporción de cicatrizados") +
  theme_minimal()

#localizacion
TES%>%
  dplyr::group_by(loca) %>%
  dplyr::summarise(tar = mean (as.numeric(cicat)-1)) %>%
  ggplot(aes(x = loca, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", col = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## `geom_smooth()` using formula = 'y ~ x'

TES %>%
  dplyr::group_by(loca) %>%
  dplyr::summarise(tar = mean(as.numeric(cicat)-1)) %>%
  ggplot(aes(x = loca, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "localizacion", y = "Proporción de cicatrizados") +
  theme_minimal()

#topo
TES%>%
  dplyr::group_by(topo) %>%
  dplyr::summarise(tar = mean (as.numeric(cicat)-1)) %>%
  ggplot(aes(x = topo, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", col = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## `geom_smooth()` using formula = 'y ~ x'

TES %>%
  dplyr::group_by(topo) %>%
  dplyr::summarise(tar = mean(as.numeric(cicat)-1)) %>%
  ggplot(aes(x = topo, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "topografia", y = "Proporción de cicatrizados") +
  theme_minimal()

#numero
TES%>%
  dplyr::group_by(nume) %>%
  dplyr::summarise(tar = mean (as.numeric(cicat)-1)) %>%
  ggplot(aes(x = nume, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", col = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## `geom_smooth()` using formula = 'y ~ x'

TES %>%
  dplyr::group_by(nume) %>%
  dplyr::summarise(tar = mean(as.numeric(cicat)-1)) %>%
  ggplot(aes(x = nume, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "numero de zonas", y = "Proporción de cicatrizados") +
  theme_minimal()

#infeccion
TES%>%
  dplyr::group_by(infe) %>%
  dplyr::summarise(tar = mean (as.numeric(cicat)-1)) %>%
  ggplot(aes(x = infe, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", col = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## `geom_smooth()` using formula = 'y ~ x'

TES %>%
  dplyr::group_by(infe) %>%
  dplyr::summarise(tar = mean(as.numeric(cicat)-1)) %>%
  ggplot(aes(x = infe, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "infeccion", y = "Proporción de cicatrizados") +
  theme_minimal()

#edema
TES%>%
  dplyr::group_by(ede) %>%
  dplyr::summarise(tar = mean (as.numeric(cicat)-1)) %>%
  ggplot(aes(x = ede, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", col = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## `geom_smooth()` using formula = 'y ~ x'

TES %>%
  dplyr::group_by(ede) %>%
  dplyr::summarise(tar = mean(as.numeric(cicat)-1)) %>%
  ggplot(aes(x = ede, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "edema", y = "Proporción de cicatrizados") +
  theme_minimal()

#neuropatia
TES%>%
  dplyr::group_by(neur) %>%
  dplyr::summarise(tar = mean (as.numeric(cicat)-1)) %>%
  ggplot(aes(x = neur, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", col = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## `geom_smooth()` using formula = 'y ~ x'

TES %>%
  dplyr::group_by(neur) %>%
  dplyr::summarise(tar = mean(as.numeric(cicat)-1)) %>%
  ggplot(aes(x = neur, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "Neuropatia", y = "Proporción de cicatrizados") +
  theme_minimal()

#tisu profundidad
TES%>%
  dplyr::group_by(tisu) %>%
  dplyr::summarise(tar = mean (as.numeric(cicat)-1)) %>%
  ggplot(aes(x = tisu, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", col = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## `geom_smooth()` using formula = 'y ~ x'

TES %>%
  dplyr::group_by(tisu) %>%
  dplyr::summarise(tar = mean(as.numeric(cicat)-1)) %>%
  ggplot(aes(x = tisu, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "Profundidad", y = "Proporción de cicatrizados") +
  theme_minimal()

#area
TES%>%
  dplyr::group_by(area) %>%
  dplyr::summarise(tar = mean (as.numeric(cicat)-1)) %>%
  ggplot(aes(x = area, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", col = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## `geom_smooth()` using formula = 'y ~ x'

TES %>%
  dplyr::group_by(area) %>%
  dplyr::summarise(tar = mean(as.numeric(cicat)-1)) %>%
  ggplot(aes(x = area, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "Area", y = "Proporción de cicatrizados") +
  theme_minimal()

#cicat
TES%>%
  dplyr::group_by(cica) %>%
  dplyr::summarise(tar = mean (as.numeric(cicat)-1)) %>%
  ggplot(aes(x = cica, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", col = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## `geom_smooth()` using formula = 'y ~ x'

TES %>%
  dplyr::group_by(cica) %>%
  dplyr::summarise(tar = mean(as.numeric(cicat)-1)) %>%
  ggplot(aes(x = cica, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "Fase de cicatrizacion", y = "Proporción de cicatrizados") +
  theme_minimal()

Categorización de las variables

#ISQuEMIA
TES <- TES %>%
  mutate(isq_factor = factor(if_else(isq %in% c(2, 3), "SI", "NO/leve")),
         isq_factor = factor(isq_factor, levels = rev(levels(isq_factor))))

table(TES$isq)
## 
##   0   1   2   3 
## 157  56  50  36
table(TES$isq_factor)
## 
##      SI NO/leve 
##      86     213
#localizacion
TES <- TES %>%
  mutate(loca_factor = factor(if_else(loca %in% c(2, 3), "Met/Tar", "Fal")),
         loca_factor = factor(loca_factor, levels = rev(levels(loca_factor))))


table(TES$loca)
## 
##   1   2   3 
## 166  99  34
table(TES$loca_factor)
## 
## Met/Tar     Fal 
##     133     166
#topog
TES <- TES %>%
  mutate(topo_factor = factor(if_else(topo %in% c(1,2),  "dor/plan/lat", "Mas de dos")))
TES <- TES %>%
  mutate(topo_factor = factor(if_else(topo %in% c(1, 2), "dor/plan/lat", "Mas de dos")),
         topo_factor = factor(topo_factor, levels = rev(levels(topo_factor))))

table(TES$topo)
## 
##   1   2   3 
## 136  60 103
table(TES$topo_factor)
## 
##   Mas de dos dor/plan/lat 
##          103          196
#nume
TES <- TES %>%
  mutate(nume_factor = factor(if_else(nume %in% c(2, 3),  "Dos o mas", "Una")))

table(TES$nume)
## 
##   1   2   3 
## 211  59  29
table(TES$nume_factor)
## 
## Dos o mas       Una 
##        88       211
#infe
TES <- TES %>%
  mutate(infe_factor = case_when(
    infe == 0 ~ "No",
    infe %in% c(1, 2) ~ "Leve/Mod",
    infe == 3 ~ "Grave",
    TRUE ~ NA_character_     # conserva los NA originales
  )) %>%
  mutate(infe_factor = factor(infe_factor, levels = c("No", "Leve/Mod", "Grave")))
TES$infe_factor <- factor (TES$infe_factor, levels= c ("Grave","Leve/Mod","No"))

table(TES$infe)
## 
##   0   1   2   3 
##  67  43 152  37
table(TES$infe_factor)
## 
##    Grave Leve/Mod       No 
##       37      195       67
#edema
TES <- TES %>%
  mutate(ede_factor = case_when(
    ede == 2 ~ "Unilateral",
    ede %in% c(0, 1) ~ "Sin/local",
    ede == 3 ~ "Bilateral",
    TRUE ~ NA_character_     # conserva los NA originales
  )) %>%
  mutate(ede_factor = factor(ede_factor, levels = c("Sin/local", "Unilateral", "Bilateral")))
TES$ede_factor <- factor (TES$ede_factor, levels= c ("Bilateral","Unilateral","Sin/local"))
table(TES$ede)
## 
##  0  1  2  3 
## 92 98 89 20
table(TES$ede_factor)
## 
##  Bilateral Unilateral  Sin/local 
##         20         89        190
#neur
TES <- TES %>%
  mutate(neur_factor = factor(
    if_else(neur %in% c(1, 2, 3), "Leve/mod/charcot", "No"),
    levels = c("No", "Leve/mod/charcot")
  ))
table(TES$neur)
## 
##   0   1   2   3 
##  10  57 216  16
table(TES$neur_factor)
## 
##               No Leve/mod/charcot 
##               10              289
#tisu (prof)
TES <- TES %>%
  mutate(tisu_factor = factor(if_else(tisu %in% c(2, 3),  "parcial/total", "Piel")))
table(TES$tisu)
## 
##   1   2   3 
##  52  98 149
table(TES$tisu_factor)
## 
## parcial/total          Piel 
##           247            52
#area
TES$area_factor <- with(TES, ifelse(area == 1, "Menor a 10",
                             ifelse(area %in% c(1, 2), "10 a 40",
                             ifelse(area == 3, "Mas de 40", NA))))
TES$area_factor <- factor(TES$area_factor, levels = c("Menor a 10", "10 a 40", "Mas de 40"))

table(TES$area)
## 
##   1   2   3 
## 169 103  27
table(TES$area_factor)
## 
## Menor a 10    10 a 40  Mas de 40 
##        169        103         27
#cica(fase)
TES <- TES %>%
  mutate(cica_factor = factor(if_else(cica %in% c(2, 3),  "Granu/infla", "Epit")))
TES$cica_factor <- factor (TES$cica_factor, levels= c ("Granu/infla","Epit"))
table(TES$cica)
## 
##   1   2   3 
##  14  48 237
table(TES$cica_factor)
## 
## Granu/infla        Epit 
##         285          14

Tabla 1

#Creo variable muerte 1 si/no para que la tabla se vea mejor
TES$cicat1 <- factor(TES$cicat,
                       levels = c(0, 1),
                       labels = c("No cicatrizo", "Cicatrizo"))

catvars = c( "isq_factor", "loca_factor", "topo_factor", "nume_factor", "infe_factor", "ede_factor", "neur_factor", "tisu_factor", "area_factor", "cica_factor", "OM","tbq", "aamen", "aamay", "dial", "ICC", "HTA", "CV", "sexo1fem")


vars = c("isq_factor", "loca_factor", "topo_factor", "nume_factor", "infe_factor", "ede_factor", "neur_factor", "tisu_factor", "area_factor", "cica_factor", "OM", "sanelian","SINBAD","tbq", "aamen", "aamay", "dial", "ICC", "HTA", "CV","Edad", "sexo1fem" )


tabla_1a <- CreateTableOne(vars = vars, strata = "cicat1", factorVars = catvars, data = TES)

print(tabla_1a)
##                                     Stratified by cicat1
##                                      No cicatrizo  Cicatrizo     p      test
##   n                                    124           175                    
##   isq_factor = NO/leve (%)              70 (56.5)    143 (81.7)  <0.001     
##   loca_factor = Fal (%)                 59 (47.6)    107 (61.1)   0.027     
##   topo_factor = dor/plan/lat (%)        70 (56.5)    126 (72.0)   0.008     
##   nume_factor = Una (%)                 62 (50.0)    149 (85.1)  <0.001     
##   infe_factor (%)                                                 0.002     
##      Grave                              23 (18.5)     14 ( 8.0)             
##      Leve/Mod                           83 (66.9)    112 (64.0)             
##      No                                 18 (14.5)     49 (28.0)             
##   ede_factor (%)                                                 <0.001     
##      Bilateral                          14 (11.3)      6 ( 3.4)             
##      Unilateral                         46 (37.1)     43 (24.6)             
##      Sin/local                          64 (51.6)    126 (72.0)             
##   neur_factor = Leve/mod/charcot (%)   123 (99.2)    166 (94.9)   0.084     
##   tisu_factor = Piel (%)                 9 ( 7.3)     43 (24.6)  <0.001     
##   area_factor (%)                                                 0.002     
##      Menor a 10                         57 (46.0)    112 (64.0)             
##      10 a 40                            49 (39.5)     54 (30.9)             
##      Mas de 40                          18 (14.5)      9 ( 5.1)             
##   cica_factor = Epit (%)                 1 ( 0.8)     13 ( 7.4)   0.017     
##   OM = 1 (%)                            65 (52.4)     74 (42.3)   0.107     
##   sanelian (mean (SD))               18.59 (4.03)  15.49 (3.96)  <0.001     
##   SINBAD (mean (SD))                  4.40 (1.12)   3.62 (1.15)  <0.001     
##   tbq (%)                                                         0.414     
##      0                                   3 ( 2.4)     11 ( 6.3)             
##      1                                  29 (23.4)     34 (19.5)             
##      2                                  61 (49.2)     86 (49.4)             
##      3                                  31 (25.0)     43 (24.7)             
##   aamen = 1 (%)                         40 (32.3)     41 (23.4)   0.119     
##   aamay = 1 (%)                         11 ( 8.9)      6 ( 3.4)   0.080     
##   dial = 1 (%)                          13 (10.5)     13 ( 7.4)   0.474     
##   ICC = 1 (%)                           20 (16.1)     20 (11.4)   0.315     
##   HTA = 1 (%)                           95 (76.6)    120 (68.6)   0.163     
##   CV = 1 (%)                            30 (24.2)     43 (24.6)   1.000     
##   Edad (mean (SD))                   59.00 (12.98) 56.98 (11.36)  0.154     
##   sexo1fem = 1 (%)                      29 (23.4)     36 (20.6)   0.660
kableone(tabla_1a)
No cicatrizo Cicatrizo p test
n 124 175
isq_factor = NO/leve (%) 70 (56.5) 143 (81.7) <0.001
loca_factor = Fal (%) 59 (47.6) 107 (61.1) 0.027
topo_factor = dor/plan/lat (%) 70 (56.5) 126 (72.0) 0.008
nume_factor = Una (%) 62 (50.0) 149 (85.1) <0.001
infe_factor (%) 0.002
Grave 23 (18.5) 14 ( 8.0)
Leve/Mod 83 (66.9) 112 (64.0)
No 18 (14.5) 49 (28.0)
ede_factor (%) <0.001
Bilateral 14 (11.3) 6 ( 3.4)
Unilateral 46 (37.1) 43 (24.6)
Sin/local 64 (51.6) 126 (72.0)
neur_factor = Leve/mod/charcot (%) 123 (99.2) 166 (94.9) 0.084
tisu_factor = Piel (%) 9 ( 7.3) 43 (24.6) <0.001
area_factor (%) 0.002
Menor a 10 57 (46.0) 112 (64.0)
10 a 40 49 (39.5) 54 (30.9)
Mas de 40 18 (14.5) 9 ( 5.1)
cica_factor = Epit (%) 1 ( 0.8) 13 ( 7.4) 0.017
OM = 1 (%) 65 (52.4) 74 (42.3) 0.107
sanelian (mean (SD)) 18.59 (4.03) 15.49 (3.96) <0.001
SINBAD (mean (SD)) 4.40 (1.12) 3.62 (1.15) <0.001
tbq (%) 0.414
0 3 ( 2.4) 11 ( 6.3)
1 29 (23.4) 34 (19.5)
2 61 (49.2) 86 (49.4)
3 31 (25.0) 43 (24.7)
aamen = 1 (%) 40 (32.3) 41 (23.4) 0.119
aamay = 1 (%) 11 ( 8.9) 6 ( 3.4) 0.080
dial = 1 (%) 13 (10.5) 13 ( 7.4) 0.474
ICC = 1 (%) 20 (16.1) 20 (11.4) 0.315
HTA = 1 (%) 95 (76.6) 120 (68.6) 0.163
CV = 1 (%) 30 (24.2) 43 (24.6) 1.000
Edad (mean (SD)) 59.00 (12.98) 56.98 (11.36) 0.154
sexo1fem = 1 (%) 29 (23.4) 36 (20.6) 0.660

UNIVARIADO

TES$cicat=TES$Cicatrización
modelo1 <- glm(cicat ~ loca_factor+ topo_factor+nume_factor+ isq_factor+ infe_factor+ ede_factor+ neur_factor+ tisu_factor+ area_factor+ cica_factor , family = "binomial", data = TES)
t1<-tab_model(modelo1, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Cicatrización")
knitr::asis_output(t1$knitr)
Cicatrización
  cicat
Predictors Odds Ratios std. Error CI p
Met/Tar Reference
Fal 1.89 0.55 1.08 – 3.35 0.027
Mas de dos Reference
dor/plan/lat 0.94 0.33 0.47 – 1.85 0.852
Dos o mas Reference
Una 5.53 2.06 2.72 – 11.75 <0.001
SI Reference
NO/leve 3.44 1.07 1.89 – 6.38 <0.001
Grave Reference
Leve/Mod 2.31 1.07 0.94 – 5.85 0.070
No 2.56 1.53 0.80 – 8.41 0.115
Bilateral Reference
Unilateral 1.98 1.19 0.63 – 6.76 0.256
Sin/local 2.45 1.40 0.82 – 7.96 0.116
No Reference
Leve/mod/charcot 0.30 0.34 0.02 – 1.90 0.284
parcial/total Reference
Piel 1.89 0.90 0.76 – 4.99 0.183
Menor a 10 Reference
10 a 40 1.75 0.64 0.87 – 3.66 0.127
Mas de 40 1.96 1.23 0.57 – 6.88 0.286
Granu/infla Reference
Epit 4.18 4.60 0.68 – 80.84 0.194
Observations 299
Deviance 321.666
log-Likelihood -160.833 
modelo2 <- glm(cicat ~ loca_factor , family = "binomial", data = TES)
modelo3 <- glm(cicat ~  topo_factor , family = "binomial", data = TES)
modelo4 <- glm(cicat ~ nume_factor , family = "binomial", data = TES)
modelo5 <- glm(cicat ~  isq_factor , family = "binomial", data = TES)
modelo6<- glm(cicat ~  infe_factor, family = "binomial", data = TES)
modelo7<- glm(cicat ~  ede_factor , family = "binomial", data = TES)
modelo8<- glm(cicat ~  neur_factor , family = "binomial", data = TES)
modelo9<- glm(cicat ~  tisu_factor , family = "binomial", data = TES)
modelo10<- glm(cicat ~  area_factor , family = "binomial", data = TES)
modelo11<- glm(cicat ~  cica_factor, family = "binomial", data = TES)
modelo12<- glm(cicat ~  evol, family = "binomial", data = TES)

library(broom)
## 
## Adjuntando el paquete: 'broom'
## The following object is masked from 'package:sjstats':
## 
##     bootstrap
library(dplyr)

# Crear una lista con los modelos
modelos <- list(
  modelo2, modelo3, modelo4, modelo5,
  modelo6, modelo7, modelo8, modelo9,
  modelo10, modelo11, modelo12
)

# Convertir a tabla ordenada
tabla_coef <- lapply(modelos, tidy) %>%
  bind_rows(.id = "modelo")

# Ver la tabla
tabla_coef
## # A tibble: 25 × 6
##    modelo term                    estimate std.error statistic  p.value
##    <chr>  <chr>                      <dbl>     <dbl>     <dbl>    <dbl>
##  1 1      (Intercept)               0.0451     0.173     0.260 7.95e- 1
##  2 1      loca_factorFal            0.550      0.237     2.32  2.05e- 2
##  3 2      (Intercept)              -0.0972     0.197    -0.492 6.22e- 1
##  4 2      topo_factordor/plan/lat   0.685      0.247     2.77  5.61e- 3
##  5 3      (Intercept)              -0.869      0.234    -3.72  2.00e- 4
##  6 3      nume_factorUna            1.75       0.278     6.27  3.52e-10
##  7 4      (Intercept)              -0.523      0.223    -2.35  1.90e- 2
##  8 4      isq_factorNO/leve         1.24       0.267     4.64  3.43e- 6
##  9 5      (Intercept)              -0.496      0.339    -1.46  1.43e- 1
## 10 5      infe_factorLeve/Mod       0.796      0.369     2.16  3.08e- 2
## # ℹ 15 more rows
#ordenar p valor
tabla_coef_ordenada <- tabla_coef %>%
  arrange(p.value)
# Ver la tabla
tabla_coef_ordenada
## # A tibble: 25 × 6
##    modelo term                    estimate std.error statistic  p.value
##    <chr>  <chr>                      <dbl>     <dbl>     <dbl>    <dbl>
##  1 3      nume_factorUna             1.75      0.278      6.27 3.52e-10
##  2 4      isq_factorNO/leve          1.24      0.267      4.64 3.43e- 6
##  3 9      (Intercept)                0.675     0.163      4.15 3.30e- 5
##  4 3      (Intercept)               -0.869     0.234     -3.72 2.00e- 4
##  5 8      tisu_factorPiel            1.43      0.388      3.67 2.38e- 4
##  6 11     (Intercept)                0.511     0.141      3.63 2.82e- 4
##  7 5      infe_factorNo              1.50      0.437      3.43 6.07e- 4
##  8 9      area_factorMas de 40      -1.37      0.439     -3.11 1.84e- 3
##  9 6      ede_factorSin/local        1.52      0.512      2.98 2.88e- 3
## 10 2      topo_factordor/plan/lat    0.685     0.247      2.77 5.61e- 3
## # ℹ 15 more rows

evaluación del ICC

TES$cicat_num <- as.numeric(as.character(TES$cicat))

mod_aleatorio_NULO <- lmer(cicat_num ~ (1 | NOMBRE), REML = T, data = TES)


tab_model(mod_aleatorio_NULO)
  cicat_num
Predictors Estimates CI p
(Intercept) 0.60 0.51 – 0.68 <0.001
Random Effects
σ2 0.24
τ00 NOMBRE 0.01
ICC 0.04
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.000 / 0.037 
performance::icc(mod_aleatorio_NULO)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.037
##   Unadjusted ICC: 0.037
library(broom.mixed)
## Warning: package 'broom.mixed' was built under R version 4.4.3
tidy(mod_aleatorio_NULO)
## # A tibble: 3 × 8
##   effect   group    term            estimate std.error statistic    df   p.value
##   <chr>    <chr>    <chr>              <dbl>     <dbl>     <dbl> <dbl>     <dbl>
## 1 fixed    <NA>     (Intercept)       0.597     0.0418      14.3  5.63   1.24e-5
## 2 ran_pars NOMBRE   sd__(Intercept)   0.0961   NA           NA   NA     NA      
## 3 ran_pars Residual sd__Observation   0.488    NA           NA   NA     NA

MODELO MIXTO 1

mod_mixto <- glmer(cicat ~ nume_factor + (1 | NOMBRE),
                   family = binomial,
                   data = TES)
## boundary (singular) fit: see help('isSingular')
tab_model(mod_mixto)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.42 0.27 – 0.66 <0.001
nume factor [Una] 5.73 3.32 – 9.89 <0.001
Random Effects
σ2 3.29
τ00 NOMBRE 0.00
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.162 / NA

MODELO MIXTO 2

mod_mixto2 <- glmer(cicat ~ nume_factor + isq_factor + (1 | NOMBRE),
                   family = binomial,
                   data = TES)


tab_model(mod_mixto2)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.17 0.09 – 0.34 <0.001
nume factor [Una] 5.79 3.24 – 10.33 <0.001
isq factor [NO/leve] 3.46 1.95 – 6.13 <0.001
Random Effects
σ2 3.29
τ00 NOMBRE 0.03
ICC 0.01
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.242 / 0.248

MODELO MIXTO 3

mod_mixto3 <- glmer(cicat ~ nume_factor + isq_factor +infe_factor +(1 | NOMBRE),
                   family = binomial,
                   data = TES)


tab_model(mod_mixto3)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.08 0.03 – 0.24 <0.001
nume factor [Una] 5.12 2.81 – 9.30 <0.001
isq factor [NO/leve] 3.73 2.06 – 6.74 <0.001
infe factor [Leve/Mod] 2.19 0.90 – 5.36 0.086
infe factor [No] 3.31 1.16 – 9.47 0.025
Random Effects
σ2 3.29
τ00 NOMBRE 0.13
ICC 0.04
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.270 / 0.297

MODELO MIXTO 4

mod_mixto4 <- glmer(cicat ~ nume_factor + isq_factor + tisu_factor + infe_factor+(1 | NOMBRE),
                   family = binomial,
                   data = TES)


tab_model(mod_mixto4)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.09 0.03 – 0.27 <0.001
nume factor [Una] 5.02 2.76 – 9.13 <0.001
isq factor [NO/leve] 3.36 1.84 – 6.12 <0.001
tisu factor [Piel] 2.26 0.91 – 5.61 0.079
infe factor [Leve/Mod] 1.96 0.82 – 4.71 0.130
infe factor [No] 2.21 0.73 – 6.71 0.161
Random Effects
σ2 3.29
τ00 NOMBRE 0.10
ICC 0.03
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.285 / 0.306

MATRIZ DE CORRELACION PARA EVALUAR COLINEALIDAD TISU INFE. Impresiona confusión

library(dplyr)

TES %>%
  select(tisu_factor, infe_factor) %>%
  mutate(across(everything(), ~as.numeric(as.factor(.)))) %>%
  cor()
##             tisu_factor infe_factor
## tisu_factor        1.00        0.45
## infe_factor        0.45        1.00
#DA COlinealidad MODERADA

MODELO MIXTO 5

mod_mixto5 <- glmer(cicat ~ nume_factor + isq_factor + tisu_factor + infe_factor+ ede_factor+(1 | NOMBRE),
                   family = binomial,
                   data = TES)


tab_model(mod_mixto5)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.04 0.01 – 0.18 <0.001
nume factor [Una] 4.83 2.59 – 9.01 <0.001
isq factor [NO/leve] 3.33 1.81 – 6.14 <0.001
tisu factor [Piel] 2.22 0.88 – 5.63 0.093
infe factor [Leve/Mod] 2.02 0.81 – 5.02 0.130
infe factor [No] 2.21 0.69 – 7.08 0.184
ede factor [Unilateral] 2.64 0.78 – 8.93 0.119
ede factor [Sin/local] 2.93 0.91 – 9.39 0.071
Random Effects
σ2 3.29
τ00 NOMBRE 0.17
ICC 0.05
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.305 / 0.338

MODELO MIXTO 6

mod_mixto6 <- glmer(cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor+ ede_factor+(1 | NOMBRE),
                   family = binomial,
                   data = TES)


tab_model(mod_mixto6)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.04 0.01 – 0.20 <0.001
topo factor
[dor/plan/lat] 		0.73 0.37 – 1.45 0.369
nume factor [Una] 5.41 2.74 – 10.67 <0.001
isq factor [NO/leve] 3.33 1.81 – 6.15 <0.001
tisu factor [Piel] 2.28 0.90 – 5.78 0.084
infe factor [Leve/Mod] 2.02 0.81 – 5.06 0.132
infe factor [No] 2.34 0.72 – 7.64 0.158
ede factor [Unilateral] 2.54 0.75 – 8.58 0.134
ede factor [Sin/local] 3.00 0.94 – 9.62 0.064
Random Effects
σ2 3.29
τ00 NOMBRE 0.17
ICC 0.05
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.307 / 0.341

MODELO MIXTO 7

mod_mixto7 <- glmer(cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor+ ede_factor+ neur_factor + (1 | NOMBRE),
                   family = binomial,
                   data = TES)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00328887 (tol = 0.002, component 1)
tab_model(mod_mixto7)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.15 0.01 – 2.19 0.166
topo factor
[dor/plan/lat] 		0.73 0.37 – 1.46 0.377
nume factor [Una] 5.35 2.70 – 10.61 <0.001
isq factor [NO/leve] 3.40 1.83 – 6.32 <0.001
tisu factor [Piel] 2.07 0.80 – 5.33 0.133
infe factor [Leve/Mod] 2.03 0.80 – 5.11 0.134
infe factor [No] 2.37 0.72 – 7.83 0.157
ede factor [Unilateral] 2.52 0.75 – 8.55 0.137
ede factor [Sin/local] 3.00 0.93 – 9.63 0.065
neur factor
[Leve/mod/charcot] 		0.26 0.03 – 2.39 0.236
Random Effects
σ2 3.29
τ00 NOMBRE 0.21
ICC 0.06
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.318 / 0.359

MODELO MIXTO 8

mod_mixto8 <- glmer(cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor+ ede_factor+ neur_factor + (1 | NOMBRE),
                   family = binomial,
                   data = TES)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00328887 (tol = 0.002, component 1)
tab_model(mod_mixto8)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.15 0.01 – 2.19 0.166
topo factor
[dor/plan/lat] 		0.73 0.37 – 1.46 0.377
nume factor [Una] 5.35 2.70 – 10.61 <0.001
isq factor [NO/leve] 3.40 1.83 – 6.32 <0.001
tisu factor [Piel] 2.07 0.80 – 5.33 0.133
infe factor [Leve/Mod] 2.03 0.80 – 5.11 0.134
infe factor [No] 2.37 0.72 – 7.83 0.157
ede factor [Unilateral] 2.52 0.75 – 8.55 0.137
ede factor [Sin/local] 3.00 0.93 – 9.63 0.065
neur factor
[Leve/mod/charcot] 		0.26 0.03 – 2.39 0.236
Random Effects
σ2 3.29
τ00 NOMBRE 0.21
ICC 0.06
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.318 / 0.359

MODELO MIXTO 9

mod_mixto9 <- glmer(cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor+ ede_factor+ neur_factor + evol + (1 | NOMBRE),
                   family = binomial,
                   data = TES)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00965086 (tol = 0.002, component 1)
tab_model(mod_mixto9)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.16 0.01 – 2.45 0.187
topo factor
[dor/plan/lat] 		0.71 0.35 – 1.46 0.355
nume factor [Una] 5.56 2.75 – 11.22 <0.001
isq factor [NO/leve] 3.52 1.85 – 6.68 <0.001
tisu factor [Piel] 1.83 0.68 – 4.89 0.230
infe factor [Leve/Mod] 2.16 0.83 – 5.62 0.114
infe factor [No] 3.24 0.92 – 11.43 0.068
ede factor [Unilateral] 2.59 0.75 – 8.99 0.134
ede factor [Sin/local] 3.30 1.00 – 10.84 0.049
neur factor
[Leve/mod/charcot] 		0.27 0.03 – 2.59 0.259
evol 0.99 0.99 – 1.00 0.012
Random Effects
σ2 3.29
τ00 NOMBRE 0.42
ICC 0.11
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.368 / 0.439

MODELO MIXTO 10

mod_mixto10 <- glmer(cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor+ ede_factor+ neur_factor + loca_factor+ (1 | NOMBRE),
                   family = binomial,
                   data = TES)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00346577 (tol = 0.002, component 1)
tab_model(mod_mixto10)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.09 0.01 – 1.49 0.093
topo factor
[dor/plan/lat] 		0.81 0.40 – 1.64 0.561
nume factor [Una] 4.78 2.39 – 9.58 <0.001
isq factor [NO/leve] 3.69 1.96 – 6.97 <0.001
tisu factor [Piel] 2.08 0.80 – 5.39 0.132
infe factor [Leve/Mod] 2.23 0.88 – 5.70 0.093
infe factor [No] 2.76 0.82 – 9.30 0.101
ede factor [Unilateral] 2.86 0.83 – 9.91 0.097
ede factor [Sin/local] 3.14 0.97 – 10.16 0.056
neur factor
[Leve/mod/charcot] 		0.25 0.03 – 2.37 0.226
loca factor [Fal] 1.85 1.04 – 3.28 0.035
Random Effects
σ2 3.29
τ00 NOMBRE 0.24
ICC 0.07
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.341 / 0.386

Modelo 11

mod_mixto11 <- glmer(cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor+ ede_factor+ neur_factor + loca_factor+ evol+ (1 | NOMBRE),
                   family = binomial,
                   data = TES)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00625735 (tol = 0.002, component 1)
tab_model(mod_mixto11)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.11 0.01 – 1.78 0.120
topo factor
[dor/plan/lat] 		0.78 0.38 – 1.62 0.506
nume factor [Una] 5.06 2.49 – 10.31 <0.001
isq factor [NO/leve] 3.74 1.95 – 7.17 <0.001
tisu factor [Piel] 1.86 0.70 – 4.99 0.215
infe factor [Leve/Mod] 2.30 0.88 – 5.99 0.090
infe factor [No] 3.52 0.99 – 12.50 0.051
ede factor [Unilateral] 2.88 0.82 – 10.15 0.100
ede factor [Sin/local] 3.41 1.03 – 11.22 0.044
neur factor
[Leve/mod/charcot] 		0.26 0.03 – 2.55 0.247
loca factor [Fal] 1.67 0.93 – 3.00 0.088
evol 0.99 0.99 – 1.00 0.020
Random Effects
σ2 3.29
τ00 NOMBRE 0.42
ICC 0.11
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.379 / 0.450

Modelo 12

mod_mixto12 <- glmer(cicat ~  nume_factor + isq_factor +  infe_factor+ ede_factor+ neur_factor +  evol+ (1 | NOMBRE),
                   family = binomial,
                   data = TES)


tab_model(mod_mixto12)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.16 0.01 – 2.50 0.191
nume factor [Una] 4.93 2.59 – 9.41 <0.001
isq factor [NO/leve] 3.77 2.01 – 7.08 <0.001
infe factor [Leve/Mod] 2.28 0.88 – 5.92 0.090
infe factor [No] 3.90 1.21 – 12.57 0.023
ede factor [Unilateral] 2.64 0.76 – 9.20 0.127
ede factor [Sin/local] 3.31 1.00 – 10.95 0.050
neur factor
[Leve/mod/charcot] 		0.22 0.02 – 2.14 0.194
evol 0.99 0.99 – 1.00 0.010
Random Effects
σ2 3.29
τ00 NOMBRE 0.45
ICC 0.12
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.359 / 0.436

modelo 13

##ELIJO ESTE COMO EFECTOS PRINCIPALES
mod_mixto13 <- glmer(cicat ~ isq_factor + infe_factor + loca_factor +  nume_factor+ ede_factor + evol +  (1 | NOMBRE),
                   family = binomial,
                   data = TES)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00214939 (tol = 0.002, component 1)
tab_model(mod_mixto13)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.02 0.00 – 0.14 <0.001
isq factor [NO/leve] 3.96 2.10 – 7.46 <0.001
infe factor [Leve/Mod] 2.46 0.94 – 6.39 0.066
infe factor [No] 4.47 1.38 – 14.42 0.012
loca factor [Fal] 1.68 0.95 – 2.99 0.076
nume factor [Una] 4.71 2.47 – 8.99 <0.001
ede factor [Unilateral] 2.93 0.83 – 10.29 0.094
ede factor [Sin/local] 3.49 1.05 – 11.52 0.041
evol 0.99 0.99 – 1.00 0.015
Random Effects
σ2 3.29
τ00 NOMBRE 0.39
ICC 0.11
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.358 / 0.427

modelo 14

mod_mixto14 <- glmer(cicat ~ topo_factor + nume_factor + isq_factor + infe_factor+ ede_factor+ neur_factor + loca_factor+ evol+ (1 | NOMBRE),
                   family = binomial,
                   data = TES)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00561179 (tol = 0.002, component 1)
tab_model(mod_mixto14)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.12 0.01 – 2.00 0.139
topo factor
[dor/plan/lat] 		0.79 0.38 – 1.65 0.532
nume factor [Una] 5.06 2.48 – 10.31 <0.001
isq factor [NO/leve] 4.04 2.12 – 7.68 <0.001
infe factor [Leve/Mod] 2.48 0.94 – 6.50 0.066
infe factor [No] 4.66 1.40 – 15.50 0.012
ede factor [Unilateral] 2.80 0.79 – 9.90 0.110
ede factor [Sin/local] 3.50 1.06 – 11.59 0.041
neur factor
[Leve/mod/charcot] 		0.22 0.02 – 2.14 0.192
loca factor [Fal] 1.65 0.92 – 2.97 0.093
evol 0.99 0.99 – 1.00 0.017
Random Effects
σ2 3.29
τ00 NOMBRE 0.48
ICC 0.13
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.370 / 0.451

modelo 15

mod_mixto15 <- glmer(cicat ~ topo_factor + nume_factor + isq_factor + infe_factor+ ede_factor+ neur_factor + loca_factor+ evol+  (1 | NOMBRE),
                   family = binomial,
                   data = TES)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00561179 (tol = 0.002, component 1)
tab_model(mod_mixto15)
  cicat
Predictors Odds Ratios CI p
(Intercept) 0.12 0.01 – 2.00 0.139
topo factor
[dor/plan/lat] 		0.79 0.38 – 1.65 0.532
nume factor [Una] 5.06 2.48 – 10.31 <0.001
isq factor [NO/leve] 4.04 2.12 – 7.68 <0.001
infe factor [Leve/Mod] 2.48 0.94 – 6.50 0.066
infe factor [No] 4.66 1.40 – 15.50 0.012
ede factor [Unilateral] 2.80 0.79 – 9.90 0.110
ede factor [Sin/local] 3.50 1.06 – 11.59 0.041
neur factor
[Leve/mod/charcot] 		0.22 0.02 – 2.14 0.192
loca factor [Fal] 1.65 0.92 – 2.97 0.093
evol 0.99 0.99 – 1.00 0.017
Random Effects
σ2 3.29
τ00 NOMBRE 0.48
ICC 0.13
N NOMBRE 15
Observations 299
Marginal R2 / Conditional R2 0.370 / 0.451

Comparación del rendimiento de los modelos

AIC(mod_mixto,  mod_mixto2, mod_mixto3, mod_mixto4, mod_mixto5, mod_mixto6, mod_mixto7, mod_mixto8, mod_mixto9, mod_mixto10, mod_mixto11, mod_mixto12,mod_mixto13, mod_mixto14)
##             df AIC
## mod_mixto    3 368
## mod_mixto2   4 351
## mod_mixto3   6 350
## mod_mixto4   7 349
## mod_mixto5   9 349
## mod_mixto6  10 350
## mod_mixto7  11 351
## mod_mixto8  11 351
## mod_mixto9  12 344
## mod_mixto10 12 348
## mod_mixto11 13 343
## mod_mixto12 10 342
## mod_mixto13 10 341
## mod_mixto14 12 342
anova(mod_mixto, mod_mixto2, mod_mixto3, mod_mixto4, mod_mixto5, mod_mixto6, mod_mixto7, mod_mixto8, mod_mixto9, mod_mixto10, mod_mixto11, mod_mixto12, mod_mixto13,mod_mixto14, test="Chisq")
## Data: TES
## Models:
## mod_mixto: cicat ~ nume_factor + (1 | NOMBRE)
## mod_mixto2: cicat ~ nume_factor + isq_factor + (1 | NOMBRE)
## mod_mixto3: cicat ~ nume_factor + isq_factor + infe_factor + (1 | NOMBRE)
## mod_mixto4: cicat ~ nume_factor + isq_factor + tisu_factor + infe_factor + (1 | NOMBRE)
## mod_mixto5: cicat ~ nume_factor + isq_factor + tisu_factor + infe_factor + ede_factor + (1 | NOMBRE)
## mod_mixto6: cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor + ede_factor + (1 | NOMBRE)
## mod_mixto12: cicat ~ nume_factor + isq_factor + infe_factor + ede_factor + neur_factor + evol + (1 | NOMBRE)
## mod_mixto13: cicat ~ isq_factor + infe_factor + loca_factor + nume_factor + ede_factor + evol + (1 | NOMBRE)
## mod_mixto7: cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor + ede_factor + neur_factor + (1 | NOMBRE)
## mod_mixto8: cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor + ede_factor + neur_factor + (1 | NOMBRE)
## mod_mixto9: cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor + ede_factor + neur_factor + evol + (1 | NOMBRE)
## mod_mixto10: cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor + ede_factor + neur_factor + loca_factor + (1 | NOMBRE)
## mod_mixto14: cicat ~ topo_factor + nume_factor + isq_factor + infe_factor + ede_factor + neur_factor + loca_factor + evol + (1 | NOMBRE)
## mod_mixto11: cicat ~ topo_factor + nume_factor + isq_factor + tisu_factor + infe_factor + ede_factor + neur_factor + loca_factor + evol + (1 | NOMBRE)
##             npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)    
## mod_mixto      3 368 379   -181      362                        
## mod_mixto2     4 351 366   -172      343 19.13  1   0.000012 ***
## mod_mixto3     6 350 372   -169      338  5.27  2     0.0717 .  
## mod_mixto4     7 349 375   -167      335  3.25  1     0.0713 .  
## mod_mixto5     9 349 383   -166      331  3.44  2     0.1793    
## mod_mixto6    10 350 387   -165      330  0.82  1     0.3647    
## mod_mixto12   10 342 379   -161      322  8.47  0               
## mod_mixto13   10 341 378   -160      321  0.98  0               
## mod_mixto7    11 351 391   -164      329  0.00  1     1.0000    
## mod_mixto8    11 351 391   -164      329  0.00  0               
## mod_mixto9    12 344 388   -160      320  8.92  1     0.0028 ** 
## mod_mixto10   12 348 393   -162      324  0.00  0               
## mod_mixto14   12 342 387   -159      318  5.76  0               
## mod_mixto11   13 343 391   -158      317  1.58  1     0.2088    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

MODELO DE EFECTOS PRINCIPALES mod_mixto13.

ranef(mod_mixto13)$NOMBRE
##                                                                      (Intercept)
## ABUDARA                                                                    0.273
## Centro de rehabilitación, Obesidad y Diabetes. Formosa                     0.493
## CER                                                                        0.512
## Clínica Centro Médico                                                      0.469
## Hosp. Dr.  J.R.Vidal                                                      -0.245
## Hospital Centro de Salud Zenón J. Santillán                               -0.628
## Hospital de Clínicas                                                       0.174
## Hospital de Día de Pie Diabético. Polo Sanitario Malvinas Argentinas      -0.584
## Hospital Güemes                                                           -0.092
## Hospital Perrupato                                                         0.444
## Hospital Posadas                                                           0.274
## Hospital san Martin                                                       -0.434
## Hospital Sommer                                                           -0.567
## HRRG                                                                       0.255
## Provincial Rosario                                                        -0.461
efectos_centros <- ranef(mod_mixto12)$NOMBRE
efectos_centros_df <- data.frame(
  Centro = rownames(efectos_centros),
  Intercepto_aleatorio = efectos_centros[, "(Intercept)"]
)
head(efectos_centros_df)
##                                                   Centro Intercepto_aleatorio
## 1                                                ABUDARA                 0.30
## 2 Centro de rehabilitación, Obesidad y Diabetes. Formosa                 0.51
## 3                                                    CER                 0.60
## 4                                  Clínica Centro Médico                 0.47
## 5                                   Hosp. Dr.  J.R.Vidal                -0.31
## 6            Hospital Centro de Salud Zenón J. Santillán                -0.71
intercepto_fijo <- fixef(mod_mixto13)["(Intercept)"]

efectos_centros_df$Intercepto_total <- intercepto_fijo + efectos_centros_df$Intercepto_aleatorio
library(ggplot2)

ggplot(efectos_centros_df, aes(x = reorder(Centro, Intercepto_total), y = Intercepto_total)) +
  geom_point() +
  coord_flip() +
  labs(x = "Centro", y = "Intercepto total", title = "Variación del intercepto por centro") +
  theme_minimal()

efectos_centros_df$Probabilidad_basal <- plogis(efectos_centros_df$Intercepto_total)
# 1️⃣ Cargar paquete
library(lme4)
library(ggplot2)
library(dplyr)

# 2️⃣ Extraer efectos aleatorios del modelo
efectos_centros <- ranef(mod_mixto13, condVar = TRUE)$NOMBRE

# 3️⃣ Convertir a data.frame
efectos_centros_df <- data.frame(
  Centro = rownames(efectos_centros),
  Intercepto_centro = efectos_centros[ , "(Intercept)"]
)

# 4️⃣ Sumar el intercepto fijo global
intercepto_fijo <- fixef(mod_mixto13)[["(Intercept)"]]
efectos_centros_df <- efectos_centros_df %>%
  mutate(
    Intercepto_total = intercepto_fijo + Intercepto_centro,
    Probabilidad_basal = plogis(Intercepto_total) # convierte log-odds a probabilidad
  )

# 5️⃣ Graficar la probabilidad basal de cicatrización por centro
ggplot(efectos_centros_df, aes(x = reorder(Centro, Probabilidad_basal), y = Probabilidad_basal)) +
  geom_point(size = 3, color = "steelblue") +
  geom_hline(yintercept = mean(efectos_centros_df$Probabilidad_basal), 
             linetype = "dashed", color = "red") +
  coord_flip() +
  labs(
    x = "Centro",
    y = "Probabilidad basal de cicatrización (ajustada)",
    title = "Variación entre centros en la probabilidad basal de cicatrización"
  ) +
  theme_minimal(base_size = 13)

#Residuos vs predichos
library(DHARMa)

sim_res <- simulateResiduals(fittedModel = mod_mixto12, n = 1000)
plot(sim_res)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully

# Test formales
testUniformity(sim_res)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.05, p-value = 0.4
## alternative hypothesis: two-sided
testDispersion(sim_res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1, p-value = 0.5
## alternative hypothesis: two.sided
testZeroInflation(sim_res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1, p-value = 1
## alternative hypothesis: two.sided
library(performance)

# Colinealidad
check_collinearity(mod_mixto13)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##         Term  VIF   VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##   isq_factor 1.08 [1.01, 1.40]     1.04      0.93     [0.72, 0.99]
##  infe_factor 1.22 [1.11, 1.45]     1.11      0.82     [0.69, 0.90]
##  loca_factor 1.07 [1.01, 1.42]     1.03      0.94     [0.70, 0.99]
##  nume_factor 1.09 [1.02, 1.38]     1.04      0.92     [0.73, 0.98]
##   ede_factor 1.24 [1.12, 1.47]     1.11      0.81     [0.68, 0.89]
##         evol 1.05 [1.01, 1.50]     1.03      0.95     [0.67, 0.99]
# Bondad de ajuste (R² marginal y condicional)
r2_nakagawa(mod_mixto13)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.427
##      Marginal R2: 0.358
# Varianza de efectos aleatorios
VarCorr(mod_mixto13)
##  Groups Name        Std.Dev.
##  NOMBRE (Intercept) 0.628

ICC

performance::icc(mod_mixto13)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.107
##   Unadjusted ICC: 0.069

Calibracion

#HyL
library(ResourceSelection)
## Warning: package 'ResourceSelection' was built under R version 4.4.3
## ResourceSelection 0.3-6   2023-06-27
pred <- predict(mod_mixto13, type = "response")
hoslem.test(TES$cicat, pred, g = 10)
## Warning in Ops.factor(1, y): '-' no es significativo para factores
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  TES$cicat, pred
## X-squared = NA, df = 8, p-value = NA
#Brier
library(DescTools)
## 
## Adjuntando el paquete: 'DescTools'
## The following objects are masked from 'package:psych':
## 
##     AUC, ICC, SD
## The following object is masked from 'package:data.table':
## 
##     %like%
pred <- predict(mod_mixto13, type = "response")
BrierScore <- BrierScore(TES$cicat, pred)
## Warning in Ops.factor(resp, (1 - pred)^2): '*' no es significativo para
## factores
## Warning in Ops.factor(1, resp): '-' no es significativo para factores
BrierScore
## [1] NA

Calibración:NO SE QUE ES ESTO EN MODELO MIXTO

TES$lp <- predict(mod_mixto13, type = "link")  # predicción en escala logit

mod_slope <- glm(cicat ~ lp, family = binomial, data = TES)
summary(mod_slope)
## 
## Call:
## glm(formula = cicat ~ lp, family = binomial, data = TES)
## 
## Coefficients:
##             Estimate Std. Error z value            Pr(>|z|)    
## (Intercept)  -0.0376     0.1495   -0.25                 0.8    
## lp            1.0881     0.1340    8.12 0.00000000000000047 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 405.76  on 298  degrees of freedom
## Residual deviance: 303.16  on 297  degrees of freedom
## AIC: 307.2
## 
## Number of Fisher Scoring iterations: 4

auc

library(pROC)
## Warning: package 'pROC' was built under R version 4.4.3
## Type 'citation("pROC")' for a citation.
## 
## Adjuntando el paquete: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
pred <- predict(mod_mixto13, type = "response")
roc_obj <- roc(TES$cicat, pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
auc(roc_obj)
## Area under the curve: 0.82
plot(roc_obj, print.auc = TRUE, col = "blue", main = "Curva ROC - Modelo mixto")