bayesglm (formula, family = gaussian, data,
weights, subset, na.action,
start = NULL, etastart, mustart,
offset, control = glm.control(...),
model = TRUE, method = "glm.fit",
x = FALSE, y = TRUE, contrasts = NULL,
prior.mean = 0, prior.scale = 2.5,
prior.scale.for.intercept = 10, prior.df = 1,
scaled = TRUE, n.iter = 100, ...)
bayesglm.fit (x, y, weights = rep(1, nobs),
start = NULL, etastart = NULL,
mustart = NULL, offset = rep(0, nobs), family = gaussian(),
control = glm.control(), intercept = TRUE,
prior.mean=0, prior.scale=2.5, prior.scale.for.intercept=10,
prior.df=1, scaled=TRUE)as.data.frame to a data frame) containing
the variables in the model. If not found in data, the
variablNULL or a numeric vector.NAs. The default is set by
the na.action setting of options, and is
NULL or a numeric vector of
length either one or equal to the number of cases.
One or glm.control
for details."glm.fit" uses iteratively reweighted
least squares (IWLS). The only current alternative is
"model.frame" which returns the model frame and does no figlm:
logical values indicating whether the response vector and model
matrix used in the fitting process should be returned as components
of the returned value.
For glm.fit: x is a designcontrasts.arg
of model.matrix.default.glm for details.glm() that uses an approximate EM
algorithm to update the betas at each step using an augmented regression
to represent the prior information.
We use Student-t prior distributions for the coefficients. The prior
distribution for the constant term is set so it applies to the value
when all predictors are set to their mean values.
If scaled=TRUE, the scales for the prior distributions of the
coefficients are determined as follows: For a predictor with only one
value, we just use prior.scale. For a predictor with two values, we use
prior.scale/range(x). For a predictor with more than two values, we use
prior.scale/(2*sd(x)).
We include all the glm() arguments but we haven't tested that all the
options (e.g., offests, contrasts,
deviance for the null model) all work.
The new arguments here are: prior.mean, prior.scale,
prior.scale.for.intercept, prior.df, and
scaled.glm,
bayespolrn <- 100
x1 <- rnorm (n)
x2 <- rbinom (n, 1, .5)
b0 <- 1
b1 <- 1.5
b2 <- 2
y <- rbinom (n, 1, invlogit(b0+b1*x1+b2*x2))
M1 <- glm (y ~ x1 + x2, family=binomial(link="logit"))
display (M1) # (using the display() function from regression.R)
M2 <- bayesglm (y ~ x1 + x2, family=binomial(link="logit"),
prior.scale=Inf, prior.df=Inf)
display (M2) # just a test: this should be identical to classical logit
M3 <- bayesglm (y ~ x1 + x2, family=binomial(link="logit"))
# default Cauchy prior with scale 2.5
display (M3)
M4 <- bayesglm (y ~ x1 + x2, family=binomial(link="logit"),
prior.scale=2.5, prior.df=1)
# Same as M3, explicitly specifying Cauchy prior with scale 2.5
display (M4)
M5 <- bayesglm (y ~ x1 + x2, family=binomial(link="logit"),
prior.scale=2.5, prior.df=7) # t_7 prior with scale 2.5
display (M5)
M6 <- bayesglm (y ~ x1 + x2, family=binomial(link="logit"),
prior.scale=2.5, prior.df=Inf) # normal prior with scale 2.5
display (M6)
# Create separation: set y=1 whenever x2=1
# Now it should blow up without the prior!
y <- ifelse (x2==1, 1, y)
M1 <- glm (y ~ x1 + x2, family=binomial(link="logit"))
display (M1)
M2 <- bayesglm (y ~ x1 + x2, family=binomial(link="logit"),
prior.scale=Inf, prior.df=Inf) # Same as M1
display (M2)
M3 <- bayesglm (y ~ x1 + x2, family=binomial(link="logit"))
display (M3)
M4 <- bayesglm (y ~ x1 + x2, family=binomial(link="logit"),
prior.scale=2.5, prior.df=1) # Same as M3
display (M4)
M5 <- bayesglm (y ~ x1 + x2, family=binomial(link="logit"),
prior.scale=2.5, prior.df=7)
display (M5)
M6 <- bayesglm (y ~ x1 + x2, family=binomial(link="logit"),
prior.scale=2.5, prior.df=Inf)
display (M6)
# bayesglm with gaussian family (bayes lm)
sigma <- 5
y2 <- rnorm (n, b0+b1*x1+b2*x2, sigma)
M7 <- bayesglm (y2 ~ x1 + x2, prior.scale=Inf, prior.df=Inf)
display (M7)Run the code above in your browser using DataLab