11# Part of the rstanarm package for estimating model parameters
22# Copyright (C) 2013, 2014, 2015, 2016, 2017 Trustees of Columbia University
3- #
3+ #
44# This program is free software; you can redistribute it and/or
55# modify it under the terms of the GNU General Public License
66# as published by the Free Software Foundation; either version 3
77# of the License, or (at your option) any later version.
8- #
8+ #
99# This program is distributed in the hope that it will be useful,
1010# but WITHOUT ANY WARRANTY; without even the implied warranty of
1111# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1212# GNU General Public License for more details.
13- #
13+ #
1414# You should have received a copy of the GNU General Public License
1515# along with this program; if not, write to the Free Software
1616# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
1717
1818# ' Bayesian beta regression models via Stan
1919# '
20- # ' \if{html}{\figure{stanlogo.png }{options: width="25" alt="https://mc-stan.org/about/logo/"}}
21- # ' Beta regression modeling with optional prior distributions for the
20+ # ' \if{html}{\figure{logo.svg }{options: width="25" alt="https://mc-stan.org/about/logo/"}}
21+ # ' Beta regression modeling with optional prior distributions for the
2222# ' coefficients, intercept, and auxiliary parameter \code{phi} (if applicable).
2323# '
2424# ' @export
2525# ' @templateVar armRef (Ch. 3-6)
2626# ' @templateVar pkg betareg
2727# ' @templateVar pkgfun betareg
28- # ' @templateVar sameargs model,offset,weights
28+ # ' @templateVar sameargs model,offset,weights
2929# ' @templateVar rareargs na.action
3030# ' @templateVar fun stan_betareg
3131# ' @templateVar fitfun stan_betareg.fit
4343# ' @template args-algorithm
4444# ' @template args-adapt_delta
4545# ' @template args-QR
46- # '
47- # ' @param link Character specification of the link function used in the model
46+ # '
47+ # ' @param link Character specification of the link function used in the model
4848# ' for mu (specified through \code{x}). Currently, "logit", "probit",
4949# ' "cloglog", "cauchit", "log", and "loglog" are supported.
50- # ' @param link.phi If applicable, character specification of the link function
51- # ' used in the model for \code{phi} (specified through \code{z}). Currently,
50+ # ' @param link.phi If applicable, character specification of the link function
51+ # ' used in the model for \code{phi} (specified through \code{z}). Currently,
5252# ' "identity", "log" (default), and "sqrt" are supported. Since the "sqrt"
5353# ' link function is known to be unstable, it is advisable to specify a
5454# ' different link function (or to model \code{phi} as a scalar parameter
5555# ' instead of via a linear predictor by excluding \code{z} from the
5656# ' \code{formula} and excluding \code{link.phi}).
57- # ' @param prior_z Prior distribution for the coefficients in the model for
57+ # ' @param prior_z Prior distribution for the coefficients in the model for
5858# ' \code{phi} (if applicable). Same options as for \code{prior}.
59- # ' @param prior_intercept_z Prior distribution for the intercept in the model
59+ # ' @param prior_intercept_z Prior distribution for the intercept in the model
6060# ' for \code{phi} (if applicable). Same options as for \code{prior_intercept}.
61- # ' @param prior_phi The prior distribution for \code{phi} if it is \emph{not}
62- # ' modeled as a function of predictors. If \code{z} variables are specified
63- # ' then \code{prior_phi} is ignored and \code{prior_intercept_z} and
61+ # ' @param prior_phi The prior distribution for \code{phi} if it is \emph{not}
62+ # ' modeled as a function of predictors. If \code{z} variables are specified
63+ # ' then \code{prior_phi} is ignored and \code{prior_intercept_z} and
6464# ' \code{prior_z} are used to specify the priors on the intercept and
6565# ' coefficients in the model for \code{phi}. When applicable, \code{prior_phi}
6666# ' can be a call to \code{exponential} to use an exponential distribution, or
6767# ' one of \code{normal}, \code{student_t} or \code{cauchy} to use half-normal,
6868# ' half-t, or half-Cauchy prior. See \code{\link{priors}} for details on these
6969# ' functions. To omit a prior ---i.e., to use a flat (improper) uniform
7070# ' prior--- set \code{prior_phi} to \code{NULL}.
71- # '
72- # ' @details The \code{stan_betareg} function is similar in syntax to
73- # ' \code{\link[betareg]{betareg}} but rather than performing maximum
74- # ' likelihood estimation, full Bayesian estimation is performed (if
75- # ' \code{algorithm} is \code{"sampling"}) via MCMC. The Bayesian model adds
71+ # '
72+ # ' @details The \code{stan_betareg} function is similar in syntax to
73+ # ' \code{\link[betareg]{betareg}} but rather than performing maximum
74+ # ' likelihood estimation, full Bayesian estimation is performed (if
75+ # ' \code{algorithm} is \code{"sampling"}) via MCMC. The Bayesian model adds
7676# ' priors (independent by default) on the coefficients of the beta regression
7777# ' model. The \code{stan_betareg} function calls the workhorse
7878# ' \code{stan_betareg.fit} function, but it is also possible to call the
7979# ' latter directly.
80- # '
80+ # '
8181# ' @seealso The vignette for \code{stan_betareg}.
8282# ' \url{https://mc-stan.org/rstanarm/articles/}
83- # '
84- # ' @references Ferrari, SLP and Cribari-Neto, F (2004). Beta regression for
83+ # '
84+ # ' @references Ferrari, SLP and Cribari-Neto, F (2004). Beta regression for
8585# ' modeling rates and proportions. \emph{Journal of Applied Statistics}.
8686# ' 31(7), 799--815.
87- # '
87+ # '
8888# ' @examples
8989# ' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") {
9090# ' ### Simulated data
9696# ' y <- rbeta(N, mu * phi, (1 - mu) * phi)
9797# ' hist(y, col = "dark grey", border = FALSE, xlim = c(0,1))
9898# ' fake_dat <- data.frame(y, x, z)
99- # '
99+ # '
100100# ' fit <- stan_betareg(
101- # ' y ~ x | z, data = fake_dat,
102- # ' link = "logit",
103- # ' link.phi = "log",
101+ # ' y ~ x | z, data = fake_dat,
102+ # ' link = "logit",
103+ # ' link.phi = "log",
104104# ' algorithm = "optimizing" # just for speed of example
105- # ' )
105+ # ' )
106106# ' print(fit, digits = 2)
107107# ' }
108108stan_betareg <-
109- function (formula ,
110- data ,
111- subset ,
112- na.action ,
113- weights ,
114- offset ,
115- link = c(" logit" , " probit" , " cloglog" , " cauchit" , " log" , " loglog" ),
116- link.phi = NULL ,
117- model = TRUE ,
118- y = TRUE ,
119- x = FALSE ,
120- ... ,
121- prior = normal(autoscale = TRUE ),
122- prior_intercept = normal(autoscale = TRUE ),
123- prior_z = normal(autoscale = TRUE ),
124- prior_intercept_z = normal(autoscale = TRUE ),
125- prior_phi = exponential(autoscale = TRUE ),
126- prior_PD = FALSE ,
127- algorithm = c(" sampling" , " optimizing" , " meanfield" , " fullrank" ),
128- adapt_delta = NULL ,
129- QR = FALSE ) {
130-
109+ function (
110+ formula ,
111+ data ,
112+ subset ,
113+ na.action ,
114+ weights ,
115+ offset ,
116+ link = c(" logit" , " probit" , " cloglog" , " cauchit" , " log" , " loglog" ),
117+ link.phi = NULL ,
118+ model = TRUE ,
119+ y = TRUE ,
120+ x = FALSE ,
121+ ... ,
122+ prior = normal(autoscale = TRUE ),
123+ prior_intercept = normal(autoscale = TRUE ),
124+ prior_z = normal(autoscale = TRUE ),
125+ prior_intercept_z = normal(autoscale = TRUE ),
126+ prior_phi = exponential(autoscale = TRUE ),
127+ prior_PD = FALSE ,
128+ algorithm = c(" sampling" , " optimizing" , " meanfield" , " fullrank" ),
129+ adapt_delta = NULL ,
130+ QR = FALSE
131+ ) {
131132 if (! requireNamespace(" betareg" , quietly = TRUE )) {
132133 stop(" Please install the betareg package before using 'stan_betareg'." )
133134 }
134135 if (! has_outcome_variable(formula )) {
135136 stop(" LHS of formula must be specified." )
136137 }
137-
138+
138139 mc <- match.call(expand.dots = FALSE )
139140 data <- validate_data(data , if_missing = environment(formula ))
140141 mc $ data <- data
141142 mc $ model <- mc $ y <- mc $ x <- TRUE
142-
143+
143144 # NULLify any Stan specific arguments in mc
144145 mc $ prior <- mc $ prior_intercept <- mc $ prior_PD <- mc $ algorithm <-
145146 mc $ adapt_delta <- mc $ QR <- mc $ sparse <- mc $ prior_dispersion <- NULL
146-
147+
147148 mc $ drop.unused.levels <- TRUE
148149 mc [[1L ]] <- quote(betareg :: betareg )
149150 mc $ control <- betareg :: betareg.control(maxit = 0 , fsmaxit = 0 )
@@ -155,60 +156,91 @@ stan_betareg <-
155156 Z <- model.matrix(br , model = " precision" )
156157 weights <- validate_weights(as.vector(model.weights(mf )))
157158 offset <- validate_offset(as.vector(model.offset(mf )), y = Y )
158-
159+
159160 # check if user specified matrix for precision model
160- if (length(grep(" \\ |" , all.names(formula ))) == 0 &&
161- is.null(link.phi ))
161+ if (
162+ length(grep(" \\ |" , all.names(formula ))) == 0 &&
163+ is.null(link.phi )
164+ ) {
162165 Z <- NULL
163-
166+ }
167+
164168 algorithm <- match.arg(algorithm )
165169 link <- match.arg(link )
166170 link_phi <- match.arg(link.phi , c(NULL , " log" , " identity" , " sqrt" ))
167-
168- stanfit <-
169- stan_betareg.fit(x = X , y = Y , z = Z ,
170- weights = weights , offset = offset ,
171- link = link , link.phi = link.phi ,
172- ... ,
173- prior = prior , prior_z = prior_z ,
174- prior_intercept = prior_intercept ,
175- prior_intercept_z = prior_intercept_z ,
176- prior_phi = prior_phi , prior_PD = prior_PD ,
177- algorithm = algorithm , adapt_delta = adapt_delta ,
178- QR = QR )
179- if (algorithm != " optimizing" && ! is(stanfit , " stanfit" )) return (stanfit )
180- if (is.null(link.phi ) && is.null(Z ))
171+
172+ stanfit <-
173+ stan_betareg.fit(
174+ x = X ,
175+ y = Y ,
176+ z = Z ,
177+ weights = weights ,
178+ offset = offset ,
179+ link = link ,
180+ link.phi = link.phi ,
181+ ... ,
182+ prior = prior ,
183+ prior_z = prior_z ,
184+ prior_intercept = prior_intercept ,
185+ prior_intercept_z = prior_intercept_z ,
186+ prior_phi = prior_phi ,
187+ prior_PD = prior_PD ,
188+ algorithm = algorithm ,
189+ adapt_delta = adapt_delta ,
190+ QR = QR
191+ )
192+ if (algorithm != " optimizing" && ! is(stanfit , " stanfit" )) {
193+ return (stanfit )
194+ }
195+ if (is.null(link.phi ) && is.null(Z )) {
181196 link_phi <- " identity"
197+ }
182198 sel <- apply(X , 2L , function (x ) ! all(x == 1 ) && length(unique(x )) < 2 )
183- X <- X [ , ! sel , drop = FALSE ]
199+ X <- X [, ! sel , drop = FALSE ]
184200 if (! is.null(Z )) {
185201 sel <- apply(Z , 2L , function (x ) ! all(x == 1 ) && length(unique(x )) < 2 )
186- Z <- Z [ , ! sel , drop = FALSE ]
202+ Z <- Z [, ! sel , drop = FALSE ]
187203 }
188- fit <-
189- nlist(stanfit , algorithm , data , offset , weights ,
190- x = X , y = Y , z = Z %ORifNULL % model.matrix(y ~ 1 ),
191- family = beta_fam(link ), family_phi = beta_phi_fam(link_phi ),
192- formula , model = mf , terms = mt , call = match.call(),
193- na.action = attr(mf , " na.action" ), contrasts = attr(X , " contrasts" ),
194- stan_function = " stan_betareg" )
204+ fit <-
205+ nlist(
206+ stanfit ,
207+ algorithm ,
208+ data ,
209+ offset ,
210+ weights ,
211+ x = X ,
212+ y = Y ,
213+ z = Z %ORifNULL % model.matrix(y ~ 1 ),
214+ family = beta_fam(link ),
215+ family_phi = beta_phi_fam(link_phi ),
216+ formula ,
217+ model = mf ,
218+ terms = mt ,
219+ call = match.call(),
220+ na.action = attr(mf , " na.action" ),
221+ contrasts = attr(X , " contrasts" ),
222+ stan_function = " stan_betareg"
223+ )
195224 out <- stanreg(fit )
196225 if (algorithm == " optimizing" ) {
197226 out $ log_p <- stanfit $ log_p
198227 out $ log_g <- stanfit $ log_g
199228 }
200- out $ xlevels <- lapply(mf [,- 1 ], FUN = function (x ) {
229+ out $ xlevels <- lapply(mf [, - 1 ], FUN = function (x ) {
201230 xlev <- if (is.factor(x ) || is.character(x )) levels(x ) else NULL
202231 xlev [! vapply(xlev , is.null , NA )]
203232 })
204233 out $ levels <- br $ levels
205- if (! x )
234+ if (! x ) {
206235 out $ x <- NULL
207- if (! y )
236+ }
237+ if (! y ) {
208238 out $ y <- NULL
209- if (! model )
239+ }
240+ if (! model ) {
210241 out $ model <- NULL
211-
242+ }
243+
212244 structure(out , class = c(" stanreg" , " betareg" ))
213245 }
214246
@@ -219,21 +251,27 @@ beta_fam <- function(link = "logit") {
219251 if (link == " loglog" ) {
220252 out <- binomial(" cloglog" )
221253 out $ linkinv <- function (eta ) {
222- 1 - pmax(pmin(- expm1(- exp(eta )), 1 - .Machine $ double.eps ),
223- .Machine $ double.eps )
254+ 1 -
255+ pmax(
256+ pmin(- expm1(- exp(eta )), 1 - .Machine $ double.eps ),
257+ .Machine $ double.eps
258+ )
224259 }
225260 out $ linkfun <- function (mu ) log(- log(mu ))
226261 } else {
227262 out <- binomial(link )
228263 }
229264 out $ family <- " beta"
230265 out $ variance <- function (mu , phi ) mu * (1 - mu ) / (phi + 1 )
231- out $ dev.resids <- function (y , mu , wt )
266+ out $ dev.resids <- function (y , mu , wt ) {
232267 stop(" 'dev.resids' function should not be called" )
233- out $ aic <- function (y , n , mu , wt , dev )
268+ }
269+ out $ aic <- function (y , n , mu , wt , dev ) {
234270 stop(" 'aic' function should not have been called" )
235- out $ simulate <- function (object , nsim )
271+ }
272+ out $ simulate <- function (object , nsim ) {
236273 stop(" 'simulate' function should not have been called" )
274+ }
237275 return (out )
238276}
239277
@@ -242,11 +280,14 @@ beta_phi_fam <- function(link = "log") {
242280 out <- poisson(link )
243281 out $ family <- " beta_phi"
244282 out $ variance <- function (mu , phi ) mu * (1 - mu ) / (phi + 1 )
245- out $ dev.resids <- function (y , mu , wt )
283+ out $ dev.resids <- function (y , mu , wt ) {
246284 stop(" 'dev.resids' function should not be called" )
247- out $ aic <- function (y , n , mu , wt , dev )
285+ }
286+ out $ aic <- function (y , n , mu , wt , dev ) {
248287 stop(" 'aic' function should not have been called" )
249- out $ simulate <- function (object , nsim )
288+ }
289+ out $ simulate <- function (object , nsim ) {
250290 stop(" 'simulate' function should not have been called" )
291+ }
251292 return (out )
252293}
0 commit comments