# NOT RUN {
# Descriptive Statistics
summary(WVwells)
# Record means of predictors:
# These are used BOTH to eliminate the intercept and to recover predictions later.
mean.logGas<-mean(WVwells$logGProd2012);mean.logGas
mean.logElevation<-mean(WVwells$logElevation);mean.logElevation
mean.RockPres<-mean(WVwells$RockPres);mean.RockPres
# Outcome Variable, De-Meaned
logGas <- matrix(WVwells$logGProd2012-mean.logGas,ncol=1)
# Explanatory Variable Matrix: DE-MEANED PREDICTORS AND NO CONSTANT TERM
# Because we deducted the mean from all predictors and the outcome,
# it is valid to do regression through the origin.
wv.pred <-cbind(WVwells$logElevation-mean.logElevation,WVwells$RockPres-mean.RockPres)
dimnames(wv.pred)[[2]] <- c("LogElevation", "RockPressure")
# OLS Model
fracking.ols<-lm(logGas~wv.pred-1)
summary(fracking.ols)
# Set Number of Iterations:
# WARNING: 50 iterations is intensive on many machines.
# This example was tuned on Amazon Web Services (EC2) over many hours
# with 5,000 iterations--unsuitable in 2020 for most desktop machines.
M<-50#00 #The comment character constrains this to 50 iterations.
set.seed(1000,kind="Mersenne-Twister")#SET SEED FOR CONSISTENCY
# }
# NOT RUN {
# Trial Run, Linear Model of Ideology with Geospatial Errors Using Metropolis-Hastings:
wv.mat <- metropolis.krige(y=logGas, X=wv.pred, east=WVwells$UTMESrf, north=WVwells$UTMNSrf,
powered.exp=0.5, mcmc.samples=M, spatial.share=0.60, range.share=0.31, beta.var=1000,
range.tol=0.1, b.tune=1, nugget.tune=1, psill.tune=30)
# Discard first 20% of Iterations as Burn-In (User Discretion Advised).
wv.mat <- wv.mat[(ceiling(0.2*M)+1):M,]
# Summarize Results
apply(wv.mat,2,quantile,c(0.5,0.05,0.95))
# Convergence Diagnostics
geweke(wv.mat)
heidel.welch(wv.mat)
# Draw Semivariogram
graph.terms<-apply(wv.mat[,1:3],2,quantile,0.5)
raw.semivar<-semivariogram(x=logGas,east=WVwells$UTMESrf,north=WVwells$UTMNSrf)
resid.semivar<-semivariogram(x=fracking.ols$residuals,east=WVwells$UTMESrf,
north=WVwells$UTMNSrf,draw.plot=FALSE)
points(resid.semivar,pch=3,col='blue')
lines(exponential.semivariogram(nugget=graph.terms[1],decay=graph.terms[2],
partial.sill=graph.terms[3],
distance=as.numeric(names(resid.semivar)),power=0.5),col='red')
# Predictive Data for Two Wells Tapped in 2013
well.1<-c(log(1110)-mean.logElevation,1020-mean.RockPres)
well.2<-c(log(643)-mean.logElevation,630-mean.RockPres)
wells.2013<-rbind(well.1,well.2)
site.1<-c(557306.0, 4345265)
site.2<-c(434515.7, 4258449)
well.locations<-rbind(site.1,site.2)
colnames(well.locations)<-c("eastings","northings")
# Make predictions from median parameter values:
median.pred<-krige.pred(pred.x=wells.2013,
pred.east=well.locations[,"eastings"],pred.north=well.locations[,"northings"],
train.y=logGas,train.x=wv.pred,train.east=WVwells$UTMESrf,
train.north=WVwells$UTMNSrf,mcmc.iter=wv.mat)
# Prediction of deviation on logged scale:
median.pred
# Prediction in thousands of cubic feet (MCF):
round(exp(median.pred+mean.logGas))
# Make predictions with 90% credible intervals:
cred.pred<-krige.pred(pred.x=wells.2013,
pred.east=well.locations[,"eastings"],pred.north=well.locations[,"northings"],
train.y=logGas,train.x=wv.pred,train.east=WVwells$UTMESrf,
train.north=WVwells$UTMNSrf,mcmc.iter=wv.mat,credible=0.90)
# Prediction of deviation on logged scale:
cred.pred
# Prediction in thousands of cubic feet (MCF) and the true yield in 2013:
Actual.Yield<-c(471171, 7211)
round(cbind(exp(cred.pred+mean.logGas),Actual.Yield))
# }
Run the code above in your browser using DataLab