## generate some data
set.seed(111)
n <- 200
## regressor
dat <- data.frame(x = runif(n, -3, 3))
## response
dat$y <- with(dat, 1.5 + sin(x) + rnorm(n, sd = 0.6))
## estimate models with
## bayesx REML and MCMC
b1 <- bayesx(y ~ sx(x), method = "REML", data = dat)
## same using mgcv syntax
b1 <- bayesx(y ~ s(x, bs = "ps", k = 20), method = "REML", data = dat)
## now with MCMC
b2 <- bayesx(y ~ sx(x), method = "MCMC",
iter = 1200, burnin = 200, data = dat)
## compare reported output
summary(c(b1, b2))
## plot the effect for both models
plot(c(b1, b2), residuals = TRUE)
## use confint
confint(b1, level = 0.99)
confint(b2, level = 0.99)
## Not run:
# ## more examples
# set.seed(111)
# n <- 500
#
# ## regressors
# dat <- data.frame(x = runif(n, -3, 3), z = runif(n, -3, 3),
# w = runif(n, 0, 6), fac = factor(rep(1:10, n/10)))
#
# ## response
# dat$y <- with(dat, 1.5 + sin(x) + cos(z) * sin(w) +
# c(2.67, 5, 6, 3, 4, 2, 6, 7, 9, 7.5)[fac] + rnorm(n, sd = 0.6))
#
# ## estimate models with
# ## bayesx MCMC and REML
# ## and compare with
# ## mgcv gam()
# b1 <- bayesx(y ~ sx(x) + sx(z, w, bs = "te") + fac,
# data = dat, method = "MCMC")
# b2 <- bayesx(y ~ sx(x) + sx(z, w, bs = "te") + fac,
# data = dat, method = "REML")
# b3 <- gam(y ~ s(x, bs = "ps") + te(z, w, bs = "ps") + fac,
# data = dat)
#
# ## summary statistics
# summary(b1)
# summary(b2)
# summary(b3)
#
# ## plot the effects
# op <- par(no.readonly = TRUE)
# par(mfrow = c(3, 2))
# plot(b1, term = "sx(x)")
# plot(b1, term = "sx(z,w)")
# plot(b2, term = "sx(x)")
# plot(b2, term = "sx(z,w)")
# plot(b3, select = 1)
# vis.gam(b3, c("z","w"), theta = 40, phi = 40)
# par(op)
#
# ## combine models b1 and b2
# b <- c(b1, b2)
#
# ## summary
# summary(b)
#
# ## only plot effect 2 of both models
# plot(b, term = "sx(z,w)")
#
# ## with residuals
# plot(b, term = "sx(z,w)", residuals = TRUE)
#
# ## same model with kriging
# b <- bayesx(y ~ sx(x) + sx(z, w, bs = "kr") + fac,
# method = "REML", data = dat)
# plot(b)
#
#
# ## now a mrf example
# ## note: the regional identification
# ## covariate and the map regionnames
# ## should be coded as integer
# set.seed(333)
#
# ## simulate some geographical data
# data("MunichBnd")
# N <- length(MunichBnd); n <- N*5
#
# ## regressors
# dat <- data.frame(x1 = runif(n, -3, 3),
# id = as.factor(rep(names(MunichBnd), length.out = n)))
# dat$sp <- with(dat, sort(runif(N, -2, 2), decreasing = TRUE)[id])
#
# ## response
# dat$y <- with(dat, 1.5 + sin(x1) + sp + rnorm(n, sd = 1.2))
#
# ## estimate models with
# ## bayesx MCMC and REML
# b1 <- bayesx(y ~ sx(x1) + sx(id, bs = "mrf", map = MunichBnd),
# method = "MCMC", data = dat)
# b2 <- bayesx(y ~ sx(x1) + sx(id, bs = "mrf", map = MunichBnd),
# method = "REML", data = dat)
#
# ## summary statistics
# summary(b1)
# summary(b2)
#
# ## plot the spatial effects
# plot(b1, term = "sx(id)", map = MunichBnd,
# main = "bayesx() MCMC estimate")
# plot(b2, term = "sx(id)", map = MunichBnd,
# main = "bayesx() REML estimate")
# plotmap(MunichBnd, x = dat$sp, id = dat$id,
# main = "Truth")
#
# ## try geosplines instead
# b <- bayesx(y ~ sx(id, bs = "gs", map = MunichBnd) + sx(x1), data = dat)
# summary(b)
# plot(b, term = "sx(id)", map = MunichBnd)
#
# ## geokriging
# b <- bayesx(y ~ sx(id, bs = "gk", map = MunichBnd) + sx(x1),
# method = "REML", data = dat)
# summary(b)
# plot(b, term = "sx(id)", map = MunichBnd)
#
# ## perspective plot of the effect
# plot(b, term = "sx(id)")
#
# ## image and contour plot
# plot(b, term = "sx(id)", image = TRUE,
# contour = TRUE, grid = 200)
#
#
# ## model with random effects
# set.seed(333)
# N <- 30
# n <- N*10
#
# ## regressors
# dat <- data.frame(id = sort(rep(1:N, n/N)), x1 = runif(n, -3, 3))
# dat$re <- with(dat, rnorm(N, sd = 0.6)[id])
#
# ## response
# dat$y <- with(dat, 1.5 + sin(x1) + re + rnorm(n, sd = 0.6))
#
# ## estimate model
# b <- bayesx(y ~ sx(x1) + sx(id, bs = "re"), data = dat)
# summary(b)
# plot(b)
#
# ## extract estimated random effects
# ## and compare with true effects
# plot(fitted(b, term = "sx(id)")$Mean ~ unique(dat$re))
#
#
# ## now a spatial example
# ## with structured and
# ## unstructered spatial
# ## effect
# set.seed(333)
#
# ## simulate some geographical data
# data("MunichBnd")
# N <- length(MunichBnd); names(MunichBnd) <- 1:N
# n <- N*5
#
# ## regressors
# dat <- data.frame(id = rep(1:N, n/N), x1 = runif(n, -3, 3))
# dat$sp <- with(dat, sort(runif(N, -2, 2), decreasing = TRUE)[id])
# dat$re <- with(dat, rnorm(N, sd = 0.6)[id])
#
# ## response
# dat$y <- with(dat, 1.5 + sin(x1) + sp + re + rnorm(n, sd = 0.6))
#
# ## estimate model
# b <- bayesx(y ~ sx(x1) +
# sx(id, bs = "mrf", map = MunichBnd) +
# sx(id, bs = "re"), method = "MCMC", data = dat)
# summary(b)
#
# ## plot all spatial effects
# plot(b, term = "sx(id):mrf", map = MunichBnd,
# main = "Structured spatial effect")
# plot(b, term = "sx(id):re", map = MunichBnd,
# main = "Unstructured spatial effect")
# plot(b, term = "sx(id):total", map = MunichBnd,
# main = "Total spatial effect", digits = 4)
#
#
# ## some experiments with the
# ## stepwise algorithm
# ## generate some data
# set.seed(321)
# n <- 1000
#
# ## regressors
# dat <- data.frame(x1 = runif(n, -3, 3), x2 = runif(n),
# x3 = runif(n, 3, 6), x4 = runif(n, 0, 1))
#
# ## response
# dat$y <- with(dat, 1.5 + sin(x1) + 0.6 * x2 + rnorm(n, sd = 0.6))
#
# ## estimate model with STEP
# b <- bayesx(y ~ sx(x1) + sx(x2) + sx(x3) + sx(x4),
# method = "STEP", algorithm = "cdescent1", CI = "MCMCselect",
# iter = 10000, step = 10, data = dat)
# summary(b)
# plot(b)
#
#
# ## a probit example
# set.seed(111)
# n <- 1000
# dat <- data.frame(x <- runif(n, -3, 3))
#
# dat$z <- with(dat, sin(x) + rnorm(n))
# dat$y <- rep(0, n)
# dat$y[dat$z > 0] <- 1
#
# b <- bayesx(y ~ sx(x), family = "binomialprobit", data = dat)
# summary(b)
# plot(b)
#
#
# ## estimate varying coefficient models
# set.seed(333)
# n <- 1000
# dat <- data.frame(x = runif(n, -3, 3), id = factor(rep(1:4, n/4)))
#
# ## response
# dat$y <- with(dat, 1.5 + sin(x) * c(-1, 0.2, 1, 5)[id] + rnorm(n, sd = 0.6))
#
# ## estimate model
# b <- bayesx(y ~ sx(x, by = id, center = TRUE),
# method = "REML", data = dat)
# summary(b)
# plot(b, resid = TRUE, cex.resid = 0.1)
# ## End(Not run)
Run the code above in your browser using DataLab