Load required package:
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 ...
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…
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
drilling_ROP_na <- na.omit(drilling_ROP)
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
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
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
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
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).
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.