# NOT RUN {
ny <- NY_subset
#data cleaning
ny$cathOrth<-ny$catholic+ny$orthodox
ny$consRelig<-ny$mormon+ny$evangelical
ny$jewMus<-ny$jewish+ny$islam
# Explanatory Variable Matrix
psrm.data <-cbind(ny$age, ny$educ, I(ny$age*ny$educ), as.numeric(ny$race==2),
as.numeric(ny$race==3), ny$female, I(as.numeric(ny$race==2)*ny$female),
I(as.numeric(ny$race==3)*ny$female), ny$cathOrth, ny$consRelig,
ny$jewMus, ny$mainline, ny$rural, ny$ownership,
as.numeric(ny$empstat==2), as.numeric(ny$empstat==3),ny$inc14)
dimnames(psrm.data)[[2]] <- c("Age", "Education", "Age.education",
"African.American", "Nonwhite.nonblack","Female",
"African.American.female", "Nonwhite.nonblack.female",
"Catholic.Orthodox", "Evang.Mormon", "Jewish.Muslim",
"Mainline","Ruralism", "Homeowner", "Unemployed",
"Not.in.workforce","Income")
# Outcome Variable
ideo <- matrix(ny$ideology,ncol=1)
# Set Number of Iterations:
# WARNING: 20 iterations is intensive on many machines.
# This example was tuned on Amazon Web Services (EC2) over many hours
# with 20,000 iterations--unsuitable in 2020 for most desktop machines.
#M<-20000
M<-100
set.seed(1,kind="Mersenne-Twister")
# Estimate the Model
ny.fit <- metropolis.krige(formula = ideo ~ psrm.data, coords = cbind(ny$eastings, ny$northings),
powered.exp=1, n.iter=M, spatial.share=0.31,range.share=0.23,beta.var=10,
range.tol=0.01, b.tune=0.1, nugget.tune=20, psill.tune=5)
# Discard first 20% of Iterations as Burn-In (User Discretion Advised).
ny.fit <- burnin(ny.fit, M/5)
# Summarize Results
summary(ny.fit)
#Convergence Diagnostics: Geweke and Heidelberger-Welch
geweke(ny.fit)
heidel.welch(ny.fit)
# Draw Semivariogram
semivariogram(ny.fit)
# }
Run the code above in your browser using DataLab