## funktion sx() returns a list
## which is then processed by function
## bayesx.construct to build the
## BayesX model term structure
sx(x)
bayesx.construct(sx(x))
bayesx.construct(sx(x, bs = "rw1"))
bayesx.construct(sx(x, bs = "factor"))
bayesx.construct(sx(x, bs = "offset"))
bayesx.construct(sx(x, z, bs = "te"))
## varying coefficients
bayesx.construct(sx(x1, by = x2))
bayesx.construct(sx(x1, by = x2, center = TRUE))
## using a map for markov random fields
data("FantasyBnd")
plot(FantasyBnd)
bayesx.construct(sx(id, bs = "mrf", map = FantasyBnd))
## random effects
bayesx.construct(sx(id, bs = "re"))
## examples using optional controlling
## parameters and objects
bayesx.construct(sx(x, bs = "ps", knots = 20))
bayesx.construct(sx(x, bs = "ps", nrknots = 20))
bayesx.construct(sx(x, bs = "ps", knots = 20, nocenter = TRUE))
## use of bs with original
## BayesX syntax
bayesx.construct(sx(x, bs = "psplinerw1"))
bayesx.construct(sx(x, bs = "psplinerw2"))
bayesx.construct(sx(x, z, bs = "pspline2dimrw2"))
bayesx.construct(sx(id, bs = "spatial", map = FantasyBnd))
bayesx.construct(sx(x, z, bs = "kriging"))
bayesx.construct(sx(id, bs = "geospline", map = FantasyBnd, nrknots = 5))
bayesx.construct(sx(x, bs = "catspecific"))
if (FALSE) {
## 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)
## increase inner knots
## decrease degree of the P-spline
b2 <- bayesx(y ~ sx(x, knots = 30, degree = 2), method = "REML", data = dat)
## compare reported output
summary(c(b1, b2))
## plot the effect for both models
plot(c(b1, b2), residuals = TRUE)
## 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 model
b <- bayesx(y ~ sx(x) + sx(z, w, bs = "te") + fac,
data = dat, method = "MCMC")
summary(b)
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
names(MunichBnd) <- 1:N
## 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
b <- bayesx(y ~ sx(x1) + sx(id, bs = "mrf", map = MunichBnd),
method = "REML", data = dat)
## summary statistics
summary(b)
## plot the effects
op <- par(no.readonly = TRUE)
par(mfrow = c(1,2))
plot(b, term = "sx(id)", map = MunichBnd,
main = "bayesx() estimate")
plotmap(MunichBnd, x = dat$sp, id = dat$id,
main = "Truth")
par(op)
## 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, bs = "psplinerw1") + 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))
}
Run the code above in your browser using DataLab