Contest

Una società brianzola lavora i propri prodotti tramite l’utilizzo di una serie di macchine a controllo numerico che necessitano di continua manutenzione, la quale comporta un costo sempre differente. Scopo dell’indagine è ricercare un pattern che preveda l’importo del costo nei prossimi mesi e pianificare di conseguenza le risorse.

library(caret)
library(ggplot2)
library(plotly)
library(gridExtra)
library(grid)
library(randomForest)
library(corrplot)

Exploratory Analysis

I dati sono contenuti in un file csv composto da 1600 osservazioni/macchine delle quali sono stati registrate 20 variabili, una delle quali è faultCost, il costo che si vuole prevedere. Alcune variabili sono di tipo numerico, altre invece sono di tipo fattoriale.

#data loading
dati <- read.csv('dati.csv')
str(dati)
## 'data.frame':    1600 obs. of  20 variables:
##  $ id                : int  968 240 819 692 420 1085 1998 365 1022 1240 ...
##  $ machineType       : Factor w/ 11 levels "DFTT","Enx","OPTX",..: 3 8 6 9 7 1 2 7 11 10 ...
##  $ startingTime      : int  220 420 120 120 220 320 240 240 100 220 ...
##  $ mediumRPM         : num  11962 11897 12043 12022 12040 ...
##  $ maxRPM            : int  17868 19775 19855 16600 18163 19665 19040 19150 19791 16939 ...
##  $ minRPM            : int  100 115 124 176 122 178 135 137 193 158 ...
##  $ rcol              : Factor w/ 20 levels "0A","0B","0C",..: 11 13 7 18 16 15 15 17 15 1 ...
##  $ wf                : int  0 1 0 0 0 1 0 1 1 1 ...
##  $ region            : Factor w/ 4 levels "East","North",..: 3 2 1 3 3 1 2 1 3 1 ...
##  $ margin            : num  21 27.9 16.8 25.7 13 ...
##  $ margin5           : num  1109 443 844 524 621 ...
##  $ margin10          : num  15643 7892 7536 10635 5863 ...
##  $ m1                : num  17861 8778 9224 11684 7105 ...
##  $ m2                : num  1151 499 877 576 647 ...
##  $ mstar             : num  1494 2731 2586 1617 999 ...
##  $ mediumRPMCorrected: num  116 121 121 117 114 ...
##  $ maxRPMCorrected   : num  139 150 150 135 138 ...
##  $ minRPMCorrected   : num  39.9 53.4 52.1 42.3 33.5 ...
##  $ faultTime         : int  1363 1689 1623 1035 1053 1287 1069 1418 289 806 ...
##  $ faultCost         : int  551 2433 1947 1245 574 2233 575 709 2411 1857 ...

Nel seguente grafico ad istogramma si osserva la distribuzione dell’importo di costo faultCost per ognuna delle 1600 macchine. Si può notare come il costo si concentri attorno ai valori 500 e 1500.

q1 <- qplot(dati$faultCost, main = 'Variabile di costo',
      ylab = 'N osservazioni',xlab = 'faultCost',
      fill=I('blue'),colour=I('white')) + 
    scale_x_continuous(breaks = seq(0,2500,500))
ggplotly(q1)

Fra le variabili registrate per ogni macchina ve ne sono 3 di tipo factor, machineType, rcol e region, la cui distribuzione è rappresentata nei seguenti grafici boxplot. Non si notano particolari differenze fra le varie classi e quindi non sembra esservi alcun pattern rispetto alla variabile di costo. In sostanza queste variabili non sembrano utili al fine di predire il costo di manutenzione macchine.

#type of machine
f1 <- ggplot(data = dati, aes(x=machineType,y=faultCost))+
    geom_boxplot(fill='lightblue') +
    theme(axis.text.x = element_text(angle = 90))
#rcol 
f2 <- ggplot(data = dati, aes(x=rcol, y=faultCost)) + 
    geom_boxplot(fill='yellow') +
    theme(axis.text.x = element_text(angle = 90))
#region
f3 <- ggplot(data = dati, aes(x=region, y=faultCost)) + 
    geom_boxplot(fill='green') +
    theme(axis.text.x = element_text(angle=45,hjust=1))
grid.arrange(f1,f2,f3,ncol=2,top='Variabili Factor')

Proseguendo con l’analisi, attraverso il calcolo della correlazione si può osservare la presenza di relazioni tra le variabili di tipo numerico. Nel seguente grafico, i coefficienti di correlazione più alti sono evidenziati da un colore più marcato, evidenziando quali variabili sono correlate tra loro. In particolare si osserva che la variabile faultCost è fortemente correlata con mediumRPMCorrected, mstar e minRPMCorrected. Queste tre variabili sono inoltre correlate fra loro.

num_columns <- setdiff(names(dati),c('id','machineType','rcol','region'))
cor_mat <- stats::cor(dati[,num_columns])
corrplot(cor_mat,method = 'number',type = 'lower',order = 'hclust',
         tl.cex = 0.7, diag = F,tl.srt = 45,tl.offset = 1,
         number.cex = 0.6, tl.col = 'black')

La forte correlazione registrata per mstar e mediumRPMCorrected può essere osservata nei seguenti grafici, dove le due variabili mostrano una positiva correlazione con faultCost. Ciò significa che esse potrebbero essere utilizzate per predire l’importo di costo delle macchine.

#high correlation variables
g1 <- ggplot(data = dati, aes(x=mediumRPMCorrected,y=faultCost)) + 
    geom_point() + geom_smooth(method = 'lm') +
    ylim(c(0,3000)) + xlab('') + ggtitle('mediumRPMCorrected')
g2 <- ggplot(data = dati, aes(x=mstar,y=faultCost)) + 
    geom_point() + geom_smooth(method = 'lm') +
    ylim(c(0,3000)) + xlab('') + ggtitle('mstar')
grid.arrange(g1,g2,ncol=2)

Modelli

La costruzione di alcuni modelli basilari può rivelare quali siano le variabili con il maggiore valore predittivo.

Un primo modello di regressione lineare, costruito su tutte le variabili numeriche dimostra per esempio che mstar e margin5 potrebbero essere utilizzate per costruire un modello.

#basic linear model
fit_lm <- lm(faultCost~., data = dati[,num_columns])
summary(fit_lm)  
## 
## Call:
## lm(formula = faultCost ~ ., data = dati[, num_columns])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -53.662 -24.024  -0.018  23.604  54.788 
## 
## Coefficients: (2 not defined because of singularities)
##                      Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)        -3.255e+03  2.713e+03   -1.200   0.2304    
## startingTime        1.639e-02  8.123e-03    2.018   0.0438 *  
## mediumRPM          -2.249e-01  1.972e-01   -1.141   0.2542    
## maxRPM             -8.786e-03  1.917e-02   -0.458   0.6469    
## minRPM              4.306e-02  3.802e-02    1.132   0.2576    
## wf                  2.368e+00  1.415e+00    1.673   0.0946 .  
## margin              1.055e+00  4.982e-02   21.184  < 2e-16 ***
## margin5            -1.002e+00  2.541e-03 -394.308  < 2e-16 ***
## margin10            9.495e-05  1.588e-04    0.598   0.5499    
## m1                         NA         NA       NA       NA    
## m2                         NA         NA       NA       NA    
## mstar               7.852e-01  1.711e-01    4.589 4.81e-06 ***
## mediumRPMCorrected  5.420e+01  4.650e+01    1.166   0.2440    
## maxRPMCorrected     2.676e+00  5.344e+00    0.501   0.6167    
## minRPMCorrected    -2.355e+00  2.513e+00   -0.937   0.3487    
## faultTime          -1.376e-03  1.717e-03   -0.801   0.4230    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.24 on 1586 degrees of freedom
## Multiple R-squared:  0.9985, Adjusted R-squared:  0.9985 
## F-statistic: 8.22e+04 on 13 and 1586 DF,  p-value: < 2.2e-16

Utilizzando quindi mstar e margin5 come variabili di input, si può ottenere un modello con R-squared pari a 0.99. Tale modello riesce a spiegare gran parte della varianza di faultCost e potrebbe essere usato per predirne il valore.

fit_lm2 <- lm(faultCost~mstar+margin5-1, data = dati)
summary(fit_lm2)
## 
## Call:
## lm(formula = faultCost ~ mstar + margin5 - 1, data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -234.06  -22.66   23.95   69.98  191.12 
## 
## Coefficients:
##          Estimate Std. Error t value Pr(>|t|)    
## mstar    1.060398   0.001671   634.5   <2e-16 ***
## margin5 -0.940084   0.005835  -161.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 70.64 on 1598 degrees of freedom
## Multiple R-squared:  0.9984, Adjusted R-squared:  0.9984 
## F-statistic: 4.897e+05 on 2 and 1598 DF,  p-value: < 2.2e-16

Passando ad un Random Forest, applicato nella sua versione più basilare, si può osservare quali siano le variabili che esso sceglie per costruire il modello ottimale. Il seguente grafico mostra le variabili ordinate per importanza e si nota come vi siano le variabili osservate precedentemente.

#random forest for feature importance
train <- createDataPartition(dati$faultCost,p=0.7,list = F)
fit_rf <- randomForest(faultCost~.,data = dati,subset = train,method='rf')
varImpPlot(fit_rf,n.var = 10, main = 'Variabili Importanti')

Per completare l’analisi si applica Recursive Feature Elimination, algoritmo di selezione automatica delle variabili, il cui risultato dimostra l’importanza delle variabili trovate nei modelli lineare e Random Forest.

#recursive feature elimination
control <- rfeControl(functions = rfFuncs, method = 'cv', number = 10)
results <- rfe(dati[,1:19],dati[,20],sizes = c(1:19),
               rfeControl = control)
results
## 
## Recursive feature selection
## 
## Outer resampling method: Cross-Validated (10 fold) 
## 
## Resampling performance over subset size:
## 
##  Variables   RMSE Rsquared RMSESD RsquaredSD Selected
##          1 849.17 0.002796  26.56   0.003289         
##          2 800.54 0.003355  22.29   0.004653         
##          3  74.25 0.991025  13.51   0.004418         
##          4  60.74 0.993104  13.58   0.003797        *
##          5  65.87 0.991841  13.84   0.004102         
##          6  63.26 0.992342  13.66   0.003884         
##          7  68.62 0.991088  14.15   0.004286         
##          8  73.78 0.989772  15.65   0.005072         
##          9  67.47 0.991436  12.81   0.003717         
##         10  70.97 0.990562  14.09   0.004341         
##         11  74.63 0.989676  13.88   0.004482         
##         12  69.28 0.991000  13.45   0.004058         
##         13  74.58 0.989624  17.37   0.005587         
##         14  82.08 0.988085  16.92   0.005868         
##         15  77.64 0.989369  15.54   0.004858         
##         16  82.15 0.988260  15.47   0.005038         
##         17  86.22 0.987199  15.61   0.005398         
##         18  84.05 0.988081  13.95   0.004845         
##         19  87.89 0.987125  14.11   0.005163         
## 
## The top 4 variables (out of 4):
##    margin5, m2, mstar, minRPMCorrected

Conclusioni

Riassumendo l’analisi svolta:

In conclusione, Exploratory analysis e modelli sembrano dimostrare che le variabili sopra menzionate potrebbero essere utilizzate per prevedere l’importo del costo delle macchine nei prossimi mesi e pianificare di conseguenza le risorse.