DATA PREPARATION

Load required package:

LOad Dataset & Data Structur

drilling <- read.csv("Loris_Drilling.csv", stringsAsFactors = T)

str(drilling)
## 'data.frame':    54196 obs. of  27 variables:
##  $ wellName   : Factor w/ 7 levels "Loris 10H L1",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ datasetName: Factor w/ 3 levels "Drilling","Drilling L1",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Depth      : num  401 402 402 402 403 ...
##  $ Bitm       : num  0.64 1.16 1.63 2.11 2.67 3.16 3.6 4.11 4.63 5.07 ...
##  $ BitRevs    : int  280 305 326 347 371 392 411 432 453 471 ...
##  $ DrillHrs   : num  0.01 0.02 0.02 0.03 0.03 0.04 0.04 0.05 0.06 0.06 ...
##  $ Dxc        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ FlowRate   : num  414 415 414 415 415 ...
##  $ MSE        : num  439 419 407 437 471 ...
##  $ MSE_DAVG   : num  1.03 1.02 1.02 1.02 1.02 ...
##  $ MW.In      : num  0.216 0.216 0.216 0.216 0.216 0.216 0.216 0.216 0.216 0.216 ...
##  $ MW.Out     : num  0.26 0.27 0.28 0.28 0.29 0.29 0.3 0.3 0.31 0.32 ...
##  $ POWER      : num  2100 2094 2090 2099 2082 ...
##  $ Pump.time  : num  75.9 78.4 80.5 81.4 84.6 ...
##  $ PVT        : num  1.99 2 2 2.01 2.02 2.02 2.03 2.03 2.04 2.04 ...
##  $ ROP        : num  61 60 60 61 61 60 60 60 61 61 ...
##  $ RotatingHrs: num  61 60 60 61 61 60 60 60 61 61 ...
##  $ RPM        : num  226 227 226 227 227 ...
##  $ RPM.SURF   : num  32.9 32.9 33 33 32.9 32.9 33 33 33 32.9 ...
##  $ SPP        : num  30.8 30.8 30.8 30.8 30.7 30.8 30.8 30.8 30.8 30.8 ...
##  $ Temp.In    : num  2 3 3 2 2 3 3 3 3 1 ...
##  $ Temp.Out   : num  1 2 2 1 2 2 2 2 3 1 ...
##  $ TORQUE     : int  3 4 4 3 3 4 4 4 3 2 ...
##  $ Torque.max.: num  412 413 414 414 414 ...
##  $ Torque_2   : num  0.1 0.9 1.2 1 0.8 1.6 0.7 1.1 1.3 0.1 ...
##  $ TVD        : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ WOB        : num  NA NA NA NA NA NA NA NA NA NA ...

Table data

rmarkdown::paged_table(drilling)

The data has 54,196 rows and 27 columns. wellName is a unique identifier for each well. Our target variable is ROP, which signifies the Rate operation Penetration, means how fast the drill well could penetrate to the TD Target. We will use other variable except the wellName.

Before we go further, do some EDA: 1. first we need to make sure that our data is clean and will be useful. If you watch closely, there are some problems with the categorical variable. We need to delete some columns are: “Torque.max.” , “Torque_2” , “MW Out” and “datasetName”.

drilling_clean <- drilling %>%
                select(-c("Torque.max." , "Torque_2" , MW.Out , datasetName))
glimpse(drilling_clean)
## Rows: 54,196
## Columns: 23
## $ wellName    <fct> Loris 12H, Loris 12H, Loris 12H, Loris 12H, Loris 12H, Lor…
## $ Depth       <dbl> 401.0, 401.5, 402.0, 402.5, 403.0, 403.5, 404.0, 404.5, 40…
## $ Bitm        <dbl> 0.64, 1.16, 1.63, 2.11, 2.67, 3.16, 3.60, 4.11, 4.63, 5.07…
## $ BitRevs     <int> 280, 305, 326, 347, 371, 392, 411, 432, 453, 471, 483, 500…
## $ DrillHrs    <dbl> 0.01, 0.02, 0.02, 0.03, 0.03, 0.04, 0.04, 0.05, 0.06, 0.06…
## $ Dxc         <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ FlowRate    <dbl> 413.94, 414.54, 413.74, 414.94, 415.34, 414.94, 415.24, 41…
## $ MSE         <dbl> 438.70, 418.52, 407.39, 437.37, 471.06, 432.84, 432.84, 37…
## $ MSE_DAVG    <dbl> 1.033, 1.025, 1.021, 1.017, 1.022, 1.023, 1.023, 1.022, 1.…
## $ MW.In       <dbl> 0.216, 0.216, 0.216, 0.216, 0.216, 0.216, 0.216, 0.216, 0.…
## $ POWER       <dbl> 2099.75, 2094.23, 2089.72, 2098.98, 2082.44, 2083.50, 2072…
## $ Pump.time   <dbl> 75.89, 78.36, 80.46, 81.38, 84.61, 85.54, 86.14, 87.28, 88…
## $ PVT         <dbl> 1.99, 2.00, 2.00, 2.01, 2.02, 2.02, 2.03, 2.03, 2.04, 2.04…
## $ ROP         <dbl> 61, 60, 60, 61, 61, 60, 60, 60, 61, 61, 61, 61, 61, 61, 61…
## $ RotatingHrs <dbl> 61, 60, 60, 61, 61, 60, 60, 60, 61, 61, 61, 61, 61, 61, 60…
## $ RPM         <dbl> 226.07, 226.73, 226.44, 226.83, 227.13, 227.35, 227.18, 22…
## $ RPM.SURF    <dbl> 32.9, 32.9, 33.0, 33.0, 32.9, 32.9, 33.0, 33.0, 33.0, 32.9…
## $ SPP         <dbl> 30.8, 30.8, 30.8, 30.8, 30.7, 30.8, 30.8, 30.8, 30.8, 30.8…
## $ Temp.In     <dbl> 2, 3, 3, 2, 2, 3, 3, 3, 3, 1, 2, 2, 1, 2, 3, 2, 2, 2, 2, 2…
## $ Temp.Out    <dbl> 1, 2, 2, 1, 2, 2, 2, 2, 3, 1, 1, 0, 1, 1, 1, 1, 2, 1, 1, 1…
## $ TORQUE      <int> 3, 4, 4, 3, 3, 4, 4, 4, 3, 2, 2, 4, 2, 3, 6, 3, 3, 2, 2, 4…
## $ TVD         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ WOB         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
  1. Remove rows where the value of ROP does not make sense
boxplot(drilling_clean$ROP)

drilling_ROP <- subset(drilling_clean, drilling_clean$ROP > (quantile(drilling_clean$ROP, .25) - 1.5*IQR(drilling_clean$ROP)) & drilling_clean$ROP < (quantile(drilling_clean$ROP, .75) + 1.5*IQR(drilling_clean$ROP)))

boxplot(drilling_ROP$ROP)

# Check missing Value

colSums(is.na(drilling_ROP))
##    wellName       Depth        Bitm     BitRevs    DrillHrs         Dxc 
##           0           0           0           0           0           0 
##    FlowRate         MSE    MSE_DAVG       MW.In       POWER   Pump.time 
##           0           0           0           0           0           0 
##         PVT         ROP RotatingHrs         RPM    RPM.SURF         SPP 
##           0           0           0           0           0           0 
##     Temp.In    Temp.Out      TORQUE         TVD         WOB 
##           0           0           0        8313        8313

Delete missing value “NA”

drilling_ROP_na <- na.omit(drilling_ROP)

looking Correlation each variable

library(GGally)
ggcorr(drilling_ROP_na, label = TRUE, label_size = 2.9, hjust = 1, layout.exp = 2)
## Warning in ggcorr(drilling_ROP_na, label = TRUE, label_size = 2.9, hjust = 1, :
## data in column(s) 'wellName' are not numeric and were ignored

By its run data correlation for realted to ROP each variable are close to low number ~ 0 (weak correlation)

boxplot(drilling_ROP_na)

### Regression Linear Model

Create Model Regression Linear with target of ROP for all variable (Predictors).

Before we make the model, we need to split the data into train dataset and test dataset. We will use the train dataset to train the linear regression model. The test dataset will be used as a comparasion and see if the model get overfit and can not predict new data that hasn’t been seen during training phase. We will 80% of the data as the training data and the rest of it as the testing data.

set.seed(123)
samplesize <- round(0.8 * nrow(drilling_ROP_na), 0)
index <- sample(seq_len(nrow(drilling_ROP_na)), size = samplesize)

data_train <- drilling_ROP_na[index, ]
data_test <- drilling_ROP_na[-index, ]

#Linear Regresssion

set.seed(123)
m_drill <- lm(ROP ~ ., data = drilling_ROP_na)

summary(m_drill)
## 
## Call:
## lm(formula = ROP ~ ., data = drilling_ROP_na)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -108.471   -9.735   -0.647    8.341  269.206 
## 
## Coefficients:
##                            Estimate     Std. Error t value             Pr(>|t|)
## (Intercept)           92.8183826729   2.3229110121  39.958 < 0.0000000000000002
## wellNameLoris 11H    -22.3088771748   1.4906828793 -14.966 < 0.0000000000000002
## wellNameLoris 12H L1   2.9675075569   0.5448414385   5.447         0.0000000516
## wellNameLoris 13H L2   9.6736126680   0.5207119245  18.578 < 0.0000000000000002
## wellNameLoris 14H L1  10.1909007594   0.5148013491  19.796 < 0.0000000000000002
## wellNameLoris 14H L2   9.5383787345   0.5195811442  18.358 < 0.0000000000000002
## Depth                  0.0099396564   0.0004559951  21.798 < 0.0000000000000002
## Bitm                   0.0278076886   0.0008490603  32.751 < 0.0000000000000002
## BitRevs                0.0001649446   0.0000068915  23.935 < 0.0000000000000002
## DrillHrs              -1.1369437722   0.0643743672 -17.661 < 0.0000000000000002
## Dxc                   -0.0002174694   0.0007743364  -0.281              0.77883
## FlowRate               0.0229900208   0.0016042725  14.330 < 0.0000000000000002
## MSE                   -0.0000226918   0.0000003127 -72.565 < 0.0000000000000002
## MSE_DAVG              -0.0001409718   0.0000027446 -51.364 < 0.0000000000000002
## MW.In                -49.0883438775   2.0203951517 -24.296 < 0.0000000000000002
## POWER                  0.0101436874   0.0018232700   5.563         0.0000000266
## Pump.time             -0.0087931459   0.0006879325 -12.782 < 0.0000000000000002
## PVT                    0.0025260587   0.0002602550   9.706 < 0.0000000000000002
## RotatingHrs           -0.3135635066   0.0057540301 -54.495 < 0.0000000000000002
## RPM                    0.3995917005   0.0052969169  75.439 < 0.0000000000000002
## RPM.SURF               0.0020918246   0.0007051931   2.966              0.00302
## SPP                   -0.0238429965   0.0004090453 -58.289 < 0.0000000000000002
## Temp.In                1.3378162816   0.0310538202  43.081 < 0.0000000000000002
## Temp.Out              -0.0023783297   0.0017229608  -1.380              0.16748
## TORQUE                -0.0134718503   0.0358476624  -0.376              0.70706
## TVD                   -0.0539420723   0.0011693056 -46.132 < 0.0000000000000002
## WOB                    0.3236265559   0.0107636470  30.067 < 0.0000000000000002
##                         
## (Intercept)          ***
## wellNameLoris 11H    ***
## wellNameLoris 12H L1 ***
## wellNameLoris 13H L2 ***
## wellNameLoris 14H L1 ***
## wellNameLoris 14H L2 ***
## Depth                ***
## Bitm                 ***
## BitRevs              ***
## DrillHrs             ***
## Dxc                     
## FlowRate             ***
## MSE                  ***
## MSE_DAVG             ***
## MW.In                ***
## POWER                ***
## Pump.time            ***
## PVT                  ***
## RotatingHrs          ***
## RPM                  ***
## RPM.SURF             ** 
## SPP                  ***
## Temp.In              ***
## Temp.Out                
## TORQUE                  
## TVD                  ***
## WOB                  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.28 on 41683 degrees of freedom
## Multiple R-squared:  0.6481, Adjusted R-squared:  0.6479 
## F-statistic:  2953 on 26 and 41683 DF,  p-value: < 0.00000000000000022
plot(drilling_ROP_na$ROP, drilling_ROP_na$"RotatingHrs")
abline(m_drill$coefficients[1],m_drill$coefficients[19], col = "red")

by its summary data our Adjusted R-squared = 0.6479

Model Performance

The performance of our model (how well our model predict the target variable) can be calculated using root mean squared error:

RMSE is better than MAE or mean absolute error, because RMSE squared the difference between the actual values and the predicted values, meaning that prediction with higher error will be penalized greatly. This metric is often used to compare two or more alternative models, even though it is harder to interpret than MAE. We can use the RMSE () functions from caret package. Below is the first model (with complete variables) performance.

lm_pred_ROP <- predict(m_drill, newdata = data_test %>% select(-ROP))

# RMSE of train dataset
RMSE(pred = m_drill$fitted.values, obs = data_train$ROP)
## Warning in pred - obs: longer object length is not a multiple of shorter object
## length
## [1] 37.28579
RMSE(pred = lm_pred_ROP, obs = data_test$ROP)
## [1] 17.38003

Model Regresi Backward

m1 <- lm(ROP ~ . -wellName, drilling_ROP_na)

summary(m1)
## 
## Call:
## lm(formula = ROP ~ . - wellName, data = drilling_ROP_na)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -109.721   -9.623   -0.553    8.124  261.210 
## 
## Coefficients:
##                   Estimate     Std. Error t value             Pr(>|t|)    
## (Intercept) 110.9354679547   2.1083006175  52.618 < 0.0000000000000002 ***
## Depth         0.0118113833   0.0004115474  28.700 < 0.0000000000000002 ***
## Bitm          0.0268187568   0.0008400912  31.924 < 0.0000000000000002 ***
## BitRevs       0.0001765254   0.0000063635  27.740 < 0.0000000000000002 ***
## DrillHrs     -1.2070700074   0.0616265295 -19.587 < 0.0000000000000002 ***
## Dxc           0.0002271037   0.0007834526   0.290              0.77191    
## FlowRate      0.0259132584   0.0015903521  16.294 < 0.0000000000000002 ***
## MSE          -0.0000219845   0.0000003085 -71.272 < 0.0000000000000002 ***
## MSE_DAVG     -0.0001462033   0.0000026865 -54.422 < 0.0000000000000002 ***
## MW.In       -58.4557469472   2.0098292822 -29.085 < 0.0000000000000002 ***
## POWER         0.0355359583   0.0012062442  29.460 < 0.0000000000000002 ***
## Pump.time    -0.0173044798   0.0005324309 -32.501 < 0.0000000000000002 ***
## PVT           0.0044569285   0.0002311082  19.285 < 0.0000000000000002 ***
## RotatingHrs  -0.3449530924   0.0056420640 -61.140 < 0.0000000000000002 ***
## RPM           0.3695260186   0.0049447627  74.731 < 0.0000000000000002 ***
## RPM.SURF      0.0019144472   0.0006783154   2.822              0.00477 ** 
## SPP          -0.0225233203   0.0004068181 -55.365 < 0.0000000000000002 ***
## Temp.In       1.1745853378   0.0289799949  40.531 < 0.0000000000000002 ***
## Temp.Out     -0.0017856727   0.0017417975  -1.025              0.30528    
## TORQUE       -0.4287778060   0.0305572951 -14.032 < 0.0000000000000002 ***
## TVD          -0.0560881776   0.0010715917 -52.341 < 0.0000000000000002 ***
## WOB           0.2942346569   0.0104668878  28.111 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.49 on 41688 degrees of freedom
## Multiple R-squared:  0.6394, Adjusted R-squared:  0.6392 
## F-statistic:  3519 on 21 and 41688 DF,  p-value: < 0.00000000000000022

Evaluasi Model

Normalitas

hist(m1$residuals, breaks = 20)

# Heteroscedasticity

library(lmtest)
plot(drilling_ROP_na$ROP, m1$residuals)
abline(h = 0, col = "red")

bptest(m1)
## 
##  studentized Breusch-Pagan test
## 
## data:  m1
## BP = 11163, df = 21, p-value < 0.00000000000000022

Variance Inflation Factor (Multicollinearity)

library(car)

vif(m1)
##       Depth        Bitm     BitRevs    DrillHrs         Dxc    FlowRate 
##   25.992583   20.521903   65.110104  119.147827    1.002782   25.661950 
##         MSE    MSE_DAVG       MW.In       POWER   Pump.time         PVT 
##    1.692650    1.815456    2.906143    3.680511   20.201856    5.824745 
## RotatingHrs         RPM    RPM.SURF         SPP     Temp.In    Temp.Out 
##    9.155503    3.043153   41.711447   25.911137    3.013968    1.040778 
##      TORQUE         TVD         WOB 
##    8.804101    5.998570    2.025363

There are Multicollinearity between variable (each predictor were dependen).

some variables has value >10 are: - Depth - Bitm - BitRevs - DrillHrs - FlowRate - Pump.time - RPM.SURF - SPP

by its case, the model could be better using classification model (C2).

Conclusion & Suggestion

Model m_drill has value of Adjusted R-squared = 0.6479. the models could not be modeled using Linear Regression. its must be using Classification data model.