2.1 Descripción del experimento

^{Triana et al. (2009)} usaron un diseño factorial 2x3² para medir el porcentaje de grado de secado (GS) del acetato de calcio dihidratado mediante un secador de tipo transportador en función de las covariables: velocidad rotacional del tornillo (N) en dos niveles 3 y 7 rpm, la temperatura (T) del sólido en tres niveles 45, 53, 61 grados centígrados y el flujo de alimentación (C) en tres niveles de 6, 8 y 10 L/min. En la Figura 1 se muestra un esquema del secador tipo tornillo transportador usado en el experimento.

Fuente: ^{Triana et al. (2009)}.

Figura 1 Ilustración sistema de secado usado en el experimento. 

Los elementos que conforman el sistema de secado fueron: un cristalizador con chaqueta C-101, un tornillo secador S-101, una bomba de doble diafragma B-101 y válvulas de recirculación y bloqueo (V-102, V-101 y V-103). La presión del aire de entrada se mantuvo constante en 40 psi y el SCD usado tenía 8.66 metros de longitud y 0.17 m de diámetro del tornillo. La temperatura [°C] del tornillo secador se controló por medio de la válvula V-105. La variable velocidad del tornillo [rpm] se controló por medio del variador de frecuencia VG-101 y el flujo de alimentación [L/min] se controló con la válvula V-104.

En la Tabla 1 se presentan los niveles de las variables de estudio (T, N y C) y la variable respuesta (GS). ^{Triana et al. (2009)} decidieron eliminar las observaciones 1 y 13 destacadas en negrita porque estos datos causan una desviación adimensional que se atribuye a una propiedad del sólido no controlada. Al igual que Triana et al. (2009), aquí se usaron solo las observaciones restantes para construir los modelos de regresión.

2.2 GAMLSS (Generalized Additive Model for Location Scale and Shape)

Los modelos de regresión son una de las herramientas más populares en estadística y sus orígenes se encuentran a finales del siglo XVIII. ^{Legendre (1805)} reportó el método de mínimos cuadrados y lo usó para analizar un conjunto de datos provenientes de una encuesta en la cual se tenían dos variables y cinco observaciones. ^{Gauss (1809)} publicó una justificación del método de mínimos cuadrados e introdujo la distribución normal para los errores del modelo.

En los modelos de regresión el objetivo es explicar la media (o valor esperado) de una variable y en función de un conjunto de p variables explicativas x ₁,⋯,x _p de la siguiente manera.

Los supuestos de este modelo de regresión son:

- se disponen de la información de p variables explicativas y una variable respuesta y para n unidades muestrales o casos,

- las variables observadas para un caso son independientes de los otros casos,

- la variable respuesta y debe tener distribución normal N(µ,σ ²) y

- la varianza σ² debe ser constante.

En la práctica los usuarios se encuentran con variables respuesta que no necesariamente tienen distribución normal, como el número de defectos en un sistema productivo (variable de conteo), el contenido neto de una bebida gaseosa (variable positiva) o el porcentaje de un medicamento que queda en el cuerpo de algún paciente luego de cierto periodo (variable entre 0 y 1). En estos casos sería un error usar un modelo de regresión clásico ya que se viola el supuesto de normalidad y por lo tanto se hace necesario constar una alternativa de modelación.

Los modelos GAMLSS (generalized additive models for location, shape and scale) fueron propuestos por ^{Rigby y Stasinopoulos (2005)} surgieron como una extensión a los modelos de regresión tradicionales porque flexibilizan el supuesto de normalidad de la variable respuesta y permiten que el usuario pueda elegir la distribución entre un listado de más de cien distribuciones. Otra ventaja de los modelos GAMLSS es que se pueden modelar todos los parámetros de la distribución seleccionada y no solo la media como en los modelos de regresión tradicionales.

Los modelos GAMLSS asumen que se dispone de n observaciones independientes y con función de densidad (o de masa) f ( y|θ ) en la cual el vector θ = (μ,σ,η,τ)^T representa el vector parámetros que se pueden modelar en función de las variables explicativas. Los primeros dos parámetros µ y σ corresponden por lo general a los parámetros de ubicación y escala de la distribución, mientras que los parámetros restantes (si aplican a la distribución) están asociados a otras características de la distribución.

Simbólicamente los modelos GAMLSS se puede escribir de la siguiente manera.

donde las funciones g _k (∙) con k = 1,2,3,4 se conocen como funciones de enlace y sirven para transformar el lado derecho (predictor lineal) de cada ecuación anterior en valores apropiados para cada parámetro. Los parámetros β _kj con j= 0,1,⋯,p corresponden a los efectos de las variables explicativas y se estiman por el método de máxima verosimilitud. Los modelos GAMLSS son muy flexibles porque en el lado derecho (predictor lineal) de las cuatro expresiones anteriores, se pueden incluir polinomios, polinomios fraccionales, splines, efectos aleatorios, entre otros.

En el lenguaje de programación R (2020) está disponible el paquete gamlss con el cual los usuarios pueden aplicar modelos GAMLSS de una forma sencilla. Este paquete cuenta con más de 100 distribuciones continuas, discretas y mixtas para seleccionar la distribución de la variable respuesta y es posible que el usuario defina su propia función de densidad, para una mayor compresión de estos modelos se recomienda consultar ^{Stasinopoulos et al. (2017)}.

2.3 Selección de la distribución para la variable respuesta

Para identificar las distribuciones de probabilidad (con soporte en R, R⁺ y (0,1)) que mejor explican el comportamiento de la variable porcentaje grado de secado, se usó un enfoque marginal, es decir, sin utilizar covariables auxiliares. Para esto se usó la función fitDist del paquete gamlss la cual suministra una lista de distribuciones que se ajustan mejor a la variable respuesta. El criterio de selección interno de fitDist es el GAIC (criterio de información de ^{Akaike generalizado) propuesto por Akaike (1974)}.

En la Figura 2 se presenta seis veces el histograma de la variable respuesta grado de secado (GS) y las curvas de densidad empíricas para las distribuciones normal, exponencial, Weibull 3, gamma, beta y logit-normal. Los valores de GAIC en esta búsqueda marginal para las distribuciones normal, exponencial, Weibull 3, gamma, beta y logit-normal fueron -55,33; -61,32; -60,96; -60,71; -60,95 y -58,28 respectivamente. Aunque la distribución normal tiene soporte en los números reales y la variable respuesta es positiva, aquí se incluyó la distribución normal como distribución de referencia ya que fue la usada por ^{Triana et al. (2009)} para analizar los datos experimentales.

Figura 2 Histograma para la variable grado de secado con las curvas de densidad empíricas para las distribuciones normal, exponencial, Weibull 3, gamma, beta y logit-normal. 

De la figura anterior se observa que las seis densidades (distribuciones) intentan reproducir la forma del histograma. De las seis distribuciones consideradas, la exponencial es la única con un parámetro, mientras que las restantes tienen dos parámetros. Estas seis distribuciones identificadas aquí fueron las distribuciones usadas para construir los modelos de regresión que se presentan en la siguiente sección.

3.1 Aplicación de los modelos de regresión

^{Triana et al. (2009)} ajustaron un modelo de regresión lineal múltiple y usaron las variables velocidad rotacional del tornillo (N), temperatura (T) y flujo de alimentación (C) para explicar el grado de secado (GS). El modelo ajustado que reportaron Triana et al. (2009) se muestra a continuación:

Este modelo se caracteriza por tener un pseudo R ² = 0,86 con una correlación lineal de Pearson entre las observaciones y _i . y las estimaciones Ê ( y _i) de 0,93. Este modelo se usará como modelo de referencia para comparar los modelos propuestos.

En el presente artículo se consideraron modelos GAMLSS con variable respuesta distribuida normal, exponencial, Weibull 3, gamma, beta y logit-normal. La distribución normal tiene soporte en R, la distribución exponencial, Weibull 3 y gamma tienen soporte en R⁺ mientras que las distribuciones beta y logit-normal tienen soporte en (0; 1).

Para cada uno de estos casos se modeló el parámetro μ de cada distribución en función de las variables velocidad rotacional del tornillo (N), temperatura (T) y flujo de alimentación (C), usando términos cuadráticos e interacciones entre estas variables. Para las distribuciones normal, Weibull 3, gamma, beta y logit-normal que tienen un segundo parámetro σ , no se utilizaron covariables en la modelación de este parámetro debido a que sólo se cuentan con 16 observaciones experimentales lo cual dificulta la estimación de efectos.

En la Tabla 2 se presenta un resumen de los modelos ajustados. En la segunda fila de la Tabla 2 está el modelo de referencia propuesto por ^{Triana et al. (2009)} y en las filas siguientes están los modelos finales luego de un proceso de selección de variables. En las tres últimas columnas de la Tabla 2 se presentan tres medidas para evaluar y comparar los modelos. La columna de correlación corresponde al coeficiente de correlación lineal de Pearson entre las observaciones y _i . y las estimaciones Ê ( y _i) dadas por el modelo ajustado.

Tabla 2 Medidas de comparación para los modelos ajustados. 

Modelo	Distribución	Variable respuesta	Correlación	GAIC	Pseudo R 2
Referencia	Normal	GS	0,93	-79,03	0,86
Mod1	Normal		0,92	-57,14	0,90
Mod2	Exponencial	GS	0,73	-66,65	0,41
Mod3	Gamma	GS	0,96	-82,84	0,85
Mod5	Weibull 3	GS	0,96	-84,06	0,86
Mod6	Beta	GS	0,98	-94,86	0,93
Mod7	Logit-normal	GS	0,83	-75,16	0,73

Una característica deseable en un modelo es que genere estimaciones Ê ( y _i) cercanas a los valores reales (y _i) de la variable de interés. De la tabla anterior se observa que los modelos gamma, Weibull 3 y beta presentaron los mayores valores de correlación y en particular el modelo beta presentó la más alta correlación con un valor de 0.98.

En la Figura 3 se muestra el diagrama de dispersión entre y ₁ . y Ê ( y _i) para los modelos gamma, Weibull 3 y beta. La línea a 45° que se observa en los diagramas de dispersión sirve como referencia para saber si las estimaciones son cercanas a los valores observados. De la figura se observa que las estimaciones Ê ( y _i) acompañan a los valores reales y _i del grado de secado observados en el experimento.

Figura 3 Diagrama de dispersión entre y. y las estimaciones Ê (y _¡) para los modelos gamma, Weibull 3 y beta. La recta en color azul a trazos es una recta de referencia a 45°. 

Los residuales RQR (randomized quantile residuais) son útiles para evaluar modelos cuya variable respuesta es diferente a la normal y éstos residuales tienen distribución normal si el modelo fue bien ajustado (^{Stasinopoulos et al., 2017}). En la Figura 4 se muestran los gráficos cuantil-cuantil para los residuales RQR de los modelos gamma, Weibull 3 y beta. De la figura se observa que los residuales están, en la mayoría de los casos, dentro de las bandas de confianza de color azul lo que indica que los residuales tiene distribución normal estándar. Adicionalmente, se aplicó la prueba Kolmogorov-Smirnov para estudiar la normalidad de los residuales y se encontraron valores-P de 0,88; 0,50 y 0,90 para los residuales de los modelos gamma, Weibull 3 y beta respectivamente, lo cual indica que podemos asumir que los residuales tienen distribución normal estándar.

Figura 4 Gráfico cuantil-cuantil para los residuales RQR modelos gamma, Weibull 3 y beta. 

Al comparar los modelos según su GAIC, correlación, Pseudo R² y normalidad de los residuales, se decidió elegir el modelo beta como modelo final porque éste presentó los mejores resultados. En la siguiente tabla se presenta el resumen de los efectos estimados para el modelo beta.

Usando la información de la tabla anterior es posible escribir las expresiones matemáticas para los parámetros µ^ y σ^ del modelo de regresión beta. En el modelo de regresión beta ambos parámetros µ^ y σ^ deben estar en el intervalo (0,1), es por esa razón es que se usa la función de enlace logit en las siguientes expresiones.

Con la parametrización de la distribución beta usada en este artículo el valor esperado y la varianza están dadas por E (X) = µ y Var (X) = σ ² µ (1 - µ) respectivamente. Esto significa que el parámetro µ coincide con el valor esperado de la variable y por lo tanto el valor estimado (promedio) de la variable GS se puede escribir así:

De igual forma la varianza estimada para la variable GS se puede escribir en forma compacta (en términos de µ^ y σ^ así:

Usando la ecuación de Ê (GS) mostrada arriba, se aplicó un método de búsqueda para obtener los valores de T, N y C que maximizan la variable grado de secado GS y lo que se encontró fue que el valor máximo grado de secado es 35,46% (o 0,3546 en fracción) cuando T =45, N =3 y C =6.

Como la expresión para Ê (GS) involucra tres variables, la variable temperatura (T) se fijó en sus cuantiles Q ₁ = 45, Q ₂ = 53 y Q ₃ = 61para poder generar los gráficos de contornos para Ê (GS) mostrados en la Figura 5.En el panel izquierdo de la figura se observa un punto de color azul que corresponde a la combinación T =45, N =3 y C= 6, la cual maximiza el grado de secado a un valor de 35,46% (o 0,3546 en fracción).

Figura 5 Gráfico de contornos para Ê (GS). En la izquierda el gráfico de contornos para T=45, en la mitad el gráfico de contornos para T=53 y a la derecha el gráfico de contornos para T=61. El punto de color azul en el gráfico de la izquierda corresponde al valor máximo de Ê (GS)=0,3546 que se obtiene con la combinación T=45, N=3 y C= 6. 

Usando la ecuación de mostrada anteriormente, se aplicó un método de búsqueda para obtener los valores de T, N y C que minimizan la varianza de variable grado de secado GS y lo que se encontró fue que la menor variabilidad del grado de secado es 1,264844e-05 y se obtiene cuando T =51,6; N =7 y C =8.

En la Figura 6 se muestran los gráficos de contornos para cuando la variable T =45; 51,6; 61. En el panel central hay un punto de color azul para identificar la combinación T =51,6; N =7 y C =8 que minimiza la varianza estimada de GS.

Figura 6 Gráfico de contornos para . En la izquierda el gráfico de contornos para T=45, en la mitad el gráfico de contornos para T=51,6 y a la derecha el gráfico de contornos para T=61. El punto de color azul en el gráfico de la mitad corresponde al valor mínimo de =1,264844e-05 que se obtiene con la combinación T=51,6; N=7 y C=8.