# Avaliação da Capacidade de Sistemas de Medição
#
# São Paulo, 25.05.2022
# Walter Ponge-Ferreira
#
# Ref.:
#
# Costa, A. F. B.; Epprecht, E. K.; Carpinetti, L. C. R.
# Controle Estatístico de Qualidade
# 2 ed., Editora Atlas, São Paulo, 2003
# Cap. 5: Avaliação de Sistema de Medição
# p.139-158
#
# Montgomery, Douglas C.
# Introdução ao Controle Estatístico da Qualidade
# 7 ed., RTC, Rio de Janeiro, 2016
# 8.7 Estudos sobre a Capacidade de um Medidor e de um Sistema de Medidas 
# p. 269-281
#
# Cano, E. L.; Maguerza, J. M.; Redchuk, A.
# Six Sigma with R
# Statiscal Engineering for Process Improvement
# Springer,  New York, 2012
# Chapter 5 - Measurement System Analysis with R
# p. 79-90
#

setwd("~/PastaTecnica/R/Capacidade de Medição")

library(qcc)
library(SixSigma)
library(UsingR)
library(readxl)
library(ggplot2)
library(formattable)

rm(list = ls())

# Dados devem estar no formato tidy
dados <- read_excel("Exemplo 1.xlsx")

class(dados)
str(dados)
is(dados)
names(dados)
dados$peça <- factor(dados$peça)
dados$operador <- factor(dados$operador)
dados$medida <- factor(dados$medida)
str(dados)
names(dados)

# Análise no agregado
attach(dados)
  mean(valor)
  sd(valor)
  fivenum(valor)
detach(dados)
hist(dados$valor, main = "Histograma Agragado", xlab = "Valor") 
qqnorm(dados$valor); grid()
plot(sort(dados$valor), ylab = "Valor"); grid()
DOTplot(dados$valor); grid()
plot(dados$valor, ylab = "Valor"); grid()
ggplot(dados, aes(x = peça, y = valor, col = operador)) + geom_point()
ggplot(dados, aes(x = peça, y = valor, col = operador)) + geom_point(size = 1) + facet_grid(~ operador)
boxplot(valor ~ operador, data = dados)
boxplot(valor ~ peça, data = dados); grid()
ggplot(dados, aes(x = valor, col = operador)) + geom_histogram() + facet_grid(~ operador)

# Número de itens, fatores e repetições
n <- nlevels(dados$peça)
m <- nlevels(dados$operador)
r <- nlevels(dados$medida)
print(paste("Nº Itens =",n,"- Nº Fatores =",m,"- Nº Repetições =",r))

tabela <- table(dados[c("operador","peça")])
print(tabela)
barplot(tabela, xlab = "Nº da Peça", ylab = "Medidas")
ggplot(dados, aes(x = peça, col = operador)) + geom_bar()

xtabs(valor/(n*r) ~ operador, data = dados)
xtabs(valor/(m*r) ~ peça, data = dados)

# -------------------------------------------------------
# Avaliação Simplificada pelo Método de Costa et al. e Montgomery

tapply(dados$valor, dados$operador, mean)
tapply(dados$valor, dados$operador, sd)
tapply(dados$valor, dados$operador, summary)
x <- split(dados$valor, dados$operador)
boxplot(x); grid(); title("Medidas por Operador")
x <- as.data.frame(x)
names(x) <- c("Operador 1", "Operador 2", "Operador 3")
formattable(x)

tapply(dados$valor, dados$peça, mean)
tapply(dados$valor, dados$peça, sd)
tapply(dados$valor, dados$peça, summary)
y <- split(dados$valor, dados$peça)
boxplot(y); grid(); title("Medidas por Peça")
y <- as.data.frame(y)
names(y) <- c("Peça 1","Peça 2","Peça 3","Peça 4","Peça 5",
              "Peça 6","Peça 7","Peça 8","Peça 9","Peça 10")
formattable(y)

# Agregando Médias e Amplitudes - Solução 1
aggregate(valor ~ operador, data = dados, mean)
aggregate(valor ~ operador, data = dados, sd)
media <- aggregate(dados$valor, list(Peça = dados$peça, Operador = dados$operador), mean)
amplitude <- aggregate(dados$valor, list(Peça = dados$peça, Operador = dados$operador), diff)
formattable(media)
Xbar <- aggregate(media$x, list(Operador = media$Operador), mean)
print(Xbar)
RXbar <- diff(range(Xbar$x))
print(RXbar)
amplitude$x <- abs(amplitude$x)
formattable(amplitude)
Rbar <- aggregate(amplitude$x, list(Operador = amplitude$Operador), mean)
print(Rbar)
Rbarbar <- mean(Rbar$x)
print(Rbarbar)

# Agregando Médias e Amplitudes - Solução 2
ind <- dados$medida == 1
tabela1 <- data.frame(peça = dados$peça[ind],
                     operador = dados$operador[ind])
nomes <- c("peça","operador")
for(i in 1:r) {
  tabela1 <- cbind(tabela1, dados$valor[dados$medida == i])
  nomes <- c(nomes, paste("Medida ",i))
  }
names(tabela1) <- nomes
head(tabela1)
formattable(tabela1)

tabela2 <- data.frame(peça = dados$peça[ind],
                      operador = dados$operador[ind])
tabela2$media <- apply(tabela1[,3:(3+r-1)],1,mean)
tabela2$amplitude <- abs(apply(tabela1[,3:(3+r-1)],1,diff))
head(tabela2)
formattable(tabela2)
summary(tabela2)

qcc(tabela1[,3:(2+r)], type = "xbar")
qcc(tabela1[,3:(2+r)], type = "R")

Rbarbar <- mean(tabela2$amplitude)
print(Rbarbar)

Xbar <- c()
for (i in 1:m) {
  Xbar[i] <- mean(tabela2$media[tabela2$operador == i])
}
print(Xbar)
RXbar <- diff(range(Xbar))
print(RXbar)

# Variabilidade de repetição
sigma_repe <- Rbarbar/ss.cc.getd2(r)
print(sigma_repe)

# Variabilidade de reprodutibilidade
sigma_repro <- sqrt((RXbar/ss.cc.getd2(m))^2 - sigma_repe^2/(r*n))
print(sigma_repro)

# Variabilidade de medição
sigma_med <- sqrt(sigma_repe^2 + sigma_repro^2)
print(sigma_med)

# Variabilidade total
sigma_total <- sd(dados$valor)
print(sigma_total)

# Métricas de avaliação

# Fração total
rho_med <- (sigma_med/sigma_total)^2
rho_proc <- 1 - rho_med
print(c(rho_med,rho_proc))

# Razão Sinal-Ruído (RSR)
RSR <- sqrt(2*rho_proc/rho_med)
print(RSR)

# Razão de Discriminação (RD)
RD <- (1+rho_proc)/(1-rho_proc)
print(RD)

# Razão da precisão para a tolerência - Porcentagem de Tolerância
USL <- 20.5
LSL <- 19.5
k <- 6
PT <- k*sigma_med/(USL - LSL)
print(PT)

# ----------------------------------------------------------------------------
# Análise de Variância
dados <- read_excel("Exemplo 1.xlsx")

modelo1 <- lm(valor ~ peça + operador + peça*operador, data = dados)
anova(modelo1)
summary(modelo1)
coefficients(modelo1)
plot(modelo1)

modelo2 <- lm(valor ~ peça + operador, data = dados)
anova(modelo2)
summary(modelo2)
plot(modelo2)

# ----------------------------------------------------------------------------
# Avaliação de Sistema de Medição pela biblioteca SixSigma

rr <- ss.rr(var = valor, part = peça, appr = operador,
            data = dados,
            main = "Avaliação do Sistema de Medição",
            sub = "Exemplo 5.1")

# ***************************************************************************
# Exemplo 5.1 - Cano, Moguerza e Redchuk

rm(list=ls())

# Dados devem estar no formato tidy
dados <- read_excel("Exemplo 5-1.xlsx")

class(dados)
dados$Voltmeter <- factor(dados$Voltmeter)
dados$Battery <- factor(dados$Battery)
dados$Run <- factor(dados$Run)
str(dados)
formattable(dados)

# Número de itens, fatores e repetições
n <- nlevels(dados$Battery)
m <- nlevels(dados$Voltmeter)
r <- nlevels(dados$Run)
print(paste("Nº Itens =",n,"- Nº Fatores =",m,"- Nº Repetições =",r))

mod <- lm(Voltage ~ Battery + Voltmeter + Battery * Voltmeter, data = dados)
res <- anova(mod)
formattable(res)
str(res)

names(res)
res$"Mean Sq"

sigma_repe <- sqrt(res$"Mean Sq"[4])
print(sigma_repe^2)

sigma_fator <- sqrt((res$"Mean Sq"[2] - res$"Mean Sq"[3])/(n*r))
print(sigma_fator^2)

if (res$"Mean Sq"[3] >= res$"Mean Sq"[4]){
  sigma_interacao <- sqrt((res$"Mean Sq"[3] - res$"Mean Sq"[4])/(r))  
} else {
  sigma_interacao <- 0.0
}
print(sigma_interacao^2)

sigma_repro <- sqrt(sigma_fator^2 + sigma_interacao^2)
print(sigma_repro^2)

sigma_med <- sqrt(sigma_repe^2 + sigma_repro^2)
print(sigma_med^2)

sigma_proc <- sqrt((res$"Mean Sq"[1] - res$"Mean Sq"[3])/(m*r))
print(sigma_proc^2)

sigma_total <- sqrt(sigma_med^2 + sigma_proc^2)
print(sigma_total^2)

var(dados$Voltage)

# Six Sigma

my.rr <- ss.rr(var = Voltage,
               part = Battery,
               appr = Voltmeter,
               data = dados,
               main = "Six Sigma Gage R&R Measure",
               sub = "Batteries Project MSA")

# FIM