Carta de control beta y rendimiento de procesos de manufactura

Rendimiento de la manufactura de comprimidos

library(tidyverse)
library(knitr)
library(fitdistrplus)
## library(bcc)

Introducción

La preparación de este docuento ha sido motivada por el artículo de Zagar y Mihelic (2022) Big data collection in pharmaceutical manufacturing and its use for product quality prediction y la disponibilidad de los datos referidos en el artículo a través de Figshare (acceso directo). La colección de datos ocupa 345 MB, y tienen un gran interés en el análisis de datos y validación de la manufactura de comprimidos, por lo que es de agradecer la disponibilidad de los autores a compartirlos con la comunidad académica.

Los datos recopilan información relevante sobre la manufactura de 1005 lotes de comprimidos recubiertos de un fármaco anticolesterilémico (no declarado) que contienen 5, 10, 20 y 40 mg de principio activo, elaborados por compresión directa, utilizando siempre la misma composición y ajustando la masa de los comprimidos para obtener las 4 formulaciones. Los lotes fueron elaborados entre noviembre de 2018 y abril de 2021. Los ficheros de datos disponibles son:

• laboratory.csv. Resume 55 parámetros de calidad clasificados con el criterio siguiente (descripción de los campos en la tabla 4 del artículo de Zagar y Mihelic):
  • genealogía del proceso de manufactura y materias primas,
  • parámetros de calidad de las materias primas,
  • parámetros de calidad de productos intermedios de la elaboración de los comprimidos,
  • parámetros de calidad del producto acabado.
• process.csv. Atributos derivados de los datos obtenidos en tiempo real durante la compresión (descripción de los campos en la tabla 6 del artículo de Zagar y Mihelic).
• 1.csv … 25.csv. Datos obtenidos en tiempo real durante la compresión de los 1005 lotes; la frecuencia de muestreo publicada es de 1 registro cada 10 segundos (10 campos por registro; descripción de los campos en la tabla 5 del artículo de Zagar y Mihelic).

Rendimiento de los procesos de manufactura

Por rendimiento del proceso se entiende el tamaño del lote obtenido (número de unidades, cantidad total, etc.) expresado como porcentaje del tamaño teórico. Las cartas de control del rendimiento del proceso interesa por las siguientes razones:

• Verificar que el proceso de encuentra en estado de control estadístico, es decir que la variabilidad observada no tiene un origen identificable y que sólo refleja la variabilidad natural del proceso.
• El rendimiento de un lote fuera de los límites de control es un indicio de la existencia de causas de fallos en el proceso de manufactura.
• Razones económicas.

Limitaciones de las cartas de control p y np

Las cartas de control p (proporción de unidades defectuosas) y np (número de unidades defectuosas) descritas habitualmente en los manuales de calidad estadística están basadas en la aproximación de la distribución binomial mediante la distribución normal (NIST). Esta aproximación se considera satisfactoria cuando \(n~p~(1-p)~>~5\) y \(0.1~<~p~<~0.9\) y no siempre es adecuada las siguientes razones:

• Si la fracción defectuosa es muy baja, el límite de control inferior puede ser negativo; de la misma forma, si el rendimiento es próximo al 100%, el límite de control superior puede ser mayor que 100%.
• Es frecuente observar que el rendimiento de los procesos de manufactura presentan distribuciones sesgadas (ver figura 2).

Cartas de control beta

Las cartas de control basadas en la distribución beta tienen por ventaja:

• La variable aleatoria con distribución beta, así como los limites de control de la carta, están definidos en el intervalo \(0 \le y \le 1\).
• La distribución beta refleja adecuadamente el sesgo observado en la distribución del rendimiento de procesos.

Objetivos

Implementación de R / RStudio la construcción de la carta beta para el control del rendimiento de un proceso de manufactura de comprimidos.

Librerías utilizadas

• Tidyverse: librerías dplyr y readr.
• fitdistrplus.

Datos

La tabla laboratory incluye dos campos relacionados con el rendimiento, tbl_yield (redimiento de la etapa de compresión) y batch_yield (reendimiento global, compresión + recubrimiento). En la tabla 7 de Zagar y Mihelic (2022) se recogen los rendimientos globales para la etapa de compresión (Tablet press yield (%), 98.3%) y para el producto final (Batch yield (%), 98.3%).

Importación de datos

Ver código para la importación de datos y los nombres de los campos de la tabla.

rm(list = ls(all = TRUE))
laboratory <- readr::read_delim(file = "./Datos/BigData/Laboratory.csv",
                              locale = locale(decimal_mark = "."),
                               col_types = list(
                                 batch = col_integer(),
                                 code = col_integer(),
                                 strength = col_factor(),
                                 size = col_integer(),
                                 start = col_character(),
                                 api_code = col_factor(),
                                 api_batch = col_factor(),
                                 smcc_batch = col_factor(),
                                 lactose_batch = col_factor(),
                                 starch_batch = col_factor(),
                                 api_water = col_double(),
                                 api_total_impurities = col_double(),
                                 api_l_impurity = col_double(),
                                 api_content = col_double(),
                                 api_ps01 = col_double(),
                                 api_ps05 = col_double(),
                                 api_ps09 = col_double(),
                                 lactose_water = col_double(),
                                 lactose_sieve0045 = col_double(),
                                 lactose_sieve015 = col_double(),
                                 lactose_sieve025 = col_double(),
                                 smcc_water = col_double(),
                                 smcc_td = col_double(),
                                 smcc_bd = col_double(),
                                 smcc_ps01 = col_double(),
                                 smcc_ps05 = col_double(),
                                 smcc_ps09 = col_double(),
                                 starch_ph = col_double(),
                                 starch_water = col_double(),
                                 tbl_min_thickness = col_double(),
                                 tbl_max_thickness = col_double(),
                                 fct_min_thickness = col_double(),
                                 fct_max_thickness = col_double(),
                                 tbl_min_weight = col_double(),
                                 tbl_max_weight = col_double(),
                                 tbl_rsd_weight = col_double(),
                                 fct_rsd_weight = col_double(),
                                 tbl_min_hardness = col_double(),
                                 tbl_max_hardness = col_double(),
                                 tbl_av_hardness = col_double(),
                                 fct_min_hardness = col_double(),
                                 fct_max_hardness = col_double(),
                                 fct_av_hardness = col_double(),
                                 tbl_max_diameter = col_double(),
                                 fct_max_diameter = col_double(),
                                 tbl_tensile = col_double(),
                                 fct_tensile = col_double(),
                                 tbl_yield = col_double(),
                                 batch_yield = col_double(),
                                 dissolution_av = col_double(),
                                 dissolution_min = col_double(),
                                 resodual_solvent = col_double(),
                                 impurities_total = col_double(),
                                 impurity_o = col_double(),
                                 impurity_l = col_double()
                               )
                               )

names(laboratory)

##  [1] "batch"                "code"                 "strength"            
##  [4] "size"                 "start"                "api_code"            
##  [7] "api_batch"            "smcc_batch"           "lactose_batch"       
## [10] "starch_batch"         "api_water"            "api_total_impurities"
## [13] "api_l_impurity"       "api_content"          "api_ps01"            
## [16] "api_ps05"             "api_ps09"             "lactose_water"       
## [19] "lactose_sieve0045"    "lactose_sieve015"     "lactose_sieve025"    
## [22] "smcc_water"           "smcc_td"              "smcc_bd"             
## [25] "smcc_ps01"            "smcc_ps05"            "smcc_ps09"           
## [28] "starch_ph"            "starch_water"         "tbl_min_thickness"   
## [31] "tbl_max_thickness"    "fct_min_thickness"    "fct_max_thickness"   
## [34] "tbl_min_weight"       "tbl_max_weight"       "tbl_rsd_weight"      
## [37] "fct_rsd_weight"       "tbl_min_hardness"     "tbl_max_hardness"    
## [40] "tbl_av_hardness"      "fct_min_hardness"     "fct_max_hardness"    
## [43] "fct_av_hardness"      "tbl_max_diameter"     "fct_max_diameter"    
## [46] "tbl_tensile"          "fct_tensile"          "tbl_yield"           
## [49] "batch_yield"          "dissolution_av"       "dissolution_min"     
## [52] "resodual_solvent"     "impurities_total"     "impurity_o"          
## [55] "impurity_l"

La tabla siguiente muestra el número de lotes en función del contenido de fármaco (campo strength) y del tamaño del lote (campo size).

T0 <- laboratory |> group_by(size, strength) |> tally() |>
  spread(strength, n)
print(T0)

## # A tibble: 9 x 5
## # Groups:   size [9]
##      size `5MG` `10M` `20M` `40M`
##     <int> <int> <int> <int> <int>
## 1  240000   218    NA    NA    NA
## 2  480000    15    NA    NA    NA
## 3  583000    NA    NA    NA    22
## 4  960000    NA    NA   437    NA
## 5 1100000    NA    NA    NA    76
## 6 1200000    NA    NA     8    NA
## 7 1920000    NA   211    NA    NA
## 8 2400000     3     2    NA    NA
## 9 4800000    13    NA    NA    NA

Rendimiento del proceso

La figura siguiente muestra los gráficos de cajas del rendimiento de la etapa de compresión en función del tamaño del lote (izquierda) y de la dosificación de fármaco (derecha). Zagar y Mihelic no informan sobre el valor mínimo requerido del rendimiento. En lo sucesivo supondremos (y por lo tanto arbitrariamente) que el valor mínimo requerido es 95% (línea roja de trazo discontínuo).

Observaciones:

• No se observa ninguna tendencia en función del tamaño del lote.
• Es apreciable un menor rendimiento de la manufactura de comprimidos de 5 mg de principio activo.

par(mfrow = c(1,2))
plot(tbl_yield ~ as.factor(size), data = laboratory, col = "lightblue", main = "compresión",
     ylim = c(85, 105), ylab = "rendimiento (%)", xlab = "tamaño lote (unds)")
abline(h = 95, lty = 2, col = "red")
plot(tbl_yield ~ strength, data = laboratory, col = "lightblue", main = "compresión",
     ylim = c(85, 105), ylab = "rendimiento (%)", xlab = "mg/comp")
abline(h = 95, lty = 2, col = "red")

Figura 1. Rendimientos en función del tamaño del lote y del contenido de fármaco

Para realizar un análisis más detallado del rendimiento de la etapa de compresión se seleccionaron lotes de igual tamaño (campo batch_size) y contenido de principio activo (campo strength). Ver la tabla siguiente.

Tabla 1. Casos seleccionados para el estudio del rendimiento del proceso de manufactura.

caso	strength	batch_size	número de lotes
A	5 mg	240,000	218
B	10 mg	1,920,000	211
C	20 mg	960,000	437
D	40 mg	1,100,000	76

 

Distribución beta

La aplicación de la distribución beta al control estadístico de la calidad ha sido revisada entre otros por Sant’Anna y ten Caten (2012). La figura 1 en el artículo de Hamasha y cols. (2025) muestra los perfiles que adopta la distribución beta en función de los valores de los parámetros parámetros de la distribución que reciben el nombre de parámetros de forma.

La función de densidad de una variable aleatoria \(Y\) con distribución beta es:

\[\begin{equation} f(y,~\theta_1,~\theta_2)~=~\frac{\Gamma(\theta_1 + \theta_2)}{\Gamma(\theta_1)~\Gamma(\theta_2)}{y^{\theta_1}~(1~-~y)^{\theta_2~-~1}} \tag{ec. 1} \end{equation}\]

donde \(y\) es la variable aleatoria (\(0~\le~y~\le~1\)), y \(\theta_1\) y \(\theta_2\) los parámetros de forma de la función (\(\theta_1~,~\theta_2 > 0\)). \(\Gamma(x)\) es la función gamma,

\[\begin{equation} \Gamma(x)~=~\int_0^\infty t^{x-1}~e^{-x}~dt \tag{ec. 2} \end{equation}\]

El valor esperado y la varianza de una variable aleatoria \(Y\) con distribución beta son:

\[\begin{equation} E(Y)~=~\frac{\theta_1}{\theta_1~+~\theta_2} \tag{ec. 3} \end{equation}\]

\[\begin{equation} Var(Y)~=~\frac{\theta_1~\theta_2}{(\theta_1+\theta_2)^2~(\theta_1+\theta_2+1)} \tag{ec. 4} \end{equation}\]

La ditribución beta está implementada en R. Para más información ejecutar ? Distributions en la consola de R.

Estimación de los factores de forma

Los métodos más comunes para la estimación de los parámetros de distribuciones estadísticas son:

• Método de los momentos, basado en la resolución de sistema formado por las ecuaciones 3 y 4 igualandolas a las estimadas muestrales de \(\bar{x}\) y \(\hat{\sigma}^2\). Este método se utiliza en la función bcc::calculate_limits() de la librería ccb
• Estimación por el método de maximización de la verisimilitud implementado en la función fitdist() de la librería fitdistrplus; esta librería también incluye también la funcion denscomp() para la representación gráfica de los resultados y pruebas para la comparación de distintas funciones de distribución.

Límites de control de la carta beta

Los límites de control inferior (\(LCL\)) y superior (\(UCL\)) se calculan a partir de la función de probabilidad de la distribución beta,

\[\begin{equation} F~=~\int_0^y~f(y,~\theta_1,~\theta_2)~dy \tag{ec. 5} \end{equation}\]

donde \(f(y,~\theta_1,~\theta_2)\) es la función de densidad (ecuación 1). Los límites de control son los cuantiles de la distribución beta calculados a partir de las estimadas de los factores de forma:

\[\begin{array}{lcl} LCL & = & F^{-1}(\alpha / 2,~\hat{\theta}_1,~\hat{\theta}_2) \\ UCL & = & F^{-1}(1-\alpha / 2,~\hat{\theta}_1,~\hat{\theta}_2) \tag{ec. 6} \end{array}\]

Para fijar los valores de \(\alpha/2\) y \(1-\alpha/2\), se utiliza el criterio habitual de las cartas de control basadas en la distribución normal, \(\pm 3 \sigma\), resultando 0.00135 y 0.99865 respectivamente.

En el código siguiente:

• Se extraen de la tabla laboratory el vector con los rendimientos de la etapa de compresión (campo tbl_yield) y se expresan como fracción.
• Se estiman los parámetros de la distribución beta para las cuatro formulaciones analizadas.
• Se representan los histogramas de los valores y la función de densidad estimada.

par(mfrow = c(2,2))

## cálculo fracción conforme para los casos seleccionados

yield_b05 <- laboratory |> dplyr::filter(size == 240000) |> dplyr::select(tbl_yield) |> mutate(fz = tbl_yield/100) |> dplyr::filter(fz <= 1)

yield_b10 <- laboratory |> dplyr::filter(size == 1920000) |> dplyr::select(tbl_yield) |> mutate(fz = tbl_yield/100) |> dplyr::filter(fz <= 1)

yield_b20 <- laboratory |> dplyr::filter(strength == "20M") |> dplyr::select(tbl_yield) |> mutate(fz = tbl_yield/100) |> dplyr::filter(fz <= 1)

yield_b40 <- laboratory |> dplyr::filter(size == 1100000) |> dplyr::select(tbl_yield) |> mutate(fz = tbl_yield/100) |> dplyr::filter(fz <= 1)

## estimación de los parámetros de la distribución beta

dist_b05f <- fitdist(yield_b05$fz, "beta", method = "mle")
dist_b10f <- fitdist(yield_b10$fz, "beta", method = "mle")
dist_b20f <- fitdist(yield_b20$fz, "beta", method = "mle")
dist_b40f <- fitdist(yield_b40$fz, "beta", method = "mle")

## gráficos

denscomp(dist_b05f, ylim = c(0, 75), xlim = c(0.90, 1), fitlw = 2,
         main = "5 MG", ylab = "densidad", xlab = "rendimiento")
denscomp(dist_b10f, ylim = c(0, 75), xlim = c(0.90, 1), fitlw = 2,
         main = "10 MG", ylab = "densidad", xlab = "rendimiento")
denscomp(dist_b20f, ylim = c(0, 75), xlim = c(0.90, 1), fitlw = 2,
         main = "20 MG", ylab = "densidad", xlab = "rendimiento")
denscomp(dist_b40f, ylim = c(0, 75), xlim = c(0.90, 1), fitlw = 2,
         main = "40 MG", ylab = "densidad", xlab = "rendimiento")

Figura 2. Histogramas y función de densidad

Parámetros estimados y límites de control

En el código siguiente:

• Se extraen las estimadas de \(\theta_1\) y \(\theta_2\) calculadas mediante la función fitdist().
• Se calculan los parámetros de la carta de control utilizando la ecuación 6.
• El gráfico siguiente muestra las cartas de control. Obsérvese la asimetría de los límites de control respecto del valor medio (líneas a trazos de color rojo).

Observaciones:

• La media de los rendimientos de los cuatro casos son muy próximo al valor informado por Zagar y Mihelic.
• El límite de control inferior es superior al 95% para las formulaciones de 10, 20 y 40 mg/comprimido, pero no para la formulación de 5 mg/comp. A falta de otra información, se deberia concluir la necesidad de revisar el proceso de manufactura de la formulación de 5 mg/comp.

k1 <- pnorm(-3)
k2 <- pnorm(3)
beta_par <- as_tibble(rbind(summary(dist_b05f)$estimate,
      summary(dist_b10f)$estimate,
      summary(dist_b20f)$estimate,
      summary(dist_b40f)$estimate))
colnames(beta_par) <- c("theta_1", "theta_2")
beta_par <- beta_par |> mutate(Ymean = theta_1/(theta_1+theta_2),
                               lcl = qbeta(k1, theta_1, theta_2),
                               ucl = qbeta(k2, theta_1, theta_2))
round(beta_par, 3)

## # A tibble: 4 x 5
##   theta_1 theta_2 Ymean   lcl   ucl
##     <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1    173.    4.19 0.976 0.928 0.997
## 2    335.    4.46 0.987 0.961 0.998
## 3    242.    3.67 0.985 0.952 0.998
## 4    252.    3.46 0.986 0.955 0.999

par(mfrow = c(2,2))

plot(yield_b05$tbl_yield, ylim = c(85, 105), main = "5 MG - 240 000 unds",
     cex = 0.5, col = "blue", ylab = "rendimiento", xlab = "índice lote")
abline(h = 100*beta_par[1, 3:5], lty = 2, col = "red")

plot(yield_b10$tbl_yield, ylim = c(85, 105), main = "10 MG - 1 920 000 unds",
     cex = 0.5, col = "blue", ylab = "rendimiento", xlab = "índice lote")
abline(h = 100*beta_par[2, 3:5], lty = 2, col = "red")

plot(yield_b20$tbl_yield, ylim = c(85, 105), main = "20 MG - 960 000 unds",
     cex = 0.5, col = "blue", ylab = "rendimiento", xlab = "índice lote")
abline(h = 100*beta_par[3, 3:5], lty = 2, col = "red")

plot(yield_b40$tbl_yield, ylim = c(85, 105), main = "40 MG - 1 100 000",
     cex = 0.5, col = "blue", ylab = "rendimiento", xlab = "índice lote")
abline(h = 100*beta_par[4, 3:5], lty = 2, col = "red")

La tabla siguiente resume algunos datos de los lotes con un rendimiento inferior al 95%.

batch_under_95 <- laboratory |> filter(tbl_yield < 95) |> dplyr::select(code, batch, strength, size, tbl_yield, code) |> arrange(strength)
batch_under_95

## # A tibble: 12 x 5
##     code batch strength    size tbl_yield
##    <int> <int> <fct>      <int>     <dbl>
##  1     1   124 5MG       240000      91.3
##  2     1   262 5MG       240000      94.3
##  3    11   332 5MG       480000      94.3
##  4     1   443 5MG       240000      92.8
##  5    21   529 5MG       240000      94.7
##  6    21   523 5MG       240000      94.3
##  7    21   557 5MG       240000      94.7
##  8    21   648 5MG       240000      91.3
##  9     7    82 20M      1200000      90.4
## 10    23   729 20M       960000      94.9
## 11    17   789 20M       960000      88.0
## 12    23   790 20M       960000      93.8

Bibliografía

Hamasha, M. M., Shawaheen, G., Mayyas, A. (2025) Exploring the shift in symmetry phenomenon in exponentially weighted moving average quality charts for statistics derived from beta distribution Statistics Optimization and Information Computing 14 1388 - 1403 https://doi.org/10.19139/soic-2310-5070-2128

Sant’Anna, A. M. O., ten Caten, C. S. (2012) Beta control chart for monitoring fraction data, Expert systems with Applications 39 10236 - 10243 https://doi.org/10.1016/j.eswa.2012.02.146

Zagar, J. and Mihelic, J. (2022) Big data collection in pharmaceutical manufacturing and its use for product quality predictions. Nature Scientific Data 9:99 https://doi.org/10.1038/s41597-022-01203-x