Reliability Intraoral Scanner

Packages

library(tidyverse)
library(janitor)
library(irr)

readxl::excel_sheets("reliability_data.xlsx")

[1] "Sheet1" "Sheet3" "Sheet2"

Original messy dataset

head(df)

Clean names and empty cols

df <- df  %>% 
        janitor::clean_names() %>% 
        janitor::remove_empty_cols()

Less messy

head(df)

Create two datasets, for min and max

min <- df %>% 
        select(min:d_21_l) %>% 
        mutate(type = "comparison") %>% 
        gather(d_13_d:d_21_l, key = "tooth", value = "difference")
colnames(min)[colnames(min) == "min"] <- "comparison"
max <- df %>% 
        select(max:d_21_l_1) %>% 
        mutate(type = "comparison") %>% 
        gather(d_13_d_1:d_21_l_1, key = "tooth", value = "difference")
colnames(max)[colnames(max) == "max"] <- "comparison"
max$tooth = substr(max$tooth,1,nchar(max$tooth) -2)

Join datasets

df2 <- bind_rows(min, max) #create a semi tidy dataset
df2 <- df2 %>% 
        select(-type)

# rm(df, min, max) #delete unused datasets

Dataset ready for work

head(df2)

Descriptive stats

df2 %>% 
        select(comparison:difference) %>% 
        group_by(comparison) %>% 
        summarise_each (funs(min, 
                      q25 = quantile(., 0.25), 
                      median = median, 
                      q75 = quantile(., 0.75), 
                      max = max,
                      mean = mean, 
                      sd = sd), 
                      difference)

`summarise_each()` is deprecated.
Use `summarise_all()`, `summarise_at()` or `summarise_if()` instead.
To map `funs` over a selection of variables, use `summarise_at()`

abs difference

df2 <- df2 %>% 
        mutate("abs_diff" = abs(difference)) #absolute difference
df2 <- df2 %>% 
        mutate("expect_diff" = 0) #expected difference = 0

Distribution

df2 %>% 
        ggplot(aes(x = abs_diff)) +
        geom_histogram(bins = 5) + 
        scale_x_log10()

There is an outlier

which.max(df2$abs_diff)

[1] 830

Delete the outlier

df2 <- df2[-830, ]

df2 %>% ggplot(aes(difference)) + 
        geom_histogram(bins = 10) + 
        labs(title = "Distribution of the differences between measurements", 
             x = "Difference", 
             y = "Counts") + 
        theme_minimal()

df3 <- df2 %>% 
        group_by(tooth, comparison) %>% 
        summarise("Average" = mean(difference)) %>% 
        spread(comparison, value = Average)

Agreement

df3 %>% 
        ungroup() %>% 
        select(1--10:9--10) %>% 
        icc(., model = "twoway", type = "agreement")

 Single Score Intraclass Correlation

   Model: twoway 
   Type : agreement 

   Subjects = 12 
     Raters = 20 
   ICC(A,1) = 0.149

 F-Test, H0: r0 = 0 ; H1: r0 > 0 
  F(11,226) = 4.47 , p = 4.22e-06 

 95%-Confidence Interval for ICC Population Values:
  0.056 < ICC < 0.377

There are significant differences between the measurements across raters

Absolute values matrix

df3_abs <- df3 %>% 
        ungroup() %>% 
        select(1--10:9--10) %>% 
        abs(.) 

df2 %>% 
        ggplot(aes(x = tooth, y = abs_diff)) +
        
        geom_boxplot() + 
        geom_jitter(alpha = 0.1) + 
        coord_flip() +
        labs( title = "Measurements per sample", 
              y = "Absolute difference", 
              x = "Sample and surface") + 
        theme_minimal()

Seems that the machine has some thresholds for detection, since there are no 0.01

set.seed(3000)
df2 <- df2 %>% 
        mutate(expected_values = rnorm(1079, mean = 0, sd = 0.01))

df2 %>% 
        select(difference, expected_values) %>% 
        gather("type", "value", 1:2) %>% 
        ggplot(aes(value, fill = type)) +
        geom_density(alpha = 0.4) + 
        labs(title = "Expected random variance vs real differences", 
             subtitle = "3 samples, 4 measurements each.", 
             x = "Difference", 
             y = "Count", 
             fill = "Key") + 
        theme_minimal()

df2 %>% 
        ggplot(aes(x = difference, y = expected_values)) +
        geom_point() +
        stat_smooth() + 
        labs(title = "Scatterplot expected vs real differences", 
             y = "Expected by random", 
             x = "Real differences") +
        theme_minimal()

library(BlandAltmanLeh)
ba_plot <- bland.altman.plot(df2$difference, df2$expected_values, graph.sys="ggplot2", conf.int=.95)
ba_plot + theme_minimal() + 
        labs(title = "Altman-Bland plot", 
             subtitle = "Expected vs real differences", 
             y = "Difference expected - real", 
             x = "Mean of measurements")

Read Understanding Bland Altman analysis

Is any statistical difference between expected and real differences?

Against expected distribution

t.test(df2$difference, df2$expected_values)

    Welch Two Sample t-test

data:  df2$difference and df2$expected_values
t = -0.26614, df = 1344.2, p-value = 0.7902
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.002066110  0.001572482
sample estimates:
   mean of x    mean of y 
0.0004355885 0.0006824026 

df2 %>% 
        select(difference, expected_values) %>% 
        gather() %>% 
        ggplot(aes(x = key, y = value)) + 
        geom_boxplot() + 
        theme_minimal()

As one sample

t.test(df2$difference, mu = 0)

    One Sample t-test

data:  df2$difference
t = 0.49828, df = 1078, p-value = 0.6184
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.001279716  0.002150893
sample estimates:
   mean of x 
0.0004355885 

Are normally distributed the differences?

shapiro.test(df2$difference)

    Shapiro-Wilk normality test

data:  df2$difference
W = 0.94936, p-value < 2.2e-16

library("ggpubr")
ggqqplot(df2$difference, 
         ylab = "Measured differences",
         ggtheme = theme_minimal())

Compare with expected values

ggqqplot(df2$expected_values, 
         ylab = "Expected values",
         ggtheme = theme_minimal())

shapiro.test(df2$expected_values)

    Shapiro-Wilk normality test

data:  df2$expected_values
W = 0.99895, p-value = 0.8007