# NOT RUN {
# }
# NOT RUN {
nyc <- NYcity_subset
#data cleaning
nyc$cathOrth<-nyc$catholic+nyc$orthodox
nyc$consRelig<-nyc$mormon+nyc$evangelical
nyc$jewMus<-nyc$jewish+nyc$islam
# Explanatory Variable Matrix
psrm.data <-cbind(nyc$age, nyc$educ, I(nyc$age*nyc$educ), as.numeric(nyc$race==2),
as.numeric(nyc$race==3), nyc$female, I(as.numeric(nyc$race==2)*nyc$female),
I(as.numeric(nyc$race==3)*nyc$female), nyc$cathOrth, nyc$consRelig,
nyc$jewMus, nyc$mainline, nyc$rural, nyc$ownership,
as.numeric(nyc$empstat==2), as.numeric(nyc$empstat==3),nyc$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(nyc$ideology,ncol=1)
# WARNING: This example was tuned on Amazon Web Services (EC2) over many hours
# with 150,000 iterations--a strain in 2020 for most desktop machines.
# A test with few iterations allows illustration.
#M<-150000
M<-150
set.seed(1,kind="Mersenne-Twister")
# Estimate the Model
nyc.fit <- metropolis.krige(formula = ideo ~ psrm.data, coords = cbind(nyc$eastings, nyc$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).
nyc.fit <- burnin(nyc.fit, M/5)
# Summarize Results
summary(nyc.fit)
#Convergence Diagnostics: Geweke and Heidelberger-Welch
geweke(nyc.fit)
heidel.welch(nyc.fit)
# Draw Semivariogram
semivariogram(nyc.fit)
# }
Run the code above in your browser using DataLab