Multiple Linear Regression of Core Permeability by using BOX-COX transformation, and Log transformation
Application of Karpur data set, West Africa
Hussein Malik Chasib, 4th year student, BUOG, Basrah, Iraq
1/8/2022
…………………………………………………………….
There are several ways to provide the data used to characterize the reservoirs, which are different from each other in the scale. The different scales of each method provide varying estimates of permeability (or other parameters), so this will affect the data distribution, and which affects the accuracy of construction geostatistical modeling. The most conventional approach for combining the core measurements, well log data into permeability modeling is the Multiple Linear Regression, so we will use it to transform data from non-normal distribution into normal distribution to obtain accurate models.
…………………………………………………………….
Loading some data:
setwd("C:/Users/Hussein Malik/Desktop/R data")
data <- read.csv("karpur.csv", header=TRUE)
data <-read.csv("C:/Users/Hussein Malik/Desktop/R data/karpur.csv",header = TRUE)Show the Loaded data as ” head, tail, range, summary, and dimension” :
head(data)## depth caliper ind.deep ind.med gamma phi.N R.deep R.med SP
## 1 5667.0 8.685 618.005 569.781 98.823 0.410 1.618 1.755 -56.587
## 2 5667.5 8.686 497.547 419.494 90.640 0.307 2.010 2.384 -61.916
## 3 5668.0 8.686 384.935 300.155 78.087 0.203 2.598 3.332 -55.861
## 4 5668.5 8.686 278.324 205.224 66.232 0.119 3.593 4.873 -41.860
## 5 5669.0 8.686 183.743 131.155 59.807 0.069 5.442 7.625 -34.934
## 6 5669.5 8.686 109.512 75.633 57.109 0.048 9.131 13.222 -39.769
## density.corr density phi.core k.core Facies phi.core.fr.
## 1 -0.033 2.205 33.9000 2442.590 F1 0.339000
## 2 -0.067 2.040 33.4131 3006.989 F1 0.334131
## 3 -0.064 1.888 33.1000 3370.000 F1 0.331000
## 4 -0.053 1.794 34.9000 2270.000 F1 0.349000
## 5 -0.054 1.758 35.0644 2530.758 F1 0.350644
## 6 -0.058 1.759 35.3152 2928.314 F1 0.353152tail(data)## depth caliper ind.deep ind.med gamma phi.N R.deep R.med SP
## 814 6080.5 8.578 683.847 672.125 24.003 0.208 1.462 1.488 -33.714
## 815 6081.0 8.590 678.002 669.171 27.855 0.214 1.475 1.494 -4.808
## 816 6081.5 8.588 668.633 661.949 32.591 0.228 1.496 1.511 -4.727
## 817 6082.0 8.588 656.328 647.918 38.547 0.243 1.524 1.543 -21.390
## 818 6082.5 8.588 643.216 628.536 45.555 0.256 1.555 1.591 -31.597
## 819 6083.0 8.588 631.098 608.163 52.244 0.265 1.584 1.644 -36.109
## density.corr density phi.core k.core Facies phi.core.fr.
## 814 0.003 2.147 26.9000 1300.790 F5 0.269000
## 815 -0.002 2.162 25.7547 1249.706 F5 0.257547
## 816 -0.008 2.158 24.5569 1196.282 F5 0.245569
## 817 -0.006 2.136 23.3592 1142.859 F5 0.233592
## 818 -0.002 2.115 22.1614 1089.435 F5 0.221614
## 819 -0.003 2.109 20.9637 1036.012 F5 0.209637range(data$depth)## [1] 5667 6083summary(data)## depth caliper ind.deep ind.med
## Min. :5667 Min. :8.487 Min. : 6.532 Min. : 9.386
## 1st Qu.:5769 1st Qu.:8.556 1st Qu.: 28.799 1st Qu.: 27.892
## Median :5872 Median :8.588 Median :217.849 Median :254.383
## Mean :5873 Mean :8.622 Mean :275.357 Mean :273.357
## 3rd Qu.:5977 3rd Qu.:8.686 3rd Qu.:566.793 3rd Qu.:544.232
## Max. :6083 Max. :8.886 Max. :769.484 Max. :746.028
## gamma phi.N R.deep R.med
## Min. : 16.74 Min. :0.0150 Min. : 1.300 Min. : 1.340
## 1st Qu.: 40.89 1st Qu.:0.2030 1st Qu.: 1.764 1st Qu.: 1.837
## Median : 51.37 Median :0.2450 Median : 4.590 Median : 3.931
## Mean : 53.42 Mean :0.2213 Mean : 24.501 Mean : 21.196
## 3rd Qu.: 62.37 3rd Qu.:0.2640 3rd Qu.: 34.724 3rd Qu.: 35.853
## Max. :112.40 Max. :0.4100 Max. :153.085 Max. :106.542
## SP density.corr density phi.core
## Min. :-73.95 Min. :-0.067000 Min. :1.758 Min. :15.70
## 1st Qu.:-42.01 1st Qu.:-0.016000 1st Qu.:2.023 1st Qu.:23.90
## Median :-32.25 Median :-0.007000 Median :2.099 Median :27.60
## Mean :-30.98 Mean :-0.008883 Mean :2.102 Mean :26.93
## 3rd Qu.:-19.48 3rd Qu.: 0.002000 3rd Qu.:2.181 3rd Qu.:30.70
## Max. : 25.13 Max. : 0.089000 Max. :2.387 Max. :36.30
## k.core Facies phi.core.fr.
## Min. : 0.42 Length:819 Min. :0.1570
## 1st Qu.: 657.33 Class :character 1st Qu.:0.2390
## Median : 1591.22 Mode :character Median :0.2760
## Mean : 2251.91 Mean :0.2693
## 3rd Qu.: 3046.82 3rd Qu.:0.3070
## Max. :15600.00 Max. :0.3630dim(data)## [1] 819 15Plot a histogram to the measured core permeabilit :
par(mfrow=c(1,1))
hist(data$k.core, main='Histogram of Permeability', xlab='Permeability, md', col='darkred')
From the above histograms, we see it has a non-normal distribution.
——————————————————————————–
Multiple Linear Regression of Core Permeability as a function of other data :
m1 <- lm(k.core ~ .-1, data = data)summary(m1)##
## Call:
## lm(formula = k.core ~ . - 1, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5585.6 -568.9 49.2 476.5 8928.4
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## depth 8.544e+00 1.785e+00 4.786 2.02e-06 ***
## caliper 1.413e+03 1.019e+03 1.387 0.165789
## ind.deep -2.418e-01 2.354e+00 -0.103 0.918220
## ind.med 1.224e+00 2.585e+00 0.473 0.636062
## gamma -4.583e+01 6.010e+00 -7.626 6.88e-14 ***
## phi.N -2.010e+03 1.476e+03 -1.362 0.173540
## R.deep -2.344e+01 6.288e+00 -3.727 0.000207 ***
## R.med 5.643e+01 9.065e+00 6.225 7.76e-10 ***
## SP -7.125e+00 3.145e+00 -2.266 0.023736 *
## density.corr -2.567e+03 4.809e+03 -0.534 0.593602
## density 2.319e+03 1.173e+03 1.976 0.048458 *
## phi.core 1.921e+02 2.282e+01 8.418 < 2e-16 ***
## FaciesF1 -6.783e+04 1.760e+04 -3.853 0.000126 ***
## FaciesF10 -6.694e+04 1.765e+04 -3.793 0.000160 ***
## FaciesF2 -6.691e+04 1.761e+04 -3.800 0.000156 ***
## FaciesF3 -6.740e+04 1.767e+04 -3.815 0.000147 ***
## FaciesF5 -6.709e+04 1.767e+04 -3.798 0.000157 ***
## FaciesF7 -6.788e+04 1.764e+04 -3.848 0.000129 ***
## FaciesF8 -6.901e+04 1.770e+04 -3.900 0.000104 ***
## FaciesF9 -7.080e+04 1.779e+04 -3.980 7.52e-05 ***
## phi.core.fr. NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1262 on 799 degrees of freedom
## Multiple R-squared: 0.8457, Adjusted R-squared: 0.8418
## F-statistic: 218.9 on 20 and 799 DF, p-value: < 2.2e-16# Plot the model (m1)
par(mfrow=c(2,2))
plot(m1)
We note from the above plot that there are some inaccurate (anomalous) or ineffective points that may cause errors in the model, so they should be reduced and then reconstruct the relationship.
m1.red <- lm(k.core ~ depth+gamma+R.deep+R.med+SP+density+phi.core+Facies-1, data = data)summary(m1.red)##
## Call:
## lm(formula = k.core ~ depth + gamma + R.deep + R.med + SP + density +
## phi.core + Facies - 1, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5461.7 -545.5 37.0 505.0 9072.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## depth 8.088 1.059 7.639 6.23e-14 ***
## gamma -51.915 4.526 -11.471 < 2e-16 ***
## R.deep -22.467 6.261 -3.588 0.000353 ***
## R.med 48.859 8.730 5.597 2.99e-08 ***
## SP -6.490 3.118 -2.082 0.037671 *
## density 2013.607 1047.054 1.923 0.054818 .
## phi.core 188.002 21.791 8.628 < 2e-16 ***
## FaciesF1 -51703.918 5802.150 -8.911 < 2e-16 ***
## FaciesF10 -50893.923 5917.391 -8.601 < 2e-16 ***
## FaciesF2 -51009.565 5854.625 -8.713 < 2e-16 ***
## FaciesF3 -51322.213 5871.731 -8.741 < 2e-16 ***
## FaciesF5 -51174.044 6008.663 -8.517 < 2e-16 ***
## FaciesF7 -52302.690 5971.976 -8.758 < 2e-16 ***
## FaciesF8 -53362.629 6033.789 -8.844 < 2e-16 ***
## FaciesF9 -54796.231 6100.716 -8.982 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1266 on 804 degrees of freedom
## Multiple R-squared: 0.8436, Adjusted R-squared: 0.8407
## F-statistic: 289.1 on 15 and 804 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(m1.red)
now, we will predict by using the (m1) model to obtain predicted k.core
k.core.pred1 <- predict(m1.red,data)by using “cbind” function we can see the difference between the “measured kcore” and “predicted k.core”
#cbind(data$k.core ,k.core.pred1 )also we can see the difference between the “measured kcore” and “predicted k.core” from their plots
par(mfrow=c(1,3))
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(data$k.core),type="l", col="red", lwd = 5, pch=17, xlab='Measured',
ylab='Depth,m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred1),type="l", col="blue", lwd = 5, pch=17, xlab='Predicted',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2,)
grid()
# Matching the two curves for ease of comparison
par(mfrow=c(1,1))
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(data$k.core),type="l", col="red", lwd = 5, xlab='Permeability',
ylab='Depth,m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
par(new = TRUE)
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred1),type="l", col="blue", lwd = 5, xlab='',
ylab='Depth,m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2, main='R-sq=8436')
grid()
legend('topright', legend=c("Observed", "Predicted"), lty=c(1,1), col=c("red", "blue"))
addition adjusted R square & root mine square error for m1 to our matching plot after we calculat them.
AdjR.sq1 <- 1-sum((k.core.pred1 - data$k.core)^2)/sum((data$k.core-mean(data$k.core))^2)
AdjR.sq1## [1] 0.6847518mspe.model1 <- sqrt(sum((k.core.pred1 - data$k.core)^2)/nrow(data))
mspe.model1## [1] 1254.46par(mfrow=c(1,1))
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(data$k.core),type="p", col="red", lwd = 5, pch=16, xlab='Permeability',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
par(new = TRUE)
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred1),type="p", col="blue", lwd = 5, pch=15, xlab='',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2, main='AdjR.sq1=0.6847 & RMSE=1254')
grid()
legend('topright', legend=c("Observed", "Predicted"), pch=c(16,15), col=c("red", "blue"))
——————————————————————————–
Construct core permeability modeling as log transformation :
hist(log10(data$k.core) ,main='Histogram of log Permeability', xlab='log Permeability, md', col='dark gray')
m2<-lm(log10(k.core) ~ .-1,data=data)summary(m2)##
## Call:
## lm(formula = log10(k.core) ~ . - 1, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5804 -0.1138 0.0322 0.1529 0.7384
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## depth 0.0007425 0.0004718 1.574 0.1160
## caliper -0.4605945 0.2693103 -1.710 0.0876 .
## ind.deep -0.0007951 0.0006222 -1.278 0.2017
## ind.med 0.0007137 0.0006833 1.044 0.2966
## gamma -0.0091269 0.0015885 -5.746 1.30e-08 ***
## phi.N -1.7628155 0.3901024 -4.519 7.16e-06 ***
## R.deep -0.0025878 0.0016620 -1.557 0.1199
## R.med 0.0044073 0.0023960 1.839 0.0662 .
## SP -0.0016935 0.0008312 -2.037 0.0419 *
## density.corr 1.4462633 1.2712045 1.138 0.2556
## density 1.6148374 0.3100921 5.208 2.44e-07 ***
## phi.core 0.0948634 0.0060329 15.724 < 2e-16 ***
## FaciesF1 -2.3461877 4.6532000 -0.504 0.6143
## FaciesF10 -2.2675417 4.6652182 -0.486 0.6271
## FaciesF2 -2.3646211 4.6544047 -0.508 0.6116
## FaciesF3 -2.3769425 4.6697172 -0.509 0.6109
## FaciesF5 -2.2367684 4.6697432 -0.479 0.6321
## FaciesF7 -2.0650257 4.6623059 -0.443 0.6579
## FaciesF8 -2.4438111 4.6776859 -0.522 0.6015
## FaciesF9 -2.7023993 4.7023453 -0.575 0.5657
## phi.core.fr. NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3335 on 799 degrees of freedom
## Multiple R-squared: 0.9891, Adjusted R-squared: 0.9888
## F-statistic: 3613 on 20 and 799 DF, p-value: < 2.2e-16# Plot the model (m2)
par(mfrow=c(2,2))
plot(m2)
We note from the above plot that there are some inaccurate (anomalous) or ineffective points that may cause errors in the model, so they should be reduced and then reconstruct the relationship.But this model”m2” is more accurate than model “m1”.
m2.red <- lm(log10(k.core) ~ caliper+gamma+phi.N+SP+density+phi.core+Facies-1,data=data)summary(m2.red)##
## Call:
## lm(formula = log10(k.core) ~ caliper + gamma + phi.N + SP + density +
## phi.core + Facies - 1, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.60108 -0.11891 0.03253 0.15792 0.70925
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## caliper -0.840136 0.182354 -4.607 4.74e-06 ***
## gamma -0.010369 0.001337 -7.754 2.68e-14 ***
## phi.N -1.596764 0.359392 -4.443 1.01e-05 ***
## SP -0.001585 0.000793 -1.998 0.04600 *
## density 1.762145 0.287056 6.139 1.31e-09 ***
## phi.core 0.093920 0.005931 15.835 < 2e-16 ***
## FaciesF1 4.981486 1.783281 2.793 0.00534 **
## FaciesF10 5.101106 1.793029 2.845 0.00455 **
## FaciesF2 4.995149 1.796855 2.780 0.00556 **
## FaciesF3 4.967152 1.799799 2.760 0.00591 **
## FaciesF5 5.158638 1.760353 2.930 0.00348 **
## FaciesF7 5.272425 1.759300 2.997 0.00281 **
## FaciesF8 4.952389 1.753512 2.824 0.00486 **
## FaciesF9 4.687264 1.762588 2.659 0.00799 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3346 on 805 degrees of freedom
## Multiple R-squared: 0.9889, Adjusted R-squared: 0.9887
## F-statistic: 5128 on 14 and 805 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(m2.red)
now, we will predict by using the (m2) model to obtain predicted k.core as log transformation.
k.core.pred2 <- predict(m2.red,data)back transformation to the predicted core permeability
k.core.pred2 <- 10^(k.core.pred2) # this back transformation to log By using “cbind” function we can see the difference between the “measured kcore” and “predicted k.core as log transformation”
# cbind(data$k.core ,k.core.pred2 )Also we can see the difference between the “measured kcore” and “predicted k.core as log transformation” from their plots
par(mfrow=c(1,3))
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(data$k.core),type="l", col="red", lwd = 5, xlab='Measured',
ylab='Depth,m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred2),type="l", col="green", lwd = 5, xlab='Predicted',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2,)
grid()
# Matching the two curves for ease of comparison
par(mfrow=c(1,1))
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(data$k.core),type="l", col="red", lwd = 5, xlab='Permeability',
ylab='Depth,m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
par(new = TRUE)
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred2),type="l", col="green", lwd = 5, xlab='',
ylab='Depth,m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2, main='R-sq=8436')
grid()
legend('topright', legend=c("Observed", "Predicted"), lty=c(1,1), col=c("red", "green"))
addition adjusted R square & root mine square error for m1 to our matching plot after we calculat them.
AdjR.sq2 <- 1-sum((k.core.pred2 - data$k.core)^2)/sum((data$k.core-mean(data$k.core))^2)
AdjR.sq2## [1] 0.6366977mspe.model2 <- sqrt(sum((k.core.pred2 - data$k.core)^2)/nrow(data))
mspe.model2## [1] 1346.681par(mfrow=c(1,1))
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(data$k.core),type="p", col="red", lwd = 5, pch=16, xlab='Permeability',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
par(new = TRUE)
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred2),type="p", col="green", lwd = 5, pch=17, xlab='',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2, main='AdjR.sq1=0.6367 & RMSE=1346')
grid()
legend('topright', legend=c("Observed", "Predicted"), pch=c(16,17), col=c("red", "green"))
——————————————————————————–
Construct core permeability modeling as BOX-COX transformation by using boxcox function :
#Loading required package: MASS
require(MASS)## Loading required package: MASSlibrary(MASS)y <- data$k.core
x <- data$phi.core
bc <- boxcox(y ~ x, data=data)
lambda <- bc$x[which.max(bc$y)]
lambda## [1] 0.2222222# calculating k.core depending on box-cox transformation "k.corebc"
k.corebc <- ((y^lambda-1)/lambda)
hist(k.corebc)
# Construct core permeability modeling as BOX-COX transformation by using box-cox function
m3 <- lm(k.corebc ~ .-1,data=data)summary(m3)##
## Call:
## lm(formula = k.corebc ~ . - 1, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6183 -0.7390 0.1217 1.0871 4.9292
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## depth -2.837e-03 2.923e-03 -0.971 0.331955
## caliper -4.670e+00 1.647e+00 -2.836 0.004688 **
## ind.deep -4.252e-03 3.800e-03 -1.119 0.263538
## ind.med 2.480e-03 4.174e-03 0.594 0.552529
## gamma -3.470e-02 1.005e-02 -3.453 0.000584 ***
## phi.N -1.029e+01 2.385e+00 -4.315 1.79e-05 ***
## R.deep 3.507e-03 1.024e-02 0.343 0.732063
## R.med -2.005e-02 1.498e-02 -1.338 0.181359
## SP -4.665e-03 5.093e-03 -0.916 0.359998
## density.corr 1.113e+01 7.766e+00 1.433 0.152239
## density 8.304e+00 1.899e+00 4.374 1.38e-05 ***
## phi.core 5.314e-01 3.845e-02 13.821 < 2e-16 ***
## k.core 1.650e-03 5.711e-05 28.886 < 2e-16 ***
## FaciesF1 4.457e+01 2.868e+01 1.554 0.120574
## FaciesF10 4.473e+01 2.875e+01 1.556 0.120154
## FaciesF2 4.413e+01 2.868e+01 1.538 0.124356
## FaciesF3 4.395e+01 2.878e+01 1.527 0.127158
## FaciesF5 4.466e+01 2.878e+01 1.552 0.121080
## FaciesF7 4.666e+01 2.874e+01 1.624 0.104864
## FaciesF8 4.510e+01 2.884e+01 1.564 0.118309
## FaciesF9 4.441e+01 2.900e+01 1.531 0.126140
## phi.core.fr. NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.037 on 798 degrees of freedom
## Multiple R-squared: 0.9891, Adjusted R-squared: 0.9888
## F-statistic: 3432 on 21 and 798 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(m3)
We note from the above plot that there are some inaccurate (anomalous) or ineffective points that may cause errors in the model, so they should be reduced and then reconstruct the relationship.
m3.red <- lm(k.corebc ~ caliper+gamma+phi.N+density+phi.core+Facies-1,data=data)summary(m3.red)##
## Call:
## lm(formula = k.corebc ~ caliper + gamma + phi.N + density + phi.core +
## Facies - 1, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.9794 -1.4494 0.2455 1.7020 8.1339
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## caliper -7.93274 1.61041 -4.926 1.02e-06 ***
## gamma -0.13535 0.01179 -11.476 < 2e-16 ***
## phi.N -10.34785 3.17510 -3.259 0.00116 **
## density 12.82088 2.53438 5.059 5.23e-07 ***
## phi.core 0.85651 0.05240 16.344 < 2e-16 ***
## FaciesF1 46.56838 15.72569 2.961 0.00315 **
## FaciesF10 48.44570 15.81731 3.063 0.00227 **
## FaciesF2 47.61675 15.84673 3.005 0.00274 **
## FaciesF3 46.59335 15.87771 2.935 0.00344 **
## FaciesF5 48.82275 15.52754 3.144 0.00173 **
## FaciesF7 47.87893 15.52847 3.083 0.00212 **
## FaciesF8 45.68235 15.47125 2.953 0.00324 **
## FaciesF9 42.39103 15.55060 2.726 0.00655 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.956 on 806 degrees of freedom
## Multiple R-squared: 0.9767, Adjusted R-squared: 0.9763
## F-statistic: 2600 on 13 and 806 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(m3.red)
now, we will predict by using the (m2) model to obtain predicted k.core as BOX-COX transformation.
k.core.pred3 <- predict(m3.red,data)
# Back transformation to the predicted k.core
k.core.pred3 <-(1+(lambda*k.corebc))^(1/lambda)by using “cbind” function we can see the difference between the “measured kcore” and “predicted k.core as BOX-COX transformation”
#cbind(data$k.core ,k.core.pred3 )also we can see the difference between the “measured kcore” and “predicted k.core as BOX-COX transformation” from their plots
par(mfrow=c(1,2))
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(data$k.core),type="l", col="red", lwd = 5, xlab='Measured',
ylab='Depth,m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred3),type="l", col="gold", lwd = 5, xlab='Predicted',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2,)
grid()
# Matching the two curves for ease of comparison
par(mfrow=c(1,1))
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(data$k.core),type="l", col="red", lwd = 5, xlab='Permeability',
ylab='Depth,m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
par(new = TRUE)
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred3),type="l", col="gold", lwd = 5, xlab='',
ylab='Depth,m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2, main='R-sq=8436')
grid()
legend('topright', legend=c("Observed", "Predicted"), lty=c(1,1), col=c("red", "gold"))
addition adjusted R square & root mine square error for m1 to our matching plot after we calculat them.
AdjR.sq3 <- 1-sum((k.core.pred3 - data$k.core)^2)/sum((data$k.core-mean(data$k.core))^2)
AdjR.sq3## [1] 1mspe.m3 <- sqrt(sum((k.core.pred3 - data$k.core)^2)/nrow(data))
mspe.m3## [1] 1.560247e-12par(mfrow=c(1,1))
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(data$k.core),type="p", col="red", lwd = 5, pch=16, xlab='Permeability',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
par(new = TRUE)
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred3),type="p", col="gold", lwd = 5, pch=11, xlab='',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2, main='AdjR.sq1=0.9999 & RMSE=0')
grid()
legend('topright', legend=c("Observed", "Predicted"), pch=c(16,15), col=c("red", "gold"))
——————————————————————————–
Plot all 3 models results together
par(mfrow=c(1,1))
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(data$k.core),type="p", col="red", lwd = 5, pch=16, xlab='Permeability',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
par(new = TRUE)
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred1),type="p", col="blue", lwd = 5, pch=15, xlab='',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
par(new = TRUE)
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred2),type="p", col="green", lwd = 5, pch=17, xlab='',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
par(new = TRUE)
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(k.core.pred3),type="p", col="gold", lwd = 5, pch=11, xlab='',
ylab='Depth, m', xlim=c(0,16000), cex=1.5, cex.lab=1.5, cex.axis=1.2)
grid()
legend('topright', legend=c("Observed", "Predicted1", "Predicted2", "Predicted3"), pch=c(16,15,17,11), col=c("red", "blue", "green", "gold"))
——————————————————————————–
Add X table
# Loading required package: x table
require(xtable)## Loading required package: xtable## Warning: package 'xtable' was built under R version 4.1.2library(xtable)