Untitled
2025-10-04
library(irr)## Cargando paquete requerido: lpSolvelibrary(expss)## Cargando paquete requerido: maditr##
## To modify variables or add new variables:
## let(mtcars, new_var = 42, new_var2 = new_var*hp) %>% head()library(Hmisc)## Warning: package 'Hmisc' was built under R version 4.4.3## Registered S3 methods overwritten by 'Hmisc':
## method from
## [.labelled expss
## print.labelled expss
## as.data.frame.labelled expss##
## Adjuntando el paquete: 'Hmisc'## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(PerformanceAnalytics)## Warning: package 'PerformanceAnalytics' was built under R version 4.4.3## Cargando paquete requerido: xts## Warning: package 'xts' was built under R version 4.4.3## Cargando paquete requerido: zoo##
## Adjuntando el paquete: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Adjuntando el paquete: 'xts'## The following objects are masked from 'package:maditr':
##
## first, last##
## Adjuntando el paquete: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legendlibrary(readxl)
#sad <- read_excel("C:/Users/Administrador/Downloads/sad.xlsx",
# sheet = "R cuestionarios")
#View(sad)
#carga de datos
library(readxl)
sad <- read_excel("C:/Users/Administrador/Downloads/TESIC.xlsx")## New names:
## • `` -> `...44`View(sad)attach(sad)
#ahora sumemos los puntos de las preguntas
#los valores 5 a 15 son las columnas correspondientes a las preguntas 1 a 10
View(sad[ , 25:34])
sad$puntaje1 <- rowSums(sad[ , 25:34])#Buscamos proporción de datos faltantes (Missingvalues)---- NO SE DETECTAN DATOS FALTANTES
#is.na(sad[ , 25:34])
#which(is.na(sad[ , 25:34]))#pasar variables categóricas a factor
names(sad)<- tolower(names (sad))
datos1=sad[ , 25:34]
sad <- sad %>%
mutate(across(19:35, as.factor))
#crear cicat_factor
sad$cicat_factor <- factor(sad$cicat, levels = c(0, 1), labels = c( "No Cicatrizo", "Cicatrizo"))
table(sad$cicat_factor)##
## No Cicatrizo Cicatrizo
## 174 232#paso variables chr a factor
sad <- sad %>% mutate_if(is.character, as.factor)
#crear ampu_factor
sad$ampu_factor <- factor(sad$ampu, levels = c(0, 1), labels = c( "No Amputado", "Amputado"))
table(sad$ampu_factor)##
## No Amputado Amputado
## 348 58#crear fallecio_factor
sad$fallecio_factor <- factor(sad$fallecio, levels = c(0, 1), labels = c( "No Fallecido", "Fallecido"))
table(sad$fallecio_factor)##
## No Fallecido Fallecido
## 392 14#variable unificada
sad <- sad %>%
mutate(
estado = tolower(trimws(estado))
)
table(sad$estado)##
## ampu cerro fallecio persiste
## 58 232 14 102sad <- sad %>%
mutate(
estado_unif = case_when(
estado %in% c("ampu") ~ "amputacion",
estado %in% c("cerro") ~ "cerro_herida",
estado %in% c("murio", "fallecio") ~ "fallecio",
estado %in% c("persiste") ~ "persiste",
estado %in% c("perdido") ~ "perdido_seguimiento",
TRUE ~ NA_character_
)
)
table(sad$estado_unif)##
## amputacion cerro_herida fallecio persiste
## 58 232 14 102sad <- sad %>%
mutate(
evento_ampu = case_when(
estado_unif %in% c("cerro_herida","fallecio") ~ 2,
estado_unif %in% c("amputacion") ~ 1,
estado_unif %in% c("perdido_seguimiento", "persiste") ~ 0,
TRUE ~ NA_real_
)
)
table(sad$evento_ampu)##
## 0 1 2
## 102 58 246tiempo_evento_am <- Surv(sad$dias_out,sad$ampu)#Distribución de frecuencia de las respuestas----
fre(c(sad[ , 25:34]))| Count | Valid percent | Percent | Responses, % | Cumulative responses, % | |
|---|---|---|---|---|---|
| 1 | 196 | 48.3 | 48.3 | 48.3 | 48.3 |
| 2 | 134 | 33.0 | 33.0 | 33.0 | 81.3 |
| 3 | 76 | 18.7 | 18.7 | 18.7 | 100.0 |
| #Total | 406 | 100 | 100 | 100 | |
| <NA> | 0 | 0.0 | |||
| 1 | 152 | 37.4 | 37.4 | 37.4 | 37.4 |
| 2 | 104 | 25.6 | 25.6 | 25.6 | 63.1 |
| 3 | 150 | 36.9 | 36.9 | 36.9 | 100.0 |
| #Total | 406 | 100 | 100 | 100 | |
| <NA> | 0 | 0.0 | |||
| 1 | 242 | 59.6 | 59.6 | 59.6 | 59.6 |
| 2 | 90 | 22.2 | 22.2 | 22.2 | 81.8 |
| 3 | 74 | 18.2 | 18.2 | 18.2 | 100.0 |
| #Total | 406 | 100 | 100 | 100 | |
| <NA> | 0 | 0.0 | |||
| 0 | 206 | 50.7 | 50.7 | 50.7 | 50.7 |
| 1 | 31 | 7.6 | 7.6 | 7.6 | 58.4 |
| 2 | 47 | 11.6 | 11.6 | 11.6 | 70.0 |
| 3 | 122 | 30.0 | 30.0 | 30.0 | 100.0 |
| #Total | 406 | 100 | 100 | 100 | |
| <NA> | 0 | 0.0 | |||
| 0 | 86 | 21.2 | 21.2 | 21.2 | 21.2 |
| 1 | 32 | 7.9 | 7.9 | 7.9 | 29.1 |
| 2 | 223 | 54.9 | 54.9 | 54.9 | 84.0 |
| 3 | 65 | 16.0 | 16.0 | 16.0 | 100.0 |
| #Total | 406 | 100 | 100 | 100 | |
| <NA> | 0 | 0.0 | |||
| 0 | 143 | 35.2 | 35.2 | 35.2 | 35.2 |
| 1 | 98 | 24.1 | 24.1 | 24.1 | 59.4 |
| 2 | 124 | 30.5 | 30.5 | 30.5 | 89.9 |
| 3 | 41 | 10.1 | 10.1 | 10.1 | 100.0 |
| #Total | 406 | 100 | 100 | 100 | |
| <NA> | 0 | 0.0 | |||
| 0 | 4 | 1.0 | 1.0 | 1.0 | 1.0 |
| 1 | 25 | 6.2 | 6.2 | 6.2 | 7.1 |
| 2 | 341 | 84.0 | 84.0 | 84.0 | 91.1 |
| 3 | 36 | 8.9 | 8.9 | 8.9 | 100.0 |
| #Total | 406 | 100 | 100 | 100 | |
| <NA> | 0 | 0.0 | |||
| 1 | 32 | 7.9 | 7.9 | 7.9 | 7.9 |
| 2 | 133 | 32.8 | 32.8 | 32.8 | 40.6 |
| 3 | 241 | 59.4 | 59.4 | 59.4 | 100.0 |
| #Total | 406 | 100 | 100 | 100 | |
| <NA> | 0 | 0.0 | |||
| 1 | 188 | 46.3 | 46.3 | 46.3 | 46.3 |
| 2 | 157 | 38.7 | 38.7 | 38.7 | 85.0 |
| 3 | 61 | 15.0 | 15.0 | 15.0 | 100.0 |
| #Total | 406 | 100 | 100 | 100 | |
| <NA> | 0 | 0.0 | |||
| 1 | 19 | 4.7 | 4.7 | 4.7 | 4.7 |
| 2 | 51 | 12.6 | 12.6 | 12.6 | 17.2 |
| 3 | 336 | 82.8 | 82.8 | 82.8 | 100.0 |
| #Total | 406 | 100 | 100 | 100 | |
| <NA> | 0 | 0.0 |
hist.data.frame(sad[ , 25:34])
#Confiabilidad I: Analisis de consistencia interna#horn. Ojo! sugiere 4 factores pero se vuelve inestable
datos<-sad[ , 25:34]
fa.parallel(datos, fa = "fa", cor = "poly")## Converted non-numeric input to numeric
## Converted non-numeric input to numeric## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric
## Converted non-numeric input to numeric## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully## Converted non-numeric input to numeric
## Parallel analysis suggests that the number of factors = 4 and the number of components = NAcolnames(datos)<- c("LOCAL","TOPO","ZON","ISQ","INF","EDE","NEU","PROF","AREA","FASE")
library(corrplot)## Warning: package 'corrplot' was built under R version 4.4.3## corrplot 0.95 loadedvars <-datos[, c("LOCAL","TOPO","ZON","ISQ","INF","EDE","NEU","PROF","AREA","FASE")]
poly <- polychoric(vars)## Converted non-numeric input to numericR_poly <- poly$rho
colnames(R_poly) <- colnames(vars)
rownames(R_poly) <- colnames(vars)
#bartlett
cortest.bartlett(R_poly, n = nrow(vars))## $chisq
## [1] 2300
##
## $p.value
## [1] 0
##
## $df
## [1] 45#KMO
KMO(R_poly)## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = R_poly)
## Overall MSA = 0.72
## MSA for each item =
## LOCAL TOPO ZON ISQ INF EDE NEU PROF AREA FASE
## 0.34 0.76 0.75 0.70 0.84 0.73 0.46 0.68 0.83 0.69kmo <- KMO(R_poly)
kmo_global <- kmo$MSA
kmo_vars <- kmo$MSAi
tabla_kmo <- data.frame(
Variable = names(kmo_vars),
KMO = round(kmo_vars, 2)
)
tabla_kmo## Variable KMO
## LOCAL LOCAL 0.34
## TOPO TOPO 0.76
## ZON ZON 0.75
## ISQ ISQ 0.70
## INF INF 0.84
## EDE EDE 0.73
## NEU NEU 0.46
## PROF PROF 0.68
## AREA AREA 0.83
## FASE FASE 0.69kmo_global## [1] 0.723#gráfico
corrplot(
R_poly,
method = "ellipse",
type = "lower",
addCoef.col = "black",
number.digits = 2, # 👈 redondeo
tl.col = "black",
tl.srt = 45, # 👈 rotación de etiquetas
diag = FALSE
)
fa(R_poly, nfactors = 3, fm = "minres", rotate = "promax")## Loading required namespace: GPArotation## Factor Analysis using method = minres
## Call: fa(r = R_poly, nfactors = 3, rotate = "promax", fm = "minres")
## Standardized loadings (pattern matrix) based upon correlation matrix
## MR1 MR3 MR2 h2 u2 com
## LOCAL -0.29 0.29 0.09 0.13 0.870 2.2
## TOPO 0.18 0.65 -0.26 0.58 0.422 1.5
## ZON -0.02 0.94 -0.11 0.86 0.140 1.0
## ISQ 0.05 0.18 -0.51 0.30 0.704 1.3
## INF 0.78 0.18 0.18 0.74 0.265 1.2
## EDE 0.20 0.15 0.41 0.24 0.760 1.8
## NEU 0.00 -0.07 0.81 0.65 0.351 1.0
## PROF 0.97 -0.01 0.07 0.92 0.076 1.0
## AREA 0.20 0.80 0.07 0.79 0.212 1.1
## FASE 0.81 0.19 -0.13 0.83 0.169 1.2
##
## MR1 MR3 MR2
## SS loadings 2.55 2.29 1.19
## Proportion Var 0.25 0.23 0.12
## Cumulative Var 0.25 0.48 0.60
## Proportion Explained 0.42 0.38 0.20
## Cumulative Proportion 0.42 0.80 1.00
##
## With factor correlations of
## MR1 MR3 MR2
## MR1 1.00 0.34 -0.12
## MR3 0.34 1.00 0.09
## MR2 -0.12 0.09 1.00
##
## Mean item complexity = 1.3
## Test of the hypothesis that 3 factors are sufficient.
##
## df null model = 45 with the objective function = 5.74
## df of the model are 18 and the objective function was 0.64
##
## The root mean square of the residuals (RMSR) is 0.04
## The df corrected root mean square of the residuals is 0.07
##
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy
## MR1 MR3 MR2
## Correlation of (regression) scores with factors 0.98 0.96 0.87
## Multiple R square of scores with factors 0.96 0.92 0.75
## Minimum correlation of possible factor scores 0.91 0.84 0.50fa(R_poly, nfactors = 4, fm = "minres", rotate = "promax")## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully## Factor Analysis using method = minres
## Call: fa(r = R_poly, nfactors = 4, rotate = "promax", fm = "minres")
## Standardized loadings (pattern matrix) based upon correlation matrix
## MR1 MR3 MR2 MR4 h2 u2 com
## LOCAL -0.03 0.09 0.05 0.66 0.45 0.54938 1.1
## TOPO -0.09 0.87 -0.19 -0.19 0.75 0.25447 1.2
## ZON -0.08 0.93 -0.03 0.16 0.81 0.19245 1.1
## ISQ 0.07 0.17 -0.50 -0.02 0.30 0.70019 1.3
## INF 0.83 0.04 0.13 0.02 0.74 0.25856 1.1
## EDE 0.04 0.26 0.47 -0.14 0.32 0.67694 1.8
## NEU 0.02 -0.10 0.74 0.13 0.59 0.41225 1.1
## PROF 0.86 -0.03 0.03 -0.26 0.88 0.12300 1.2
## AREA 0.23 0.70 0.08 0.27 0.79 0.21106 1.5
## FASE 1.05 -0.08 -0.24 0.13 1.00 -0.00061 1.1
##
## MR1 MR3 MR2 MR4
## SS loadings 2.57 2.20 1.15 0.69
## Proportion Var 0.26 0.22 0.12 0.07
## Cumulative Var 0.26 0.48 0.59 0.66
## Proportion Explained 0.39 0.33 0.17 0.10
## Cumulative Proportion 0.39 0.72 0.90 1.00
##
## With factor correlations of
## MR1 MR3 MR2 MR4
## MR1 1.00 0.60 0.03 -0.21
## MR3 0.60 1.00 0.00 -0.01
## MR2 0.03 0.00 1.00 0.04
## MR4 -0.21 -0.01 0.04 1.00
##
## Mean item complexity = 1.2
## Test of the hypothesis that 4 factors are sufficient.
##
## df null model = 45 with the objective function = 5.74
## df of the model are 11 and the objective function was 0.12
##
## The root mean square of the residuals (RMSR) is 0.01
## The df corrected root mean square of the residuals is 0.03
##
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy
## MR1 MR3 MR2 MR4
## Correlation of (regression) scores with factors 0.99 0.95 0.87 0.85
## Multiple R square of scores with factors 0.98 0.91 0.76 0.72
## Minimum correlation of possible factor scores 0.97 0.82 0.52 0.430.7 → muy bien
0.4–0.7 → aceptable < 0.3 → problemático
🔀 com = complejidad
👉 cuántos factores explican ese ítem
📌 ideal ≈ 1
Ejemplo:
ZON → com = 1.03 👉 carga en UN solo factor → perfecto ✔️ LOCAL → com = 2.19 👉 carga en varios factores → problema ⚠️ Variables fuertes (muy buenas) PROF (h2 0.92) FASE (0.83) ZON (0.86) AREA (0.79) INF (0.73)
👉 sostienen bien la estructura
⚠️ Variables intermedias NEU (0.65) TOPO (0.57)
👉 aceptables
❌ Variables problemáticas LOCAL (h2 0.13, com 2.19) EDE (h2 0.24) ISQ (h2 0.29)
h2 → qué tan bien el modelo explica cada variable u2 → lo que queda sin explicar com → si el ítem pertenece a un solo factor o a varios Un factor representa:
una dimensión latente que explica por qué ciertas variables se correlacionan entre sí.
👉 Entonces:
las variables de un mismo factor tienden a variar juntas pero eso no siempre significa que “aumentan todas en el mismo sentido” Las variables de un mismo factor están asociadas entre sí, ya sea de manera directa o inversa, porque reflejan una misma dimensión subyacente Factor isquémico (tu MR2) NEU ↑ ISQ ↓ (carga negativa)
👉 esto te dice algo interesante:
💥 ese factor no es “puro” 👉 es una dimensión donde:
neuropatía y isquemia no van necesariamente juntas
👉 esto es clínicamente coherente
#pruebo 2 factores
Fa2<-fa(datos1, nfactors = 2,rotate = "promax" , fm ="minres" , cor = "poly" , correct = 0)
fa.diagram(Fa2)
fa.plot(Fa2, cut = 0.3, labels = names(datos1))
#pruebo 3 factores
Fa3<-fa(datos1, nfactors = 3,rotate = "promax" , fm ="minres" , cor = "poly" , correct = 0)
fa.diagram(Fa3)
fa.plot(Fa3, cut = 0.3, labels = names(datos1))
#Rpoly por dominio con 3 factores
#Dominio 1
varsd1 <-datos[, c("INF","PROF","FASE")]
polyd1 <- polychoric(varsd1)## Converted non-numeric input to numericR_polyd1 <- polyd1$rho
colnames(R_polyd1) <- colnames(varsd1)
rownames(R_polyd1) <- colnames(varsd1)
#bartlett
cortest.bartlett(R_polyd1, n = nrow(varsd1))## $chisq
## [1] 927
##
## $p.value
## [1] 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000108
##
## $df
## [1] 3#gráfico
corrplot(
R_polyd1,
method = "ellipse",
type = "lower",
addCoef.col = "black",
number.digits = 2, # 👈 redondeo
tl.col = "black",
tl.srt = 45, # 👈 rotación de etiquetas
diag = FALSE
)
#Dominio 2
varsd2 <-datos[, c("ZON","AREA","TOPO")]
polyd2 <- polychoric(varsd2)## Converted non-numeric input to numericR_polyd2 <- polyd2$rho
colnames(R_polyd2) <- colnames(varsd2)
rownames(R_polyd2) <- colnames(varsd2)
#bartlett
cortest.bartlett(R_polyd2, n = nrow(varsd2))## $chisq
## [1] 637
##
## $p.value
## [1] 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000104
##
## $df
## [1] 3#gráfico
corrplot(
R_polyd2,
method = "ellipse",
type = "lower",
addCoef.col = "black",
number.digits = 2, # 👈 redondeo
tl.col = "black",
tl.srt = 45, # 👈 rotación de etiquetas
diag = FALSE
)
#DOMINIO 3
varsd3 <-datos[, c("ISQ","EDE","NEU")]
polyd3 <- polychoric(varsd3)## Converted non-numeric input to numericR_polyd3 <- polyd3$rho
colnames(R_polyd3) <- colnames(varsd3)
rownames(R_polyd3) <- colnames(varsd3)
#bartlett
cortest.bartlett(R_polyd3, n = nrow(varsd3))## $chisq
## [1] 114
##
## $p.value
## [1] 0.00000000000000000000000188
##
## $df
## [1] 3#gráfico
corrplot(
R_polyd3,
method = "ellipse",
type = "lower",
addCoef.col = "black",
number.digits = 2, # 👈 redondeo
tl.col = "black",
tl.srt = 45, # 👈 rotación de etiquetas
diag = FALSE
)
#pruebo 4 factores
Fa4<-fa(datos1, nfactors = 4,rotate = "promax" , fm ="minres" , cor = "poly" , correct = 0)## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefullyfa.diagram(Fa4)
fa.plot(Fa4, cut = 0.3, labels = names(datos1))
library(tableone)
#Las variables van directamente sin el nombre del dataframe
catvars = c( "sexo","tbq", "iam", "hta", "dial", "local","topo", "inf", "ede", "neu", "prof", "area", "zon", "fase")
#continuas
vars = c("dias_ev","san_elian", "sexo","tbq", "iam", "hta", "dial", "local","topo", "inf", "ede", "neu", "prof", "area", "zon", "fase" )
tabla_1 <- CreateTableOne(vars = vars, factorVars = catvars, data = sad)
print(tabla_1)##
## Overall
## n 406
## dias_ev (mean (SD)) 50.26 (112.85)
## san_elian (mean (SD)) 18.26 (4.09)
## sexo = M (%) 324 (79.8)
## tbq (%)
## 0 169 (41.8)
## 1 80 (19.8)
## 2 155 (38.4)
## iam = 1 (%) 83 (20.5)
## hta = 1 (%) 273 (67.4)
## dial = 1 (%) 19 ( 4.7)
## local (%)
## 1 196 (48.3)
## 2 134 (33.0)
## 3 76 (18.7)
## topo (%)
## 1 152 (37.4)
## 2 104 (25.6)
## 3 150 (36.9)
## inf (%)
## 0 86 (21.2)
## 1 32 ( 7.9)
## 2 223 (54.9)
## 3 65 (16.0)
## ede (%)
## 0 143 (35.2)
## 1 98 (24.1)
## 2 124 (30.5)
## 3 41 (10.1)
## neu (%)
## 0 4 ( 1.0)
## 1 25 ( 6.2)
## 2 341 (84.0)
## 3 36 ( 8.9)
## prof (%)
## 1 32 ( 7.9)
## 2 133 (32.8)
## 3 241 (59.4)
## area (%)
## 1 188 (46.3)
## 2 157 (38.7)
## 3 61 (15.0)
## zon (%)
## 1 242 (59.6)
## 2 90 (22.2)
## 3 74 (18.2)
## fase (%)
## 1 19 ( 4.7)
## 2 51 (12.6)
## 3 336 (82.8)kableone(tabla_1)| Overall | |
|---|---|
| n | 406 |
| dias_ev (mean (SD)) | 50.26 (112.85) |
| san_elian (mean (SD)) | 18.26 (4.09) |
| sexo = M (%) | 324 (79.8) |
| tbq (%) | |
| 0 | 169 (41.8) |
| 1 | 80 (19.8) |
| 2 | 155 (38.4) |
| iam = 1 (%) | 83 (20.5) |
| hta = 1 (%) | 273 (67.4) |
| dial = 1 (%) | 19 ( 4.7) |
| local (%) | |
| 1 | 196 (48.3) |
| 2 | 134 (33.0) |
| 3 | 76 (18.7) |
| topo (%) | |
| 1 | 152 (37.4) |
| 2 | 104 (25.6) |
| 3 | 150 (36.9) |
| inf (%) | |
| 0 | 86 (21.2) |
| 1 | 32 ( 7.9) |
| 2 | 223 (54.9) |
| 3 | 65 (16.0) |
| ede (%) | |
| 0 | 143 (35.2) |
| 1 | 98 (24.1) |
| 2 | 124 (30.5) |
| 3 | 41 (10.1) |
| neu (%) | |
| 0 | 4 ( 1.0) |
| 1 | 25 ( 6.2) |
| 2 | 341 (84.0) |
| 3 | 36 ( 8.9) |
| prof (%) | |
| 1 | 32 ( 7.9) |
| 2 | 133 (32.8) |
| 3 | 241 (59.4) |
| area (%) | |
| 1 | 188 (46.3) |
| 2 | 157 (38.7) |
| 3 | 61 (15.0) |
| zon (%) | |
| 1 | 242 (59.6) |
| 2 | 90 (22.2) |
| 3 | 74 (18.2) |
| fase (%) | |
| 1 | 19 ( 4.7) |
| 2 | 51 (12.6) |
| 3 | 336 (82.8) |
#variables categorizadas
#Las variables van directamente sin el nombre del dataframe#ROC
#install.packages("timeROC")
library(timeROC)## Warning: package 'timeROC' was built under R version 4.4.3table(sad$evento_ampu, sad$dias_out <= 90)##
## FALSE TRUE
## 0 60 42
## 1 8 50
## 2 113 133roc_90 <- timeROC(
T = sad$dias_out,
delta = sad$evento_ampu,
marker = sad$san_elian, # 👈 ACÁ está San Elián
cause = 1,
times = 90,
iid = TRUE
)
table(sad$evento_ampu)##
## 0 1 2
## 102 58 246roc_90$AUC## NULLplot(roc_90, time = 90)
roc_120 <- timeROC(
T = sad$dias_out,
delta = sad$evento_ampu,
marker = sad$san_elian,
cause = 1,
times = 120,
iid = TRUE
)
roc_120$AUC## NULLplot(roc_120, time = 120)
roc_180 <- timeROC(
T = sad$dias_out,
delta = sad$evento_ampu,
marker = sad$san_elian,
cause = 1,
times = 180,
iid = TRUE
)
roc_180$AUC## NULLplot(roc_180, time = 180)
sens <- roc_90$TP[,1] # sensibilidad
spec <- 1 - roc_90$FP[,1] # especificidad
cutoffs <- roc_90$cutoffs
youden <- sens + spec - 1
best <- which.max(youden)
#La capacidad discriminativa del puntaje de San Elián para predecir amputación mayor a 90 días se evaluó mediante curvas ROC dependientes del tiempo, considerando la presencia de eventos competitivos. El área bajo la curva fue de 0.829.#consistencia interna omega de mcdonald
vars <-datos[, c("LOCAL","TOPO","ZON","ISQ","INF","EDE","NEU","PROF","AREA","FASE")]
str(vars)## tibble [406 × 10] (S3: tbl_df/tbl/data.frame)
## $ LOCAL: Factor w/ 3 levels "1","2","3": 2 1 1 1 2 1 1 2 2 2 ...
## $ TOPO : Factor w/ 3 levels "1","2","3": 1 3 1 3 2 3 1 1 1 1 ...
## $ ZON : Factor w/ 3 levels "1","2","3": 1 1 1 3 1 3 1 1 1 1 ...
## $ ISQ : Factor w/ 4 levels "0","1","2","3": 3 4 3 1 3 4 4 2 1 1 ...
## $ INF : Factor w/ 4 levels "0","1","2","3": 1 3 3 3 3 3 3 1 1 1 ...
## $ EDE : Factor w/ 4 levels "0","1","2","3": 3 1 2 3 1 2 3 1 1 1 ...
## $ NEU : Factor w/ 4 levels "0","1","2","3": 3 3 3 3 3 3 3 3 3 3 ...
## $ PROF : Factor w/ 3 levels "1","2","3": 2 3 3 3 3 3 3 2 2 2 ...
## $ AREA : Factor w/ 3 levels "1","2","3": 1 2 2 2 2 2 1 2 1 1 ...
## $ FASE : Factor w/ 3 levels "1","2","3": 2 3 3 3 3 3 3 3 2 2 ...omega_total <- omega(
vars,
nfactors = 3,
poly = TRUE,
labels = colnames(vars)
)## Converted non-numeric input to numeric
omega_total## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.77
## G.6: 0.87
## Omega Hierarchical: 0.55
## Omega H asymptotic: 0.63
## Omega Total 0.88
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* F2* F3* h2 h2 u2 p2 com
## V1- 0.22 -0.25 0.13 0.87 0.01 2.55
## V2 0.54 0.47 -0.22 0.58 0.58 0.42 0.51 2.39
## V3 0.58 0.72 0.86 0.86 0.14 0.39 1.92
## V4 -0.50 0.30 0.30 0.70 0.12 1.38
## V5 0.64 0.55 0.74 0.74 0.26 0.55 2.13
## V6 0.21 0.40 0.24 0.24 0.76 0.18 1.98
## V7- -0.80 0.65 0.65 0.35 0.02 1.03
## V8 0.66 0.70 0.92 0.92 0.08 0.47 2.01
## V9 0.63 0.61 0.79 0.79 0.21 0.50 2.12
## V10 0.69 0.57 0.83 0.83 0.17 0.57 2.10
##
## With Sums of squares of:
## g F1* F2* F3* h2
## 2.4 1.2 1.2 1.2 4.4
##
## general/max 0.55 max/min = 3.7
## mean percent general = 0.33 with sd = 0.22 and cv of 0.68
## Explained Common Variance of the general factor = 0.4
##
## The degrees of freedom are 18 and the fit is 0.64
## The number of observations was 406 with Chi Square = 255 with prob < 0.000000000000000000000000000000000000000000072
## The root mean square of the residuals is 0.04
## The df corrected root mean square of the residuals is 0.07
## RMSEA index = 0.18 and the 90 % confidence intervals are 0.161 0.2
## BIC = 147
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 35 and the fit is 3.07
## The number of observations was 406 with Chi Square = 1230 with prob < 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000094
## The root mean square of the residuals is 0.18
## The df corrected root mean square of the residuals is 0.21
##
## RMSEA index = 0.29 and the 90 % confidence intervals are 0.277 0.304
## BIC = 1020
##
## Measures of factor score adequacy
## g F1* F2* F3*
## Correlation of scores with factors 0.79 0.76 0.81 0.86
## Multiple R square of scores with factors 0.62 0.58 0.66 0.75
## Minimum correlation of factor score estimates 0.25 0.16 0.32 0.49
##
## Total, General and Subset omega for each subset
## g F1* F2* F3*
## Omega total for total scores and subscales 0.88 0.93 0.74 0.37
## Omega general for total scores and subscales 0.55 0.51 0.42 0.08
## Omega group for total scores and subscales 0.20 0.42 0.32 0.28#pruebo para 4 factores
vars <-datos[, c("LOCAL","TOPO","ZON","ISQ","INF","EDE","NEU","PROF","AREA","FASE")]
str(vars)## tibble [406 × 10] (S3: tbl_df/tbl/data.frame)
## $ LOCAL: Factor w/ 3 levels "1","2","3": 2 1 1 1 2 1 1 2 2 2 ...
## $ TOPO : Factor w/ 3 levels "1","2","3": 1 3 1 3 2 3 1 1 1 1 ...
## $ ZON : Factor w/ 3 levels "1","2","3": 1 1 1 3 1 3 1 1 1 1 ...
## $ ISQ : Factor w/ 4 levels "0","1","2","3": 3 4 3 1 3 4 4 2 1 1 ...
## $ INF : Factor w/ 4 levels "0","1","2","3": 1 3 3 3 3 3 3 1 1 1 ...
## $ EDE : Factor w/ 4 levels "0","1","2","3": 3 1 2 3 1 2 3 1 1 1 ...
## $ NEU : Factor w/ 4 levels "0","1","2","3": 3 3 3 3 3 3 3 3 3 3 ...
## $ PROF : Factor w/ 3 levels "1","2","3": 2 3 3 3 3 3 3 2 2 2 ...
## $ AREA : Factor w/ 3 levels "1","2","3": 1 2 2 2 2 2 1 2 1 1 ...
## $ FASE : Factor w/ 3 levels "1","2","3": 2 3 3 3 3 3 3 3 2 2 ...omega_total2 <- omega(
vars,
nfactors = 4,
poly = TRUE,
labels = colnames(vars)
)## Converted non-numeric input to numeric## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully## Warning in cov2cor(t(w) %*% r %*% w): diag(V) had non-positive or NA entries;
## the non-finite result may be dubious
omega_total2## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.77
## G.6: 0.87
## Omega Hierarchical: 0.66
## Omega H asymptotic: 0.73
## Omega Total 0.9
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* F2* F3* F4* h2 h2 u2 p2 com
## V1- -0.63 0.45 0.45 0.55 0.02 1.11
## V2 0.42 0.73 -0.22 0.75 0.75 0.25 0.24 1.94
## V3 0.40 0.79 0.81 0.81 0.19 0.20 1.53
## V4 -0.49 0.30 0.30 0.70 0.12 1.47
## V5 0.83 0.74 0.74 0.26 0.94 1.12
## V6 0.22 0.49 0.32 0.32 0.68 0.10 2.01
## V7- -0.75 0.59 0.59 0.41 0.02 1.08
## V8 0.88 -0.26 0.88 0.88 0.12 0.89 1.20
## V9 0.56 0.62 0.21 0.79 0.79 0.21 0.39 2.34
## V10 0.99 1.00 1.00 0.00 0.98 1.09
##
## With Sums of squares of:
## g F1* F2* F3* F4* h2
## 3.19 0.00 1.62 1.16 0.63 4.89
##
## general/max 0.65 max/min = Inf
## mean percent general = 0.39 with sd = 0.39 and cv of 1.01
## Explained Common Variance of the general factor = 0.48
##
## The degrees of freedom are 11 and the fit is 0.12
## The number of observations was 406 with Chi Square = 46.4 with prob < 0.0000027
## The root mean square of the residuals is 0.01
## The df corrected root mean square of the residuals is 0.03
## RMSEA index = 0.089 and the 90 % confidence intervals are 0.064 0.116
## BIC = -19.6
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 35 and the fit is 2.63
## The number of observations was 406 with Chi Square = 1050 with prob < 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000072
## The root mean square of the residuals is 0.18
## The df corrected root mean square of the residuals is 0.2
##
## RMSEA index = 0.267 and the 90 % confidence intervals are 0.254 0.282
## BIC = 840
##
## Measures of factor score adequacy
## g F1* F2* F3* F4*
## Correlation of scores with factors 1 0 0.95 0.87 0.86
## Multiple R square of scores with factors 1 0 0.90 0.77 0.75
## Minimum correlation of factor score estimates 1 -1 0.80 0.53 0.49
##
## Total, General and Subset omega for each subset
## g F1* F2* F3* F4*
## Omega total for total scores and subscales 0.90 NA 0.90 0.70 0.69
## Omega general for total scores and subscales 0.66 NA 0.26 0.64 0.37
## Omega group for total scores and subscales 0.18 NA 0.63 0.07 0.32#en pacientes con isquemia
table(sad$isq)##
## 0 1 2 3
## 206 31 47 122sad_isq <- subset(sad, isq %in% c("2","3"))
roc_90i <- timeROC(
T = sad_isq$dias_out,
delta = sad_isq$evento_ampu,
marker = sad_isq$san_elian, # 👈 ACÁ está San Elián
cause = 1,
times = 90,
iid = TRUE
)
table(sad_isq$evento_ampu)##
## 0 1 2
## 49 46 74roc_90i$AUC## NULLplot(roc_90i, time = 90)
roc_120i <- timeROC(
T = sad_isq$dias_out,
delta = sad_isq$evento_ampu,
marker = sad_isq$san_elian,
cause = 1,
times = 120,
iid = TRUE
)
roc_120i$AUC## NULLplot(roc_120i, time = 120)
roc_180i <- timeROC(
T = sad_isq$dias_out,
delta = sad_isq$evento_ampu,
marker = sad_isq$san_elian,
cause = 1,
times = 180,
iid = TRUE
)
roc_180i$AUC## NULLplot(roc_180i, time = 180)
sad_inf <- subset(sad, inf %in% c("2","3"))
roc_90in <- timeROC(
T = sad_inf$dias_out,
delta = sad_inf$evento_ampu,
marker = sad_inf$san_elian, # 👈 ACÁ está San Elián
cause = 1,
times = 90,
iid = TRUE
)
table(sad_inf$evento_ampu)##
## 0 1 2
## 68 48 172roc_90in$AUC## NULLplot(roc_90in, time = 90)
roc_120in <- timeROC(
T = sad_inf$dias_out,
delta = sad_inf$evento_ampu,
marker = sad_inf$san_elian,
cause = 1,
times = 120,
iid = TRUE
)
roc_120in$AUC## NULLplot(roc_120in, time = 120)
roc_180in <- timeROC(
T = sad_inf$dias_out,
delta = sad_inf$evento_ampu,
marker = sad_inf$san_elian,
cause = 1,
times = 180,
iid = TRUE
)
roc_180in$AUC## NULLplot(roc_180in, time = 180)
#gráfico 3d
#scores <- factor.scores(datos1, Fa3, method = "tenBerge")
#df <- as.data.frame(scores$scores)
#install.packages("plotly")
#library(plotly)
#plot_ly(
# df,
# x = ~ML1,
# y = ~ML2,
# z = ~ML3,
# type = "scatter3d",
# mode = "markers",
# marker = list(size = 5)
#) %>%
# layout(
# scene = list(
# xaxis = list(title = "Factor 1"),
# yaxis = list(title = "Factor 2"),
# zaxis = list(title = "Factor 3")
# )
# )table(sad$estado_unif)##
## amputacion cerro_herida fallecio persiste
## 58 232 14 102describeBy(sad$san_elian, sad$estado_unif)##
## Descriptive statistics by group
## group: amputacion
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 58 22.5 3.22 23 22.6 2.97 15 29 14 -0.29 -0.33 0.42
## ------------------------------------------------------------
## group: cerro_herida
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 232 17.2 3.76 17 17.2 4.45 8 26 18 -0.19 -0.5 0.25
## ------------------------------------------------------------
## group: fallecio
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 14 19.3 4.75 20.5 19.3 5.19 11 27 16 -0.17 -1.32 1.27
## ------------------------------------------------------------
## group: persiste
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 102 18.2 3.58 18 18.3 2.97 9 26 17 -0.28 -0.28 0.35hist(sad$san_elian[sad$estado_unif=="amputacion"])
hist(sad$san_elian[sad$estado_unif=="cerro_herida"])
hist(sad$san_elian[sad$estado_unif=="fallecio"])
hist(sad$san_elian[sad$estado_unif=="persiste"])
shapiro.test(sad$san_elian[sad$estado_unif=="amputacion"])##
## Shapiro-Wilk normality test
##
## data: sad$san_elian[sad$estado_unif == "amputacion"]
## W = 1, p-value = 0.2shapiro.test(sad$san_elian[sad$estado_unif=="cerro_herida"])##
## Shapiro-Wilk normality test
##
## data: sad$san_elian[sad$estado_unif == "cerro_herida"]
## W = 1, p-value = 0.01shapiro.test(sad$san_elian[sad$estado_unif=="fallecio"])##
## Shapiro-Wilk normality test
##
## data: sad$san_elian[sad$estado_unif == "fallecio"]
## W = 0.9, p-value = 0.5shapiro.test(sad$san_elian[sad$estado_unif=="persiste"])##
## Shapiro-Wilk normality test
##
## data: sad$san_elian[sad$estado_unif == "persiste"]
## W = 1, p-value = 0.1qqnorm(sad$san_elian[sad$estado_unif=="amputacion"])
qqnorm(sad$san_elian[sad$estado_unif=="cerro_herida"])
qqnorm(sad$san_elian[sad$estado_unif=="fallecio"])
qqnorm(sad$san_elian[sad$estado_unif=="persiste"])
kruskal.test(sad$san_elian, sad$estado_unif)##
## Kruskal-Wallis rank sum test
##
## data: sad$san_elian and sad$estado_unif
## Kruskal-Wallis chi-squared = 76, df = 3, p-value <0.0000000000000002pairwise.wilcox.test(sad$san_elian, sad$estado_unif,"bonferroni")##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: sad$san_elian and sad$estado_unif
##
## amputacion cerro_herida fallecio
## cerro_herida < 0.0000000000000002 - -
## fallecio 0.1 0.5 -
## persiste 0.0000000001 0.1 1.0
##
## P value adjustment method: bonferronitapply(sad$san_elian,
sad$estado_unif,
quantile,
probs = c(0.25, 0.75),
na.rm = TRUE)## $amputacion
## 25% 75%
## 20.2 24.0
##
## $cerro_herida
## 25% 75%
## 15 20
##
## $fallecio
## 25% 75%
## 15 22
##
## $persiste
## 25% 75%
## 16 20#en pacientes isq
table(sad_isq$estado_unif)##
## amputacion cerro_herida fallecio persiste
## 46 66 8 49describeBy(sad_isq$san_elian, sad_isq$estado_unif)##
## Descriptive statistics by group
## group: amputacion
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 46 22.4 3.42 23 22.5 4.45 15 29 14 -0.18 -0.54 0.5
## ------------------------------------------------------------
## group: cerro_herida
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 66 18.9 2.96 19 19 2.97 12 24 12 -0.09 -0.82 0.36
## ------------------------------------------------------------
## group: fallecio
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 8 19.9 4.32 20.5 19.9 2.22 14 27 13 -0.04 -1.21 1.53
## ------------------------------------------------------------
## group: persiste
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 49 19.8 2.83 19 19.8 2.97 14 26 12 0.15 -0.51 0.4kruskal.test(sad_isq$san_elian, sad_isq$estado_unif)##
## Kruskal-Wallis rank sum test
##
## data: sad_isq$san_elian and sad_isq$estado_unif
## Kruskal-Wallis chi-squared = 27, df = 3, p-value = 0.000005pairwise.wilcox.test(sad_isq$san_elian, sad_isq$estado_unif,"bonferroni")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: sad_isq$san_elian and sad_isq$estado_unif
##
## amputacion cerro_herida fallecio
## cerro_herida 0.000004 - -
## fallecio 0.5 1.0 -
## persiste 0.000780 0.9 1.0
##
## P value adjustment method: bonferronitapply(sad_isq$san_elian,
sad_isq$estado_unif,
quantile,
probs = c(0.25, 0.75),
na.rm = TRUE)## $amputacion
## 25% 75%
## 20 24
##
## $cerro_herida
## 25% 75%
## 17 21
##
## $fallecio
## 25% 75%
## 17.8 22.0
##
## $persiste
## 25% 75%
## 18 22#infectados
table(sad_inf$estado_unif)##
## amputacion cerro_herida fallecio persiste
## 48 164 8 68describeBy(sad_inf$san_elian, sad_inf$estado_unif)##
## Descriptive statistics by group
## group: amputacion
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 48 23.1 3.04 23 23.2 2.97 15 29 14 -0.38 -0.17 0.44
## ------------------------------------------------------------
## group: cerro_herida
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 164 18.9 2.72 19 18.8 2.97 13 26 13 0.18 -0.56 0.21
## ------------------------------------------------------------
## group: fallecio
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 8 20.9 4.52 21.5 20.9 3.71 14 27 13 -0.29 -1.46 1.6
## ------------------------------------------------------------
## group: persiste
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 68 19.6 2.64 19 19.5 2.97 14 26 12 0.31 -0.5 0.32kruskal.test(sad_inf$san_elian, sad_inf$estado_unif)##
## Kruskal-Wallis rank sum test
##
## data: sad_inf$san_elian and sad_inf$estado_unif
## Kruskal-Wallis chi-squared = 58, df = 3, p-value = 0.000000000002pairwise.wilcox.test(sad_inf$san_elian, sad_inf$estado_unif,"bonferroni")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: sad_inf$san_elian and sad_inf$estado_unif
##
## amputacion cerro_herida fallecio
## cerro_herida 0.000000000001 - -
## fallecio 1.0 0.7 -
## persiste 0.000000102104 0.4 1.0
##
## P value adjustment method: bonferronitapply(sad_inf$san_elian,
sad_inf$estado_unif,
quantile,
probs = c(0.25, 0.75),
na.rm = TRUE)## $amputacion
## 25% 75%
## 21 25
##
## $cerro_herida
## 25% 75%
## 17 21
##
## $fallecio
## 25% 75%
## 18.8 23.5
##
## $persiste
## 25% 75%
## 18 22