Descripción de Datos

La descripción gráfica de los datos nos permite inspeccionar los datos para:

detectar carencias (lectura de datos incorrecta/incompleta, datos faltantes, registros erróneos, …) y resolverlas
identificar el tipo de datos que se han leído y recodificar si es oportuno
inspeccionar la naturaleza de las variables
inspeccionar las posibles relaciones entre las variables.

Datasets

Tienes varias bases de datos con las que practicar en la web de OpenMV.net. En particular, trabajamos a continuación con la BD blender-efficiency.csv:

library(ggplot2)
blender=read.csv("blender-efficiency.csv",header=TRUE, sep=",")
summary(blender)

##   ParticleSize MixerDiameter MixerRotation  BlendingTime  
##  Min.   :2     Min.   :4     Min.   :175   Min.   :45.00  
##  1st Qu.:5     1st Qu.:6     1st Qu.:200   1st Qu.:48.75  
##  Median :5     Median :6     Median :200   Median :60.00  
##  Mean   :5     Mean   :6     Mean   :200   Mean   :58.33  
##  3rd Qu.:5     3rd Qu.:6     3rd Qu.:200   3rd Qu.:60.00  
##  Max.   :8     Max.   :8     Max.   :225   Max.   :75.00  
##  BlendingEfficiency
##  Min.   :61.80     
##  1st Qu.:70.75     
##  Median :74.70     
##  Mean   :75.17     
##  3rd Qu.:80.03     
##  Max.   :92.20

Reconvertimos en factores las variables que han de actuar como tales

blender$ParticleSize=as.factor(blender$ParticleSize)
blender$MixerDiameter=as.factor(blender$MixerDiameter)
blender$MixerRotation=as.factor(blender$MixerRotation)
blender$BlendingTime=as.factor(blender$BlendingTime)

Con el fin de no repetir continuamente el nombre completo de las variables, simplificaremos a continuación con las siguientes siglas:

BlendingEfficiency: BE
BlendingTime: BT
MixerRotation: MR
MixerDiameter: MD
ParticleSize: PS

Descriptivos univariantes

Revisamos en primer lugar los datos disponibles.
Identificamos la variable respuesta y posibles predictoras.
Graficamos individualmente la respuesta.

datos=blender
summary(datos)

##  ParticleSize MixerDiameter MixerRotation BlendingTime BlendingEfficiency
##  2: 4         4: 4          175: 4        45: 5        Min.   :61.80     
##  5:10         6:10          200:10        60:10        1st Qu.:70.75     
##  8: 4         8: 4          225: 4        75: 3        Median :74.70     
##                                                        Mean   :75.17     
##                                                        3rd Qu.:80.03     
##                                                        Max.   :92.20

ggplot(datos,aes(x=BlendingEfficiency))+geom_histogram(bins=6)

La variable respuesta es BT.
Todas las variables explicativas son factores, con 3 niveles cada uno.
No disponemos de datos para todos los cruces de niveles en los factores (requeriría \(3^4\) datos en total).
No podemos afirmar nada sobre la normalidad de la respuesta, puesto que ésta se exige para cada uno de los grupos en que se claasifica tras el ajuste de un modelo. Sin embargo, a priori se entrevée cierta simetría, que en principio no nos haría recelar de la posibilidad del ajuste de un modelo ANOVA.

Descriptivos bivariantes: relaciones con la respuesta

Investigamos a continuación la relación entre las variables potencialmente predictoras y la variable respuesta. ### BlendingTime

ggplot(datos,aes(BlendingTime,BlendingEfficiency))+geom_boxplot()

Hay más variabilidad en los niveles 45 y 75 que en el nivel 60, si bien en éste todavía hay alunos datos (muy pocos) muy distanciados del resto (caja), tanto por arriba como por abajo (valores altos y valores bajos de BE).
La tendencia que se muestra en la relación con BE es al alza: niveles más altos de BT dan valores más altos de BE.
No parece que vayan a poder identificarse diferencias entre los niveles 45 y 60 de BT respecto de cómo influyen en la eficiencia BE (la caja 60 está encajada en el rango de variación de 45), aunque las medianas difieren.
Es algo más dudoso que puedan detectarse diferencias entre los niveles 60 y 75 de BT respecto de cómo influyen en la eficiencia BE, puesto que a pesar de que la caja de 60 está encajada en el rango de variación de 75, las cajas (el 75% de los datos entre los cuartiles 1 y 3) están prácticamente desencajadas.
Sí parece que puedan encontrarse diferencias entre los niveles 45 y 75 a la hora de provocar diferente eficiencia BE, puesto que las cajas no están demasiado encajadas entre sí y las medianas difieren considerablemente (la mediana en BT=75 supera el percentil 75 en el nivel BT=45).
Al haberse apreciado alguna diferencia entre algunos niveles, es posible que esta variable resulte significativa a la hora de explicar la eficiencia BE, esto es, que el hecho de utilizar uno u otro valor de BT implique variaciones en la eficiencia BE.

MixerRotation

ggplot(datos,aes(MixerRotation,BlendingEfficiency))+geom_boxplot()

Las medianas de eficiencia BE son prácticamente similares para los diferentes niveles de MR.
El rango de variación de la eficiencia BE para los niveles MR=175 y MR=225 está contenido en el rango de variación para el nivel MR=200.
La variabilidad de BE para MR=200 es considerablemente mayor que para los otros niveles. MR=200 es el nivel que mayor variabilidad genera en la respuesta.
La variabilidad de BE para MR=175 y MR=225 es considerablemente pequeña.
No parece que vaya a poderse demostrar una afectación estadística de MR sobre la eficiencia BE.

MixerDiameter

ggplot(datos,aes(MixerDiameter,BlendingEfficiency))+geom_boxplot()

A través de las medianas se aprecia una leve tendencia al alza en el efecto que tiene MD sobre la eficiencia BE: a mayor valor de MD, mayor la eficiencia ME.
Sin embargo, el rango de variación en los tres niveles de MD es prácticamente similar (ver rango de las cajas), lo que a priori no da indicios remarcables de que MD pueda provocar diferencias significativas en la eficiencia BE.
No parece que vaya a poderse demostrar una afectación estadística de MD sobre la eficiencia BE.

ParticleSize

ggplot(datos,aes(ParticleSize,BlendingEfficiency))+geom_boxplot()

Hay diferencias considerables en la eficiencia BE para los tres niveles de PS: tanto en las medianas como en el rango de variación (cajas desencajadas).
La tendencia observada es a la baja: mayor valor de PS provoca menor eficiencia BE.
La variabilidad en BE es algo mayor para PS=2 que para el resto, siendo mínima (salvo un valor) para PS=8.
Prácticamente sin ninguna duda PS será una variable significativa para explicar la eficiencia BE.

Descriptivos trivariantes: interacciones entre los factores

El objetivo básico de estos gráficos es indagar si existen interacciones entre los factores explicativos, esto es, si el hecho de que un factor varíe en un sentido afecta a la relación entre otro factor y la variable respuesta. Recordar que tendencias similares en la relación entre un factor A y la respuesta R para los diferentes niveles de un factor B, implicará que NO hay interacción. En caso contrario tendremos indicios de interacción.

Interacción BT con MR

ggplot(datos,aes(BlendingTime,BlendingEfficiency))+geom_boxplot()+facet_grid(.~MixerRotation)

Observando la relación BT-BE para los niveles MR=175 y MR=200 (en particular las medianas y cajas), apreciamos el mismo comportamiento creciente: a mayor BT, mayor BE.
La relación para el nivel MR=225 no está clara, pues si bien parece seguir un comportamiento contrario (a mayora BT, menor BE), faltan datos para BT=75.
La variabilidad es mucho mayor para MR=200 que para el resto de niveles, como ya se comentó en el gráfico de BE versus MR.
No se aprecia una interacción clara entre BT y MR a la hora de explicar BE.

Interacción BT con MD

ggplot(datos,aes(BlendingTime,BlendingEfficiency))+geom_boxplot()+facet_grid(.~MixerDiameter)

Faltan datos, de modo que no es posible identificar tendencias similares ni distintas a la hora de explicar BE con BT variando los niveles de MD.
No es posible discernir interacción alguna entre MD y BT para explicar BE.

Interacción BT con PS

ggplot(datos,aes(BlendingTime,BlendingEfficiency))+geom_boxplot()+facet_grid(.~ParticleSize)

La tendencia que se aprecia en la relación BT-BE es similar para los distintos niveles de PS: a mayor BT, mayor BE.
En principio no parece que vaya a existir interacción alguna entre BT y PS para explicar la eficiencia BE.

Tú mismo/a puedes realizar el análisis gráfico pertinente y la interpretación para el resto de interacciones.

Conclusiones del análisis descriptivo

Una vez realizada la inspección gráfica, las conclusiones que derivamos y que nos orientan a la hora de plantear el análisis estadístico son las siguientes:

Parece que BT afecta a la eficiencia BE.
Parece que MR no afecta a la eficiencia BE.
Parece que MD no afecta a la eficiencia BE.
PS parece afectar de forma considerable a la eficiencia BE.
No se aprecia una interacción clara entre BT y MR a la hora de explicar BE, si bien en MR=200 la relación BT-BE parece invertirse respecto de la observada para los otros niveles de MR.
No es posible discernir interacción alguna entre MD y BT para explicar BE.
No parece que vaya a existir interacción alguna entre BT y PS para explicar la eficiencia BE.

Análisis estadístico: modelización

Planteamos pues una modelización lineal y procedemos con la selección de variables, que puede responder a dos criterios básicos:

la significatividad estadística: ANOVA
el AIC: un modelo es mejor si es menor su AIC (combinación de verosimilitud y simplicidad del modelo para explicar los datos).

Podríamos utilizar las conclusiones derivadas del análisis descriptivo para plantear un modelo base de partida en la búsqueda de una buena modelización de la respuesta. Sin embargo, para evitar sesgos en la modelización provocados por el modelizador, es habitual tomar como partida un modelo completo con todas las interacciones.

Recordar que el modelo sólo será capaz de estimar efectos (medias) en los niveles y combinaciones de niveles para los que existan datos. Muchas interacciones serán inestimables a priori por la ausencia de datos en esta base de datos.

fit.b=lm(BlendingEfficiency~(BlendingTime*MixerRotation*MixerDiameter*ParticleSize),data=blender)
anova(fit.b)

## Analysis of Variance Table
## 
## Response: BlendingEfficiency
##                             Df  Sum Sq Mean Sq F value  Pr(>F)  
## BlendingTime                 2  135.17   67.59  2.3984 0.29425  
## MixerRotation                2    3.37    1.68  0.0597 0.94364  
## MixerDiameter                2   22.78   11.39  0.4043 0.71212  
## ParticleSize                 1 1026.05 1026.05 36.4104 0.02638 *
## BlendingTime:MixerRotation   3   14.30    4.77  0.1692 0.90893  
## MixerRotation:MixerDiameter  2    0.28    0.14  0.0050 0.99503  
## BlendingTime:ParticleSize    2    7.78    3.89  0.1381 0.87867  
## MixerDiameter:ParticleSize   1    8.70    8.70  0.3088 0.63427  
## Residuals                    2   56.36   28.18                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comprueba tú mismo qué efectos se han conseguido estimar (y cuáles no) ejecutando el código:

summary(fit.b)

Selección automática AIC

Buscamos a continuación el mejor modelo por selección automática basado en el AIC

fit.aic=step(fit.b)

anova(fit.aic)

## Analysis of Variance Table
## 
## Response: BlendingEfficiency
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## BlendingTime   2  135.17   67.59  8.2829  0.006392 ** 
## MixerDiameter  2   22.84   11.42  1.3994  0.287437    
## ParticleSize   2 1027.03  513.51 62.9313 9.508e-07 ***
## Residuals     11   89.76    8.16                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El modelo ajustado por AIC es el que contiene como factores predictores: BT, MD y PS. MR ha sido descartado. Notar que sin embargo, MD no es significativo (p-valor=0.2874).

Selección manual por significatividad

Dados los resultados del ANOVA en el modelo completo, reajustamos la modelización a un modelo base SIN interacciones. A partir de él vamos descartando variables no significativas.

fit=lm(BlendingEfficiency~(BlendingTime+MixerRotation+MixerDiameter+ParticleSize),data=blender)
anova(fit)

## Analysis of Variance Table
## 
## Response: BlendingEfficiency
##               Df  Sum Sq Mean Sq  F value    Pr(>F)    
## BlendingTime   2  135.17   67.59   7.7305  0.009346 ** 
## MixerRotation  2    3.37    1.68   0.1925  0.827869    
## MixerDiameter  2   22.78   11.39   1.3030  0.314131    
## ParticleSize   1 1026.05 1026.05 117.3569 7.598e-07 ***
## Residuals     10   87.43    8.74                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La primera candidata a salir del modelo es MR.

fit.s=update(fit,.~.-MixerRotation,data=blender)
anova(fit.s)

## Analysis of Variance Table
## 
## Response: BlendingEfficiency
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## BlendingTime   2  135.17   67.59  8.2829  0.006392 ** 
## MixerDiameter  2   22.84   11.42  1.3994  0.287437    
## ParticleSize   2 1027.03  513.51 62.9313 9.508e-07 ***
## Residuals     11   89.76    8.16                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La siguiente candidata a salir es MD.

fit.s=update(fit.s,.~.-MixerDiameter,data=blender)
anova(fit.s)

## Analysis of Variance Table
## 
## Response: BlendingEfficiency
##              Df  Sum Sq Mean Sq F value    Pr(>F)    
## BlendingTime  2  135.17   67.59  7.7879  0.005979 ** 
## ParticleSize  2 1026.80  513.40 59.1579 2.962e-07 ***
## Residuals    13  112.82    8.68                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Este modelo ya es aceptable. Lo comparamos con el que obtuvimos con el criterio AIC:

AIC(fit.s)

## [1] 96.11948

summary(fit.s)

## 
## Call:
## lm(formula = BlendingEfficiency ~ BlendingTime + ParticleSize, 
##     data = blender)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3312 -1.7039 -0.3219  0.5687  7.2688 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      84.278      1.885  44.713 1.27e-15 ***
## BlendingTime60    1.950      1.614   1.209  0.24837    
## BlendingTime75    8.187      2.185   3.748  0.00244 ** 
## ParticleSize5   -11.747      1.766  -6.652 1.58e-05 ***
## ParticleSize8   -22.650      2.083 -10.873 6.75e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.946 on 13 degrees of freedom
## Multiple R-squared:  0.9115, Adjusted R-squared:  0.8843 
## F-statistic: 33.47 on 4 and 13 DF,  p-value: 9.898e-07

AIC(fit.aic)

## [1] 96.00342

summary(fit.aic)

## 
## Call:
## lm(formula = BlendingEfficiency ~ BlendingTime + MixerDiameter + 
##     ParticleSize, data = blender)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2983 -1.0541 -0.2830  0.3391  7.3017 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      82.200      2.965  27.721 1.57e-11 ***
## BlendingTime60    2.302      2.460   0.936  0.36946    
## BlendingTime75    8.166      2.122   3.848  0.00271 ** 
## MixerDiameter6    2.127      2.578   0.825  0.42688    
## MixerDiameter8    3.375      2.020   1.671  0.12292    
## ParticleSize5   -11.829      1.769  -6.686 3.44e-05 ***
## ParticleSize8   -22.650      2.020 -11.213 2.33e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.857 on 11 degrees of freedom
## Multiple R-squared:  0.9296, Adjusted R-squared:  0.8912 
## F-statistic:  24.2 on 6 and 11 DF,  p-value: 9.9e-06

A pesar de tener una variable más que no es significativa, MD, resulta más eficiente el modelo seleccionado por el AIC. Es lógico que al incluir una variable más, se reduzcan más los residuos y se explique mejor la variabilidad en la respuesta: ver \(R^2\), Residual standard error, rango de variación de los residuos.

Hacemos a continuación una revisión rápida sobre la validación del modelo elegido (AIC):

par(mfrow=c(2,2))
plot(fit.aic)

Hay básicamente tres observaciones, la 3, 6 y especialmente la 8, que provocan residuos especialmente altosy desviaciones de normalidad. Cabría considerar la inspección de estos valores.
En principio damos por válido el modelo: no hay evidencias claras de tendencias sin explicar, cambios de variabilidad, desviaciones de normalidad en el global de los datos, …

Conclusiones

Se ratifican prácticamente todas las conclusiones del análisis descriptivo preliminar:

efecto de BT, MD y PS sobre la eficiencia BE
no es posible inferir sobre interacciones entre los factores
hay diferencias significativas entre BT=75 Y BT=45 (p-valor 0.00271 en la tabla summary) a la hora de provocar diferente eficiencia
hay diferencias significativas entre cualquiera de los dos niveles PT=5 y PT=8 respecto del nivel base PT=2 (p-valores \(3.44e-05\) y \(2.33e-07\) respectivamente en la tabla summary) a la hora de provocar diferente eficiencia.

Predicción

Predecimos a continuación la eficiencia BE para cada una de las combinaciones de los factores implicados en el ajuste:

newdata=expand.grid(BlendingTime=factor(levels(datos$BlendingTime)), MixerDiameter=factor(levels(datos$MixerDiameter)), ParticleSize=factor(levels(datos$ParticleSize)))
pred=predict(fit.aic,newdata=newdata,se.fit=T,interval="prediction")
pred.full=cbind(newdata,round(pred$fit,2),round(pred$se.fit,2))
library(DT)
datatable(pred.full,rownames=FALSE,colnames=c("BT","MD","PT", "pred","ic.lwr", "ic.upr", "se"), filter="top",options=list(searching=TRUE,paging=TRUE,info=FALSE,pageLength=9,autoWidth=TRUE),
caption="Tabla Predicción. Predicción de la eficiencia BE para los niveles y factores considerados como relevantes (junto con el error de la predicción -se- y un intervalo de confianza al 95%).")

A partir de los resultados en la predicción, concluye tú mismo/a sobre cuál sería la configuración óptima en los factores explicativos para provocar máxima eficiencia.

Referencias

Datos obtenidos de https://openmv.net/
RStudio Team (2016). RStudio: Integrated Development for R. RStudio, Inc., Boston, MA URL http://www.rstudio.com/.
Librería ggplot2. http://ggplot2.tidyverse.org/reference/index.html
Librería DT. An interface to the DataTables library https://rstudio.github.io/DT/

Descripción de Datos

Mejora de Procesos 1479 (@asunmayoral)