Requiere los paquetes “psych”" y “ggplot2”.
#Especificar los parametros: x (primera medicion), y (segunda medicion), title.name (titulo del grafico, por defecto "Grafico de Bland-Altman"), x.name (titulo eje X, por defecto "Media"), y (titulo eje Y, por defecto "Diferencia")
FiabilidadRubenFMat <- function(x, y, title.name = "Grafico de Bland-Altman", x.name = "Media", y.name = "Diferencia"){
#Funciones iniciales necesarias para los analisis
library(psych)
library(ggplot2)
DataFrameFiabilidadRubenFMat <- as.data.frame(cbind(x, y))
#Descriptivos de los datos para usos posteriores
media1 <- mean(DataFrameFiabilidadRubenFMat$x)
media2 <- mean(DataFrameFiabilidadRubenFMat$y)
#Estadisticos para los graficos de Bland-Altman
diferencia <- x - y
media <- (x + y)/2
diferenciamedia <- mean(diferencia)
LimiteSupBA <- diferenciamedia + 2*sd(diferencia)
LimiteInfBA <- diferenciamedia - 2*sd(diferencia)
DataFrameBlandAltman <- as.data.frame(cbind(media, diferencia))
CorrelacionMediaDiferencia <- cor.test(diferencia, media, alternative = "two.sided", conf.level = 0.95)
CorrelacionMediaDiferenciaAbs <- cor.test(abs(diferencia), media, alternative = "two.sided", conf.level = 0.95)
#Calculo de los ICC con el paquete psych
ICCFMat <- ICC(DataFrameFiabilidadRubenFMat, missing = TRUE, alpha = .05, lmer = FALSE, check.keys = FALSE)
#Calculo del SEM (con la formula de raiz quadrada del cuadrado medio del error del ANOVA) y MDC (al 95% de confianza) y del porcentaje de SEM y MDC con respecto a la media muestral
SEM <- sqrt(ICCFMat$stats[3, "Residual"])
PorcentajeSEM <- SEM*100/mean(media)
MDC <- SEM*sqrt(2)*1.96
PorcentajeMDC <- MDC*100/mean(media)
BlandAltmanPlot <- ggplot(DataFrameBlandAltman, aes(media, diferencia)) + geom_point(size = 2) + geom_smooth(method = "lm", size = 1.5, colour = "turquoise4") + geom_hline(yintercept = diferenciamedia, size = 1.5) + geom_hline(yintercept = LimiteSupBA, linetype = "dashed", size = 1.5, colour = "tomato") + geom_hline(yintercept = LimiteInfBA, linetype = "dashed", size = 1.5, colour = "tomato") + ggtitle(title.name) + theme(plot.title = element_text(lineheight = .8, face = "bold", hjust = 0.5, size = 20)) + theme(axis.title.x = element_text(face = "bold", size = 15), axis.text.x = element_text(size = 13, face = "bold")) + theme(axis.title.y = element_text(face = "bold", size = 15), axis.text.y = element_text(size = 13, face = "bold")) + labs(x = x.name, y = y.name)
Lista <- list(SEM, PorcentajeSEM, MDC, PorcentajeMDC, LimiteSupBA, LimiteInfBA, diferenciamedia, CorrelacionMediaDiferencia, CorrelacionMediaDiferenciaAbs, ICCFMat, ICCFMat$stats, BlandAltmanPlot)
names(Lista) <- c("SEM", "Porcentaje SEM", "MDC", "Porcentaje MDC", "Limite Superior Acuerdo Bland Altman", "Limite Inferior Acuerdo Bland Altman", "Error Sistematico", "Correlacion Diferencia y Media", "Correlacion Diferencia (Absoluta) y Media", "ICC", "ANOVA", "Grafico Bland Altman")
Lista
}
Ejemplo con datos simulados.
Medicion1 <- rnorm(n = 50, mean = 40, sd = 8)
Medicion2 <- Medicion1 + rnorm(n = 50, mean = 0, sd = 4)
FiabilidadRubenFMat(Medicion1, Medicion2, title.name = "Prueba Grafico de Bland-Altman", x.name = "Prueba titulo eje X", y.name = "Prueba titulo eje Y")
##
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
##
## %+%, alpha
## $SEM
## [1] 2.873367
##
## $`Porcentaje SEM`
## [1] 7.149656
##
## $MDC
## [1] 7.964566
##
## $`Porcentaje MDC`
## [1] 19.81783
##
## $`Limite Superior Acuerdo Bland Altman`
## [1] 7.993222
##
## $`Limite Inferior Acuerdo Bland Altman`
## [1] -8.260993
##
## $`Error Sistematico`
## [1] -0.1338855
##
## $`Correlacion Diferencia y Media`
##
## Pearson's product-moment correlation
##
## data: diferencia and media
## t = -3.1374, df = 48, p-value = 0.002911
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6197069 -0.1515719
## sample estimates:
## cor
## -0.4125157
##
##
## $`Correlacion Diferencia (Absoluta) y Media`
##
## Pearson's product-moment correlation
##
## data: abs(diferencia) and media
## t = 1.4822, df = 48, p-value = 0.1448
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.07342491 0.46071847
## sample estimates:
## cor
## 0.2091983
##
##
## $ICC
## Call: ICC(x = DataFrameFiabilidadRubenFMat, missing = TRUE, alpha = 0.05,
## lmer = FALSE, check.keys = FALSE)
##
## Intraclass correlation coefficients
## type ICC F df1 df2 p lower bound upper bound
## Single_raters_absolute ICC1 0.88 16 49 50 2.0e-18 0.81 0.93
## Single_random_raters ICC2 0.88 16 49 49 5.8e-18 0.80 0.93
## Single_fixed_raters ICC3 0.88 16 49 49 5.8e-18 0.80 0.93
## Average_raters_absolute ICC1k 0.94 16 49 50 2.0e-18 0.89 0.97
## Average_random_raters ICC2k 0.94 16 49 49 5.8e-18 0.89 0.97
## Average_fixed_raters ICC3k 0.94 16 49 49 5.8e-18 0.89 0.96
##
## Number of subjects = 50 Number of Judges = 2
## $ANOVA
## subjects Judges Residual
## df 4.900000e+01 1.00000000 49.000000
## SumSq 6.472056e+03 0.44813335 404.555526
## MS 1.320828e+02 0.44813335 8.256235
## F 1.599794e+01 0.05427817 NA
## p 5.836990e-18 0.81674993 NA
##
## $`Grafico Bland Altman`
