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
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
|
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/brglmControl.R, R/zzz_conventions.R
\name{brglmControl}
\alias{brglmControl}
\alias{brglm_control}
\title{Auxiliary function for \code{\link[=glm]{glm()}} fitting using the \code{\link[=brglmFit]{brglmFit()}}
method.}
\usage{
brglmControl(
epsilon = 1e-06,
maxit = 100,
check_aliasing = TRUE,
trace = FALSE,
type = c("AS_mixed", "AS_mean", "AS_median", "correction", "MPL_Jeffreys", "ML"),
transformation = "identity",
slowit = 1,
response_adjustment = NULL,
max_step_factor = 12,
a = 1/2,
...
)
brglm_control(
epsilon = 1e-06,
maxit = 100,
check_aliasing = TRUE,
trace = FALSE,
type = c("AS_mixed", "AS_mean", "AS_median", "correction", "MPL_Jeffreys", "ML"),
transformation = "identity",
slowit = 1,
response_adjustment = NULL,
max_step_factor = 12,
a = 1/2,
...
)
}
\arguments{
\item{epsilon}{positive convergence tolerance epsilon. Default is
\code{1e-06}.}
\item{maxit}{integer giving the maximal number of iterations
allowed. Default is \code{100}.}
\item{check_aliasing}{logical indicating where a QR decomposition
of the model matrix should be used to check for
aliasing. Default is \code{TRUE}. See Details.}
\item{trace}{logical indicating if output should be produced for
each iteration. Default is \code{FALSE}.}
\item{type}{the type of fitting method to be used. The options are
\code{"AS_mean"} (mean-bias reducing adjusted scores), \code{"AS_median"}
(median-bias reducing adjusted scores), \code{"AS_mixed"} (bias
reduction using mixed score adjustments; default),
\code{"correction"} (asymptotic bias correction), \code{"MPL_Jeffreys"}
(maximum penalized likelihood with powers of the Jeffreys prior
as penalty) and \code{"ML"} (maximum likelihood).}
\item{transformation}{the transformation of the dispersion to be
estimated. Default is \code{"identity"}. See Details.}
\item{slowit}{a positive real used as a multiplier for the
stepsize. The smaller it is the smaller the steps are. Default
is \code{1}.}
\item{response_adjustment}{a (small) positive constant or a vector
of such. Default is \code{NULL}. See Details.}
\item{max_step_factor}{the maximum number of step halving steps to
consider. Default is \code{12}.}
\item{a}{power of the Jeffreys prior penalty. See Details.}
\item{...}{further arguments passed to \code{\link[=brglmControl]{brglmControl()}}. Currently
ignored in the output.}
}
\value{
a list with components named as the arguments, including
symbolic expressions for the dispersion transformation
(\code{Trans}) and its inverse (\code{inverseTrans})
}
\description{
Typically only used internally by \code{\link[=brglmFit]{brglmFit()}}, but may be used to
construct a \code{control} argument.
}
\details{
\code{\link[=brglmControl]{brglmControl()}} provides default values and sanity checking for
the various constants that control the iteration and generally the
behaviour of \code{\link[=brglmFit]{brglmFit()}}.
When \code{trace = TRUE}, calls to \code{\link[=cat]{cat()}} produce the output for each
iteration. Hence, \verb{options(digits = *)} can be used to increase
the precision.
When \code{check_aliasing = TRUE} (default), a QR decomposition of the
model matrix is computed to check for aliasing. If the model matrix
is known to be of full rank, then \code{check_aliasing = FALSE} avoids
the extra computational overhead of an additional QR decomposition,
which can be substantial for large model matrices. However, setting
\code{check_aliasing = FALSE} tells \code{\link[=brglmFit]{brglmFit()}} that the model matrix
is full rank, and hard to trace back errors will result if it is
rank deficient.
\code{transformation} sets the transformation of the dispersion
parameter for which the bias reduced estimates are computed. Can be
one of \code{"identity"}, \code{"sqrt"}, \code{"inverse"}, \code{"log"} and
\code{"inverseSqrt"}. Custom transformations are accommodated by
supplying a list of two expressions (transformation and inverse
transformation). See the examples for more details.
The value of \code{response_adjustment} is only relevant if \code{\link[=brglmFit]{brglmFit()}}
is called with \code{start = NULL}, and \code{family} is \code{\link[=binomial]{binomial()}} or
\code{\link[=poisson]{poisson()}}. For those models, an initial maximum likelihood fit is
obtained on adjusted data to provide starting values for the
iteration in \code{\link[=brglmFit]{brglmFit()}}. The value of \code{response_adjustment}
governs how the data is adjusted. Specifically, if \code{family} is
\code{\link[=binomial]{binomial()}}, then the responses and totals are adjusted by
\code{response_adjustment} and \code{2 * response_adjustment}, respectively;
if \code{family} is \code{\link[=poisson]{poisson()}}, then the responses are adjusted by and
\code{response_adjustment}. \code{response_adjustment = NULL} (default)
is equivalent to setting it to
\code{"number of parameters" / "number of observations"}.
When \code{type = "AS_mixed"} (default), mean bias reduction is used for
the regression parameters, and median bias reduction for the
dispersion parameter, if that is not fixed. This adjustment has
been developed based on equivariance arguments (see, Kosmidis et
al, 2020, Section 4) in order to produce regression parameter
estimates that are invariant to arbitrary contrasts, and estimates
for the dispersion parameter that are invariant to arbitrary
non-linear transformations. \code{type = "AS_mixed"} and \code{type = "AS_mean"} return the same results if \code{\link[=brglmFit]{brglmFit()}} is called with
\code{family} \code{\link[=binomial]{binomial()}} or \code{\link[=poisson]{poisson()}} (i.e. families with fixed
dispersion).
When \code{type = "MPL_Jeffreys"}, \code{\link[=brglmFit]{brglmFit()}} will maximize the
penalized log-likelihood \deqn{l(\beta, \phi) + a\log \det i(\beta,
\phi)}{l(beta, phi) + a log det i(beta, phi)} where \eqn{i(\beta,
\phi)}{i(beta, phi)} is the expected information matrix about the
regression parameters \eqn{\beta} and the dispersion parameter
\eqn{\phi}. See, \code{vignette("iteration", "brglm2")} for more
information. The argument \code{a} controls the amount of penalization
and its default value is \code{a = 1/2}, corresponding to maximum
penalized likelihood using a Jeffreys-prior penalty. See, Kosmidis
& Firth (2021) for proofs and discussion about the finiteness and
shrinkage properties of the maximum penalized likelihood estimators
for binomial-response generalized linear models.
The estimates from \code{type = "AS_mean"} and \code{type = "MPL_Jeffreys"} with \code{a = 1/2} (default) are identical for
Poisson log-linear models and logistic regression models, i.e. for
binomial and Poisson regression models with canonical links. See,
Firth (1993) for details.
\code{\link[=brglm_control]{brglm_control()}} is an alias to \code{\link[=brglmControl]{brglmControl()}}.
}
\examples{
data("coalition", package = "brglm2")
## The maximum likelihood fit with log link
coalitionML <- glm(duration ~ fract + numst2, family = Gamma, data = coalition)
## Bias reduced estimation of the dispersion parameter
coalitionBRi <- glm(duration ~ fract + numst2, family = Gamma, data = coalition,
method = "brglmFit")
coef(coalitionBRi, model = "dispersion")
## Bias reduced estimation of log(dispersion)
coalitionBRl <- glm(duration ~ fract + numst2, family = Gamma, data = coalition,
method = "brglmFit", transformation = "log")
coef(coalitionBRl, model = "dispersion")
## Just for illustration: Bias reduced estimation of dispersion^0.25
my_transformation <- list(expression(dispersion^0.25), expression(transformed_dispersion^4))
coalitionBRc <- update(coalitionBRi, transformation = my_transformation)
coef(coalitionBRc, model = "dispersion")
}
\references{
Kosmidis I, Firth D (2021). Jeffreys-prior penalty, finiteness
and shrinkage in binomial-response generalized linear
models. \emph{Biometrika}, \strong{108}, 71-82. \doi{10.1093/biomet/asaa052}.
Kosmidis I, Kenne Pagui E C, Sartori N (2020). Mean and median bias
reduction in generalized linear models. \emph{Statistics and Computing},
\strong{30}, 43-59. \doi{10.1007/s11222-019-09860-6}.
Firth D (1993). Bias reduction of maximum likelihood estimates.
Biometrika, \strong{80}, 27-38. \doi{10.2307/2336755}.
}
\seealso{
\code{\link[=brglm_fit]{brglm_fit()}} and \code{\link[=glm.fit]{glm.fit()}}
}
\author{
Ioannis Kosmidis \verb{[aut, cre]} \email{ioannis.kosmidis@warwick.ac.uk}
}
|
