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.

Upload Data.

data<-read.csv("D:/Almurtaza Alaa Askar/karpur2.csv",header=T)
summary(data)
##        X             depth         caliper         ind.deep      
##  Min.   :  1.0   Min.   :5667   Min.   :8.487   Min.   :  6.532  
##  1st Qu.:205.5   1st Qu.:5769   1st Qu.:8.556   1st Qu.: 28.799  
##  Median :410.0   Median :5872   Median :8.588   Median :217.849  
##  Mean   :410.0   Mean   :5873   Mean   :8.622   Mean   :275.357  
##  3rd Qu.:614.5   3rd Qu.:5977   3rd Qu.:8.686   3rd Qu.:566.793  
##  Max.   :819.0   Max.   :6083   Max.   :8.886   Max.   :769.484  
##     ind.med            gamma            phi.N            R.deep       
##  Min.   :  9.386   Min.   : 16.74   Min.   :0.0150   Min.   :  1.300  
##  1st Qu.: 27.892   1st Qu.: 40.89   1st Qu.:0.2030   1st Qu.:  1.764  
##  Median :254.383   Median : 51.37   Median :0.2450   Median :  4.590  
##  Mean   :273.357   Mean   : 53.42   Mean   :0.2213   Mean   : 24.501  
##  3rd Qu.:544.232   3rd Qu.: 62.37   3rd Qu.:0.2640   3rd Qu.: 34.724  
##  Max.   :746.028   Max.   :112.40   Max.   :0.4100   Max.   :153.085  
##      R.med               SP          density.corr          density     
##  Min.   :  1.340   Min.   :-73.95   Min.   :-0.067000   Min.   :1.758  
##  1st Qu.:  1.837   1st Qu.:-42.01   1st Qu.:-0.016000   1st Qu.:2.023  
##  Median :  3.931   Median :-32.25   Median :-0.007000   Median :2.099  
##  Mean   : 21.196   Mean   :-30.98   Mean   :-0.008883   Mean   :2.102  
##  3rd Qu.: 35.853   3rd Qu.:-19.48   3rd Qu.: 0.002000   3rd Qu.:2.181  
##  Max.   :106.542   Max.   : 25.13   Max.   : 0.089000   Max.   :2.387  
##     phi.core         k.core            Facies             k.corel       
##  Min.   :15.70   Min.   :    0.42   Length:819         Min.   :  11.26  
##  1st Qu.:23.90   1st Qu.:  657.33   Class :character   1st Qu.: 861.89  
##  Median :27.60   Median : 1591.22   Mode  :character   Median :1833.75  
##  Mean   :26.93   Mean   : 2251.91                      Mean   :2251.91  
##  3rd Qu.:30.70   3rd Qu.: 3046.82                      3rd Qu.:3409.26  
##  Max.   :36.30   Max.   :15600.00                      Max.   :6280.63

At First , we should find the correlation between phi.core and phi.log

model1 <- lm(data$phi.core/100 ~ data$phi.N)
summary(model1)
## 
## Call:
## lm(formula = data$phi.core/100 ~ data$phi.N)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.135237 -0.030779  0.009432  0.033563  0.104025 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.30962    0.00485  63.846   <2e-16 ***
## data$phi.N  -0.18207    0.02080  -8.753   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04368 on 817 degrees of freedom
## Multiple R-squared:  0.08573,    Adjusted R-squared:  0.08462 
## F-statistic: 76.61 on 1 and 817 DF,  p-value: < 2.2e-16
plot(data$phi.N,data$phi.core,xlab="phi.log",ylab="phi.core",axes = F)
axis(2,col = "darkgreen",col.axis="black")
axis(1,col = "darkgreen",col.axis="red")
abline(model1, lwd=3, col='green')

Plot a histogram to the measured core permeabilit, From the below histograms, we see it has a non-normal distribution.

par(mfrow=c(1,1))
hist(data$k.core, main='Histogram of Permeability', xlab='Permeability, md', col='darkblue')

Multiple Linear Regression of Core Permeability as a function of other data

model1<-lm(k.core~.-1 ,data=data)
par(mfrow=c(2,2))
plot (model1)

phi.corel <- predict(model1,data)
## Warning in predict.lm(model1, data): prediction from a rank-deficient fit may be
## misleading
#cbind(karpur$phi.core/100,phi.corel)

Construction a relationship between permeability calculated from core and core porosity corrected to the log scale in order to get core premeability corrected to the log scale

model2<-lm(k.core~phi.corel+Facies-1,data=data)
summary(model2)
## 
## Call:
## lm(formula = k.core ~ phi.corel + Facies - 1, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5585.6  -574.1    45.5   483.5  8923.9 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## phi.corel  1.000e+00  4.690e-02   21.32   <2e-16 ***
## FaciesF1  -1.164e-10  1.617e+02    0.00        1    
## FaciesF10 -3.738e-11  1.026e+02    0.00        1    
## FaciesF2  -3.472e-10  4.433e+02    0.00        1    
## FaciesF3  -3.951e-11  1.766e+02    0.00        1    
## FaciesF5   3.586e-13  2.795e+02    0.00        1    
## FaciesF7  -9.498e-11  4.261e+02    0.00        1    
## FaciesF8   9.805e-11  1.851e+02    0.00        1    
## FaciesF9   4.034e-11  1.046e+02    0.00        1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1253 on 810 degrees of freedom
## Multiple R-squared:  0.8457, Adjusted R-squared:  0.844 
## F-statistic: 493.2 on 9 and 810 DF,  p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(model2)

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-16
par(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.6847518
mspe.model1 <- sqrt(sum((k.core.pred1 - data$k.core)^2)/nrow(data))
mspe.model1
## [1] 1254.46
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, 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='red')

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.58310 -0.11654  0.03049  0.14948  0.72583 
## 
## Coefficients: (1 not defined because of singularities)
##                Estimate Std. Error t value Pr(>|t|)    
## X             1.325e-02  1.336e-02   0.992   0.3213    
## depth        -2.517e-02  2.611e-02  -0.964   0.3355    
## caliper      -2.632e-01  3.348e-01  -0.786   0.4321    
## ind.deep     -8.036e-04  6.223e-04  -1.291   0.1969    
## ind.med       7.310e-04  6.835e-04   1.069   0.2852    
## gamma        -9.196e-03  1.590e-03  -5.784 1.05e-08 ***
## phi.N        -1.787e+00  3.909e-01  -4.572 5.61e-06 ***
## R.deep       -2.238e-03  1.699e-03  -1.317   0.1882    
## R.med         3.788e-03  2.476e-03   1.530   0.1264    
## SP           -1.694e-03  8.312e-04  -2.038   0.0419 *  
## density.corr  1.289e+00  1.281e+00   1.007   0.3145    
## density       1.558e+00  3.153e-01   4.942 9.41e-07 ***
## phi.core      9.453e-02  6.042e-03  15.645  < 2e-16 ***
## FaciesF1      1.429e+02  1.464e+02   0.976   0.3295    
## FaciesF10     1.429e+02  1.463e+02   0.976   0.3292    
## FaciesF2      1.428e+02  1.464e+02   0.976   0.3295    
## FaciesF3      1.428e+02  1.463e+02   0.976   0.3296    
## FaciesF5      1.429e+02  1.464e+02   0.977   0.3291    
## FaciesF7      1.431e+02  1.463e+02   0.978   0.3285    
## FaciesF8      1.427e+02  1.464e+02   0.975   0.3298    
## FaciesF9      1.424e+02  1.463e+02   0.973   0.3307    
## k.corel              NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3335 on 798 degrees of freedom
## Multiple R-squared:  0.9891, Adjusted R-squared:  0.9888 
## F-statistic:  3441 on 21 and 798 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-16
par(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.6366977
mspe.model2 <- sqrt(sum((k.core.pred2 - data$k.core)^2)/nrow(data))
mspe.model2
## [1] 1346.681
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.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: MASS
library(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.646  -0.746   0.135   1.045   4.796 
## 
## Coefficients: (1 not defined because of singularities)
##                Estimate Std. Error t value Pr(>|t|)    
## X             1.397e-01  8.148e-02   1.714 0.086896 .  
## depth        -2.758e-01  1.593e-01  -1.732 0.083724 .  
## caliper      -2.590e+00  2.044e+00  -1.267 0.205404    
## ind.deep     -4.342e-03  3.796e-03  -1.144 0.253095    
## ind.med       2.662e-03  4.170e-03   0.638 0.523372    
## gamma        -3.541e-02  1.005e-02  -3.525 0.000448 ***
## phi.N        -1.055e+01  2.387e+00  -4.418 1.13e-05 ***
## R.deep        7.209e-03  1.045e-02   0.690 0.490596    
## R.med        -2.660e-02  1.545e-02  -1.722 0.085441 .  
## SP           -4.663e-03  5.087e-03  -0.917 0.359580    
## density.corr  9.477e+00  7.816e+00   1.213 0.225670    
## density       7.706e+00  1.928e+00   3.997 7.01e-05 ***
## phi.core      5.278e-01  3.846e-02  13.724  < 2e-16 ***
## k.core        1.650e-03  5.704e-05  28.930  < 2e-16 ***
## FaciesF1      1.575e+03  8.931e+02   1.763 0.078259 .  
## FaciesF10     1.574e+03  8.928e+02   1.763 0.078228 .  
## FaciesF2      1.574e+03  8.929e+02   1.763 0.078342 .  
## FaciesF3      1.574e+03  8.928e+02   1.762 0.078376 .  
## FaciesF5      1.574e+03  8.928e+02   1.763 0.078241 .  
## FaciesF7      1.576e+03  8.926e+02   1.766 0.077862 .  
## FaciesF8      1.575e+03  8.928e+02   1.764 0.078160 .  
## FaciesF9      1.574e+03  8.928e+02   1.763 0.078291 .  
## k.corel              NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.034 on 797 degrees of freedom
## Multiple R-squared:  0.9891, Adjusted R-squared:  0.9888 
## F-statistic:  3284 on 22 and 797 DF,  p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(m3)

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-16
par(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="gold", 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="red", 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="gold", 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="red", 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("gold", "red"))

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] 1
mspe.m3 <- sqrt(sum((k.core.pred3 - data$k.core)^2)/nrow(data))
mspe.m3
## [1] 1.560247e-12
par(mfrow=c(1,1))
plot(y=y<-(data$depth),ylim=rev(range(data$depth)),x=x<-(data$k.core),type="p", col="gold", 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="red", 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("gold", "red"))

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="gold", 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="red", 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("gold", "blue", "green", "red"))