Proyecto 3 - Estadística avanzada.

Introducción.

Para la realización de este proyecto, se pretende reproducir de manera parcial (de acuerdo a los datos disponibles) los resultados del artículo:

Giavarina D. Understanding Bland Altman analysis. Biochem Med (Zagreb). 2015;25(2):141-151. Published 2015 Jun 5. doi:10.11613/BM.2015.015

---

El análisis de diferencias: El método de Bland y Altman.

Primero, cargamos los datos y los visualizamos:

# Cargamos los datos
df <- read.csv('datos.csv')

# Visualizamos estructura del dataset
head(df, n = 8)

# Extraemos variables 
x <- df$Method.A
y <- df$Method.B

# Graficamos los datos
plot(x, y, type = 'p', 
     main = 'Scatter plot of measurements', 
     col = 'blue', xlab = 'Method A',
     ylab = 'Method B')

# Etiquetas de líneas.
legend("topleft", c("Data"), fill=c("blue"))

# Agregamos una cuadrícula.
grid()

Construimos una regresión lineal, para ver si desde ahí podemos observar si los métodos de medición están de acuerdo:

# Realizamos la regresion lineal
modelo <- lm(Method.B ~ Method.A, data = df)

# El resultado
summary(modelo)

Call:
lm(formula = Method.B ~ Method.A, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-93.615 -11.908   0.052  12.314  55.504 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.0163     8.8667   1.355    0.186    
Method.A      1.0416     0.0181  57.561   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 32.49 on 28 degrees of freedom
Multiple R-squared:  0.9916,    Adjusted R-squared:  0.9913 
F-statistic:  3313 on 1 and 28 DF,  p-value: < 2.2e-16

Gráficamos este resultado de regresión con 99% de confianza:

# Definimos intervalo de interes
interval <- seq(min(x), max(x), length.out = 184)

# Calculamos las estimaciones para cada punto en el intervalo
mu <- predict(modelo, newdata = data.frame('Method A' = interval), 
              interval = "confidence", level = 0.99)

# Creamos gráfica
plot(x, y, type = 'p', 
     main = 'Confidence interval of 99% for linear regression', 
     col = 'blue', xlab = 'Method A',
     ylab = 'Method B')

# Coloreamos zona de confianza
polygon(c(rev(interval), interval),
        c(rev(mu[ ,3]), mu[ ,2]), col = 'grey',
        border = NA)

# Graficamos los datos
lines(x, y, type = 'p', col = 'blue')

# Graficamos mejor linea
abline(modelo)

# Etiquetas de líneas.
legend("topleft", c("Data"), fill=c("blue"))

# Agregamos una cuadrícula.
grid()

Ahora vamos a agregar las columnas necesarias al dataframe para poder construir las gráficas A&B. Primero, la columna de las diferencias:

df$diff <- df$Method.A - df$Method.B
head(df, 10)

NA

En teoría, si dos métodos estuvieran completamente de acuerdo, podríamos esperar que la diferencia entre sus mediciones sean siempre cere, y estén siempre de acuerdo). El problema de hacerlo de esta manera, es que no se considera el error de medición que pudiera existir en cada método.

Para considerar esto, debemos obtener el promedio de las diferencias de ambos métodos, ya que si la diferencia de cada medición es únicamente al error del método, entonces este promedio debería ser cercano a 0:

mean <- mean(df$diff)
mean

[1] -27.16667

En nuestro caso, el promedio de las diferencias no es cero, por lo que, en promedio, el método B mide \(27.17\) unidades más que el método A, en cada meidición. A esto le consideramos sesgo.

Si estamos seguros de que uno de los métodos está correcto, entonces con este promedio de diferencias podemos “ajustar” el otro método, pero en caso de que no haya ninguno que funcione como referencia, debemos investigar mas al respecto, obteniendo el promedio de las mediciones.

df$mean <- (df$Method.A + df$Method.B)/2
head(df, 8)

A partir de aquí, podemos construir la gráfica A&B, que es simplemente las diferencias entre los métodos contra el promedio entre los métodos. Esto nos permite investigar cualquier relación posible entre el error de medición y el valor real.

# Extraemos variables 
x <- df$mean
y <- df$diff

# Graficamos los datos
plot(x, y, type = 'p', 
     main = 'A&B plot', 
     col = 'blue', xlab = 'Mean of methods',
     ylab = 'Difference of methods')

abline(h = mean, col = 'red')

# Etiquetas de líneas.
legend("top", c("Data", "Mean (Bias)"), fill=c("blue", "red"))

# Agregamos una cuadrícula.
grid()

La línea roja, que corresponde con el promedio de las diferencias, es el bias respecto a la no diferencia entre los métodos. Se puede observar que las mediciones arriba de 200 son las responsables de este bias negativo.

Podemos hacer otra regresión lineal pero esta vez considerando estos ejes y datos, para encontrar alguna relación entre la diferencia de las mediciones y el valor real (que estamos estimando con el promedio de los métodos):

# Realizamos la regresion lineal
modelo <- lm(diff ~ mean, data = df)

# El resultado
summary(modelo)

Call:
lm(formula = diff ~ mean, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-54.536 -12.483  -1.058  10.294  94.298 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -10.14690    8.74881   -1.16   0.2559  
mean         -0.04505    0.01733   -2.60   0.0147 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 31.79 on 28 degrees of freedom
Multiple R-squared:  0.1945,    Adjusted R-squared:  0.1657 
F-statistic:  6.76 on 1 and 28 DF,  p-value: 0.01472

Visualizamos esta refresión con el intervalo de confianza al 0.95:

# Definimos intervalo de interes
interval <- seq(0, 1000, length.out = 184)

# Calculamos las estimaciones para cada punto en el intervalo
mu <- predict(modelo, newdata = data.frame("mean" = interval), 
              interval = "confidence", level = 0.95)

# Creamos gráfica
plot(x, y, type = 'p', 
     main = 'Confidence interval of 95% for linear regression \nin A&B plot', 
     col = 'blue', xlab = 'Mean of methods',
     ylab = 'Differences of methods')

# Coloreamos zona de confianza
polygon(c(rev(interval), interval),
        c(rev(mu[ ,3]), mu[ ,2]), col = 'grey',
        border = NA)

# Graficamos los datos
lines(x, y, type = 'p', col = 'blue')

# Graficamos mejor linea
abline(modelo)

# Gráfica de promedio (bias)
abline(h = mean, col = 'red')

# Etiquetas de líneas.
legend("top", c("Data", "Mean (Bias)"), fill=c("blue", "red"))

# Agregamos una cuadrícula.
grid()

Podemos ver que existe una tendencia negativa en la diferencia entre los métodos cuando el valor real estimado crece.

Nos interesa ahora obtener límites de acuerdo, esto es, un intervalo en el que podrémos decir que las diferencias entre las mediciones de ambos métodos se encuentran. Para esto, ya conocemos el promedio de las diferencias, obtenemos también la desviación estándar:

std <- sd(df$dif)
std

[1] 34.80595

Ahora consruimos un intervalo de confianza. Podríamos asumir que la distribución de la diferncia de los métodos es normal; sin embargo, debemos revisar que en efecto los datos sean normales. Para esto, se puede utilizar alguno de los múltiples test existentes. En este caso, mostraremos el histograma correspondiente y realizaremos una inspección visual y también con el test de Shapiro-Wilk:

differences <- df$diff
hist(differences)

shapiro.test(differences)

    Shapiro-Wilk normality test

data:  differences
W = 0.97956, p-value = 0.8137

Por lo tanto, podemos asumir que es normal (p = 0.814). Finalmente, construimos los intervalos de confianza al 95%:

up <- mean + 1.96*std
low <- mean - 1.96*std

print(paste(low, up))

[1] "-95.3863249384042 41.0529916050708"

Finalmente, nuestra gráfica A&B completa es:

# Extraemos variables 
x <- df$mean
y <- df$diff

# Graficamos los datos
plot(x, y, type = 'p', 
     main = 'A&B plot', 
     col = 'blue', xlab = 'Mean of methods',
     ylab = 'Difference of methods',
     ylim = c(-110, 50))

abline(h = mean, col = 'red')

abline(h = up, col = 'green')
abline(h = low, col = 'green')

# Etiquetas de líneas.
legend("top", c("Data", "Mean (Bias", "Limit"), fill=c("blue", "red", "green"))

# Agregamos una cuadrícula.
grid()

Por ultimo, es importante mencionar que el sistema de gráfica A&B no dice si el acuerdo es suficente, únicamente cuantifica el sesgo del rango de acuerdo, donde el 95% de las dfierences entre los métodos se encontrará. Se pueden realizar pruebas de significancia sobre este sesgo.

Referencias.

1. Gallagher R, McKinley S, Dracup K. Predictors of women’s attendance at cardiac rehabilitation programs. Prog Cardiovasc Nurs. 2003 Summer;18(3):121-6. doi: 10.1111/j.0889-7204.2003.02129.x. PMID: 12893973.
2. Wayne W. Daniel y Chad L. Cross, Biostatistics: A Foundation for Analysis in the Health Sciences, 10° edition, Wiley.

---

Anexos.

Diapositivas de presentación de proyecto.

Se pueden encontrar en este link.