knitr::opts_chunk$set(comment = NA, warning = F, message = F, fig.retina = 3)Ronda Mayo 2023
1 Análisis de % de absorción de agua
Importante antes de empezar
Crear el Proyecto de análisis y guardar el archivo Excel con los datos en la misma carpeta
1.1 Setup
No queremos ver los mensajes, comentarios ni los warning. La opción fig.retina = 3 es para obtener gráficos de alta resolución (archivo final más pesado).
1.2 Importar datos
library(readxl)
datos.ronda <- read_excel('datos.xlsx', sheet = 'ronda')
datos.ronda# A tibble: 8 × 5
Laboratorio DENSIDAD_REAL_ARIDO_SATURADO DENSIDAD_REAL_ARIDO_S…¹ DENSIDAD_NETA
<dbl> <dbl> <dbl> <dbl>
1 110 2620 2580 2670
2 120 2670 2620 2760
3 130 2650 2630 2700
4 140 2620 2580 2700
5 150 2620 2580 2690
6 160 2610 2550 2710
7 170 2630 2590 2700
8 180 2585 2518 2700
# ℹ abbreviated name: ¹DENSIDAD_REAL_ARIDO_SECO
# ℹ 1 more variable: ABSORCION_DE_AGUA <dbl>Si queremos darle formato a la tabla usamos el comando kable()del package knitr
library(knitr)
kable(datos.ronda)| Laboratorio | DENSIDAD_REAL_ARIDO_SATURADO | DENSIDAD_REAL_ARIDO_SECO | DENSIDAD_NETA | ABSORCION_DE_AGUA |
|---|---|---|---|---|
| 110 | 2620 | 2580 | 2670 | 1.30 |
| 120 | 2670 | 2620 | 2760 | 1.90 |
| 130 | 2650 | 2630 | 2700 | 1.10 |
| 140 | 2620 | 2580 | 2700 | 1.65 |
| 150 | 2620 | 2580 | 2690 | 1.65 |
| 160 | 2610 | 2550 | 2710 | 2.20 |
| 170 | 2630 | 2590 | 2700 | 1.65 |
| 180 | 2585 | 2518 | 2700 | 2.69 |
1.3 Selección de la variable a analizar
En este caso analizaremos el % de absorción de agua y lo guardaremos como x
x <- datos.ronda$ABSORCION_DE_AGUA
x[1] 1.30 1.90 1.10 1.65 1.65 2.20 1.65 2.691.4 EDA: Exploratory Data Analysis
1.4.1 Boxplot
boxplot(x)
1.4.2 Stripchart
stripchart(x, pch = 19, method = 'stack')
1.4.3 Density plot
plot(density(x), main = 'kernel density plot',
xlab = 'kg/m3')
1.4.4 QQ-plot de normalidad
qqnorm(x)
qqline(x)
1.4.5 Test de Shapiro de Normalidad
shapiro.test(x)
Shapiro-Wilk normality test
data: x
W = 0.94562, p-value = 0.66711.4.6 Test de Rosner para evaluar outliers
library(EnvStats)
rosnerTest(x)$all.stats i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
1 0 1.767500 0.5015903 2.69 8 1.839150 2.126645 FALSE
2 1 1.635714 0.3625308 2.20 6 1.556518 2.019969 FALSE
3 2 1.541667 0.2888194 1.10 3 1.529214 1.887145 FALSE1.4.7 Conclusiones del EDA
No hay evidencias de datos anómalos (outliers)
No es concluyente el test de normalidad debido al reducido conjunto de datos disponibles (asimetría en boxplot + p \(>\) 0,05 en test de Shapiro + QQ-plot no concluyente)
Datos repetidos (n = 3 resultados 1,65 %) ¿De nuevo los mismos laboratorios?
Kernel density plot muestra una moda relevante cercana a 1,65 %
1.5 Estimación del valor asignado \(x_{pt}\), su incertidumbre \(u(x_{pt})\) y el \(\sigma_{pt}\)
Generamos la tabla con los distintos algoritmos.
Recordemos de la sección Section 1.3 que hemos identificado la variable a analizar y que la hemos guardado como x
1.5.1 Cargar librerías y script
library(metRology)
library(revss)
library(robustbase)
source('QHampel.r')1.5.2 Número de datos
p <- length(x)1.5.3 Valor asignado \(x_{pt}\)
promedio <- mean(x)
mediana <- median(x)
algA.mu <- algA(x)$mu
QHampel.mu <- QHampel(x, 1:p)$x.star1.5.4 Dispersión
s <- sd(x)
s.mad <- mad(x)
algA.s <- algA(x)$s
nIQR <- 0.7413*IQR(x)
QHampel.s <- QHampel(x, 1:p)$s.star
Qn.s <- Qn(x)
algB.s <- robScale(x)1.5.5 Tabla con todos los estimadores
estimadores <- data.frame(estimador = c('Promedio-DS', 'Mediana-MAD',
'Algoritmo A', 'Q/Hampel',
'Qn', 'Algoritmo B', 'nIQR'),
centralidad = c(promedio,
mediana,
algA.mu,
QHampel.mu,
NA,
NA,
NA),
dispersion = c(s,
s.mad,
algA.s,
QHampel.s,
Qn.s,
algB.s,
nIQR))estimadores estimador centralidad dispersion
1 Promedio-DS 1.767500 0.5015903
2 Mediana-MAD 1.650000 0.4447800
3 Algoritmo A 1.747859 0.5234076
4 Q/Hampel 1.757888 0.5701429
5 Qn NA 0.5203340
6 Algoritmo B NA 0.3897328
7 nIQR NA 0.3057863Le damos formato con kable()
library(knitr) # ya la habíamos cargado
options(knitr.kable.NA = '') # Para que no aparezcan los NA
kable(estimadores, digits = 2) # con digits controlamos los decimales| estimador | centralidad | dispersion |
|---|---|---|
| Promedio-DS | 1.77 | 0.50 |
| Mediana-MAD | 1.65 | 0.44 |
| Algoritmo A | 1.75 | 0.52 |
| Q/Hampel | 1.76 | 0.57 |
| Qn | 0.52 | |
| Algoritmo B | 0.39 | |
| nIQR | 0.31 |
Exportamos la tabla con el comando write_xlsx del package writexl
library(writexl)
write_xlsx(estimadores, 'ronda_mayo_2023_absorcion_agua.xlsx')El archivo ronda_mayo_2023_absorcion_agua.xlsx quedó guarado en la misma carpeta del proyecto donde está guardado este análisis.
Importante
Siempre cerrar el archivo Excel antes de volver a compilar el archivo Quarto
1.5.6 Elección del valor asignado
Nuevamente, al igual que el análisis de densidad real del árido saturado, se observa que la mediana (línea roja) coincide con el máximo del kernel density:
plot(density(x), main = 'kernel density plot',
xlab = 'kg/m3')
abline(v = mediana, col = 'red')
Como referencia se grafica el promedio (línea negra):
plot(density(x), main = 'kernel density plot',
xlab = 'kg/m3')
abline(v = mediana, col = 'red')
abline(v = promedio) # el color por defecto es negro
Podemos asignar como valor asignado \(x_{pt}\) a la mediana de los datos, calculada en la seccion Section 1.5.3
x.pt <- mediana1.5.7 Elección de la dispersión \(\sigma_{pt}\)
Nuevamente nos guiaremos por la Norma Chilena NCh 1239-2009 y la ASTM C128 para obtener el \(\sigma_{pt}\) que usaremos en la evaluación de desempeño:
De acuerdo a la sección 10.4 de NCh 1239-2009, el criterio de aceptación entre duplicados relizados por el mismo analista en el mismo laboratorio es 0,30 % de agua.
En el caso de la ASTM C128 (norma de referencia de la NCh 1239-2009), se indican valores de desviación estándar intra-laboratorio y entre-laboratorios para Absorption, como se muestra en la siguiente tabla:
En este caso, el dato de interés es la *Standard Deviation Multilaboratory Precision$ de 0.23 %.
Se define el límite de reproducibilidad como \(R\):
Límite de reproducibilidad
La máxima diferencia tolerable entre duplicados obtenidos por dos laboratorios distintos utilizando el mismo método
\[R = 2.8\cdot \sigma_{R}\] De esta definición podemos obtener 2 estimadores de \(\sigma_{pt}\)
- De la NCh 1239-2009 \(\sigma_{pt} = 2*0.3/2.8 = 0.21\%\)
- De la ASTM C128 \(\sigma_{pt} = 0.23 \%\)
Ambos estimadores son muy consistentes.
Warning
El factor 2 que está multiplicando en la NCh es debido a que esta norma sólo indica límite de repetibilidad, es decir, 1 analista en el mismo laboratorio.
En este caso, para ser consistentes con el análisis de densidades, tomaremos como \(\sigma_{pt}\) el dato indicado por la NCh 1239-2009 de 0,21 %.
sigma.pt <- 2*0.3/2.8 # el dato se calcula con todos los decimales1.5.8 Estimación de la incertidumbre del valor asignado \(u(x_{pt})\)
Ecuación (6) de la guía ISO 13528:2022
\[u(x_{pt}) = 1.25\times \frac{s^{*}}{\sqrt{p}}\]
donde \(s^{*}\) es la desviación estándar robusta de los resultados y \(p\) es el número de datos.
Importante
No confundir la incertidumbre del valor asignado \(u(x_{pt})\) con el \(\sigma_{pt}\) ¡son cosas distintas!
De acuerdo al punto D.1.4.2 Nota 3 usaremos el estimador obtenido por el algoritmo B \(s^{*} = 0.39\,\,\%\)
\[\begin{eqnarray}u(x_{pt}) &=& 1.25\times \frac{s^{*}}{\sqrt{p}}\\ &=& 0.17\,\,\%\end{eqnarray}\]
el cual fue obtenido con el siguiente código:
s.star <- algB.s
u.x.pt <- 1.25*s.star/sqrt(p)
u.x.pt[1] 0.1722392Por lo tanto, la incertidumbre expandida al 95% de confianza es U.x.pt <- 2*u.x.pt
\(U(x_{pt}) = 2u(x_{pt}) = 0.34\,\, \%\)
Valor asignado \(\pm\) incertidumbre expandida al 95% de confianza
\(1.65 \pm 0.34\,\, \%\)
1.5.9 Criterio de incertidumbre límite del valor asignado 9.2.1
Ecuación (10):
\[u(x_{pt}) < 0.3\sigma_{pt}\]
En este caso \(u(x_{pt}) = 0.17 < 0.21 = 0.3\sigma_{pt}\), por lo tanto, se cumple con el criterio especificado en 9.2.1 de la ISO 13528:2022
1.6 Evaluación del desempeño con \(z_{\text{score}}\)
Utilizando como valor asignado la mediana \(x_{pt} = 1.65\,\, \%\) y como \(\sigma_{pt} = 0.21 \,\, \%\) calculamos el estadístico de desempeño \(z_{\text{score}}\) para cada laboratorio:
\[z_{\text{score}} = \frac{\overline{x} - x_{pt}}{\sigma_{pt}}\]
z.score <- (x - x.pt)/sigma.pt
z.score <- data.frame(Laboratorio = datos.ronda$Laboratorio, z = z.score)
kable(z.score, digits = 2) # cargando el package knitr con library(knitr)| Laboratorio | z |
|---|---|
| 110 | -1.63 |
| 120 | 1.17 |
| 130 | -2.57 |
| 140 | 0.00 |
| 150 | 0.00 |
| 160 | 2.57 |
| 170 | 0.00 |
| 180 | 4.85 |
Graficamos el \(z_\text{score}\) considerando como límites de proficiencia \(\pm 3\)
lim.sup <- 3 # límite superior de proficiencia
lim.inf <- -3 # límite inferior de proficiencia
plot(z.score$Laboratorio, # eje X
z.score$z, # eje Y
pch = 19, # tipo de puntos
ylim = c(-5, 5), # límites de y
ylab = 'Z-score', # etiqueta eje Y
xlab = 'Laboratorio', # etiqueta eje X
main = 'Z-score % de absorción de agua') # título prinicpal
abline(h = lim.sup, col = 'red')
abline(h = lim.inf, col = 'red')
abline(h = 0, lty = 2)