-
Notifications
You must be signed in to change notification settings - Fork 16
/
coefficients.R
179 lines (173 loc) · 6.17 KB
/
coefficients.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
#' Coefficients of a Bayesian Model Average object
#'
#' Extract conditional posterior means and standard deviations, marginal
#' posterior means and standard deviations, posterior probabilities, and
#' marginal inclusions probabilities under Bayesian Model Averaging from an
#' object of class 'bas'
#'
#' Calculates posterior means and (approximate) standard deviations of the
#' regression coefficients under Bayesian Model averaging using g-priors and
#' mixtures of g-priors. Print returns overall summaries. For fully Bayesian
#' methods that place a prior on g, the posterior standard deviations do not
#' take into account full uncertainty regarding g. Will be updated in future
#' releases.
#'
#' @aliases coef.bas coef coefficients coefficients.bas print.coef.bas
#' @param object object of class 'bas' created by BAS
#' @param x object of class 'coef.bas' to print
#' @param n.models Number of top models to report in the printed summary, for
#' coef the default is to use all models. To extract summaries for the Highest
#' Probability Model, use n.models=1 or estimator="HPM".
#' @param estimator return summaries for a selected model, rather than using
#' BMA. Options are 'HPM' (highest posterior probability model) ,'MPM' (median
#' probability model), and 'BMA'
#' @param digits number of significant digits to print
#' @param ... other optional arguments
#' @return \code{coefficients} returns an object of class coef.bas with the
#' following:
#' \item{conditionalmeans}{a matrix with conditional posterior means
#' for each model}
#' \item{conditionalsd}{ standard deviations for each model }
#' \item{postmean}{marginal posterior means of each regression coefficient
#' using BMA}
#' \item{postsd}{marginal posterior standard deviations using BMA}
#' \item{postne0}{vector of posterior inclusion probabilities, marginal
#' probability that a coefficient is non-zero}
#' @note With highly correlated variables, marginal summaries may not be
#' representative of the joint distribution. Use \code{\link{plot.coef.bas}} to
#' view distributions. The value reported for the intercept is
#' under the centered parameterization. Under the Gaussian error
#' model it will be centered at the sample mean of Y.
#' @author Merlise Clyde \email{clyde@@duke.edu}
#' @seealso \code{\link{bas}}, \code{\link{confint.coef.bas}}
#' @references Liang, F., Paulo, R., Molina, G., Clyde, M. and Berger, J.O.
#' (2005) Mixtures of g-priors for Bayesian Variable Selection. Journal of the
#' American Statistical Association. 103:410-423. \cr
#' \doi{10.1198/016214507000001337}
#' @keywords regression
#' @examples
#'
#' data("Hald")
#' hald.gprior = bas.lm(Y~ ., data=Hald, n.models=2^4, alpha=13,
#' prior="ZS-null", initprobs="Uniform", update=10)
#' coef.hald.gprior = coefficients(hald.gprior)
#' coef.hald.gprior
#' plot(coef.hald.gprior)
#' confint(coef.hald.gprior)
#'
#' #Estimation under Median Probability Model
#' coef.hald.gprior = coefficients(hald.gprior, estimator="MPM")
#' coef.hald.gprior
#' plot(coef.hald.gprior)
#' plot(confint(coef.hald.gprior))
#'
#'
#' coef.hald.gprior = coefficients(hald.gprior, estimator="HPM")
#' coef.hald.gprior
#' plot(coef.hald.gprior)
#' confint(coef.hald.gprior)
#'
#' # To add estimation under Best Predictive Model
#'
#'
#' @rdname coef
#' @family bas methods
#' @export
coef.bas <- function(object, n.models, estimator = "BMA", ...) {
if (estimator == "BPM") {
stop("Extracting coefficients for the BPM is not implemented yet")
}
if (estimator == "MPM") {
nvar <- object$n.vars - 1
bestmodel <- (0:nvar)[object$probne0 > .5]
best <- 1
models <- rep(0, nvar + 1)
models[bestmodel + 1] <- 1
if (sum(models) > 1) {
object <- bas.lm(
eval(object$call$formula),
data = eval(object$call$data),
weights = eval(object$call$weights),
n.models = 1,
alpha = object$g,
initprobs = object$probne0,
prior = object$prior,
modelprior = object$modelprior,
update = NULL,
bestmodel = models,
prob.local = 0.0
)
}
}
postprobs <- object$postprobs
if (estimator == "MPM" | estimator == "HPM") {
n.models <- 1
}
if (missing(n.models)) {
n.models <- length(postprobs)
}
topm <- order(-postprobs)[1:n.models]
postprobs <- postprobs[topm] / sum(postprobs[topm])
shrinkage <- object$shrinkage[topm]
conditionalmeans <- list2matrix.bas(object, "mle")[topm, , drop = F]
conditionalmeans[, -1] <- sweep(conditionalmeans[, -1, drop = F], 1,
shrinkage,
FUN = "*"
)
postmean <- as.vector(postprobs %*% conditionalmeans)
conditionalsd <- list2matrix.bas(object, "mle.se")[topm, , drop = F]
if (!(object$prior == "AIC" || object$prior == "BIC")) {
conditionalsd[, -1] <- sweep(conditionalsd[, -1, drop = F], 1,
sqrt(shrinkage),
FUN = "*"
)
}
postsd <- sqrt(postprobs %*% conditionalsd^2 +
postprobs %*% ((sweep(conditionalmeans, 2, postmean, FUN = "-")
)^2))
postsd <- as.vector(postsd)
if (is.null(object$df[topm])) {
df <- rep(object$n, length(postprobs))
if (object$prior == "BIC" | object$prior == "AIC") {
df <- df - object$rank
} else {
df <- df - 1
}
} else {
df <- object$df[topm]
}
out <- list(
postmean = postmean,
postsd = postsd,
probne0 = object$probne0,
conditionalmeans = conditionalmeans,
conditionalsd = conditionalsd,
namesx = object$namesx,
postprobs = postprobs,
n.vars = object$n.vars,
n.models = n.models,
df = df,
estimator = estimator
)
class(out) <- "coef.bas"
return(out)
}
#' Print coefficients generated from coef.bas
#' @rdname coef
#' @aliases print.coef.bas
#' @family bas coefs
#' @method print coef.bas
#' @export
print.coef.bas <- function(x,
digits = max(3, getOption("digits") - 3), ...) {
out <- cbind(x$postmean, x$postsd, x$probne0)
dimnames(out) <- list(x$namesx, c("post mean", "post SD", "post p(B != 0)"))
cat("\n Marginal Posterior Summaries of Coefficients: \n")
cat("\n Using ", x$estimator, "\n")
cat("\n Based on the top ", x$n.models, "models \n")
print.default(format(out, digits = digits),
print.gap = 2,
quote = FALSE, ...
)
invisible()
}