Get the profiles of arbitrary objectives for arbitrary ‘glm’-like models
Calculates the profiles of arbitrary objectives (inference functions
in the terminology of Lindsay and Qu, 2003) for the parameters of
arbitrary glm-like models with linear
predictor. It provides a variety of options such as profiling over a
pre-specified grid, profiling until the profile of the objective
reaches the values of a quantile, calculating the profile traces along
with the profiled objectives, and others.
profileModel(fitted, gridsize = 20, stdn = 5, stepsize = 0.5,
grid.bounds = NULL, quantile = NULL,
objective = stop("'objective' is missing."),
agreement = TRUE, verbose = TRUE, trace.prelim = FALSE,
which = 1:length(coef(fitted)), profTraces = TRUE,
zero.bound = 1e-08, scale = FALSE,
stdErrors = NULL, ...)
prelim.profiling(fitted, quantile = qchisq(0.95, 1),
objective = stop("'objective' is missing."),
verbose = TRUE, which = 1:length(coef(fitted)),
stepsize = 0.5, stdn = 5, agreement = TRUE,
trace.prelim = FALSE,
stdErrors = NULL, ...)
profiling(fitted, grid.bounds, gridsize = 20, verbose = TRUE,
objective = stop("'objective' is missing."),
agreement = TRUE, which = 1:length(coef(fitted)),
profTraces = TRUE, zero.bound = 1e-08, ...)fitted |
a |
which |
which parameters should be profiled? Has to be a vector
of integers for |
grid.bounds |
a matrix of dimension |
gridsize |
The number of equidistant parameter values to be taken
between the values specified in the entries of |
stepsize |
a positive integer that is used in
|
stdn |
in |
quantile |
a quantile, indicating the range that the profile must
cover. The default value in |
objective |
the function to be profiled. It is a function of the
restricted fitted object and other arguments (see
|
agreement |
logical indicating whether the fitting method used
for |
verbose |
logical. If |
trace.prelim |
logical. If |
profTraces |
logical indicating whether the profile traces should
be returned. The default is |
zero.bound |
a small positive constant. The difference of the
objective at the restricted fit from the objective at
|
scale |
logical. The
default is |
stdErrors |
The vector estimated asymptotic standard errors reported
from the fitting procedure. The default is |
... |
further arguments passed to the specified objective. |
mf <- model.frame(fitted$terms,data=eval(fitted$call$data)) ;
model.matrix(fitted$terms,mf,contrasts = fitted$contrasts)
(or just model.matrix(fitted), for fitted objects that
support the model.matrix method.)
Exception to this are objects returned by BTm of the
BradleyTerry package, where some special handling of the required
objects takes place.
Note that any or both of data and contrasts could be
NULL. This depends whether the data argument has been
supplied to the procedure and whether fitted$contrast exists.
If the fitting procedure that resulted fitted supports an
etastart argument (see glm) and
fitted$linear.predictor contains the estimated linear
predictors then during profiling, the appropriate starting values
are supplied to the fitting procedure. In this way, the iteration is
accelerated and is more stable, numerically. However, it is not necessary
that etastart is supported. In the latter case no starting
values are supplied to the fitting procedure during profiling.
Support for a summary method is
optional. summary is only used for obtaining the
estimated asymptotic standard errors associated to the coefficients in
fitted. If stdErrors=NULL the standard errors are taken
to be summary(fitted)$coefficients[,2] which is the place where
the estimated asymptotic standard errors usually are for
glm-like objects. If this this is not the case then
stdErrors should be set appropriately.
profiling is the workhorse function that does the basic operation of
profiling objectives over a user-specified grid of values. For a given
parameter β, the restricted fit
F(b) is calculated by constraining
β to a point b of the grid. Then the difference
D(F(b)) = P(F(b)) - P(G),
is calculated, where P is the objective specified by the user
and G is the original fit (fitted). For convex
objectives that are minimized at the estimates of G (see
agreement), D(G)=0.
prelim.profiling refers only to convex objectives and searches for
and returns the grid bounds (grid.bounds) for each
profiled parameter that should be used in order the profile to cover
quantile. For a given parameter β,
prelim.profiling also checks whether such enclosure can be
found and returns a logical matrix intersects of dimension
length(which) by 2 that indicates if the profile covers the
quantile to the left and to the right of the estimate in
fitted. At step i of the search a value b_i is
proposed for β and D(F(b_i)) is calculated. If
D(F(b_i))<q, where q is quantile, the next
proposed value is
b_{i+1} = b_i +- (i+1) C min(s, 30) /abs(L) ,
where C is stepsize, s is the
estimated asymptotic standard error of β from G and
L is the slope of the line segment connecting the points
[b_i, D(F(b_i))] and
[b_{i-1},
D(F(b_{i-1}))]. +- is + if the search is on the
right of the estimate of β and - on the left. If an
increase of D is expected then the step slows down. If
abs(L)<1 then abs(L) is set to 1 and if
abs(L)>500 then abs(L) is set to 500. In
this way the iteration is conservative by avoiding very small steps
but not over-conservative by avoiding very large steps.
If the maximum number of steps stdn/stepsize (call this M)
was taken and the quantile was not covered by the profile but the three
last absolute slopes where positive then the iteration is restarted
form b_{M-1} with 2C instead of C in the step
calculation. If the three last slopes were less than 1e-8 in
absolute value then the iteration stops and it is considered that
D has an asymptote at the corresponding direction (left or right).
Note that when the latter takes place the iteration has already moved
6 C min(s, 30) units on the scale of β,
since the first value of b were a slope of 1e-8 in absolute value
was detected. Thus we could safely say that an asymptote has been
detected and avoid calculation of F(beta=b) for
extremely large b's.
Very small values of stepsize make prelim.profiling take
very small steps with the effect of slowing down the execution
time. Large values of stepsize are only recommended when the
estimated asymptotic standard errors are very small in fitted.
profileModel is a wrapper function that collects and combines
the capabilities of profiling and prelim.profiling by
providing a unified interface for their functions, as well as appropriateness
checks on the arguments. When both quantile and
grid.bounds are NULL then profiling is called and
profiling takes place for stdn estimated asymptotic standard
errors on the left and on the right of the estimates in
fitted. This could be used for taking a quick look of the
profiles around the estimate. With only the quantile being
NULL, profiling is performed on the the specified grid of
values. When quantile is specified and grid.bounds is
NULL, prelim.profiling is called and its result is
passed to profiling. If both quantile and
grid.bounds then grid.bounds prevails and profiling is
performed on the specified grid.
profiling returns a list of profiles, with one named component
for each parameter profiled. Each component of
the list contains the profiled parameter values and the corresponding
differences of the objective at the restricted fit from the
objective at fitted. When profTraces=TRUE the corresponding
profile traces are cbind'ed to each component of the
list.
prelim.profiling returns a list with components
intersects and grid.bounds.
profileModel returns an object of class "profileModel"
that has the attribute includes.traces corresponding to the
value of the profTraces argument. The "profileModel"
object is a list of the following components:
profiles |
the result of |
fit |
the |
quantile |
the |
gridsize |
the |
intersects |
if |
profiled.parameters |
a vector of integers indicating which parameters were profiled. |
profiled.objective |
the profiled objective with any additional
arguments passed through |
isNA |
a logical vector indicating which of the parameters in
|
agreement |
the |
zero.bound |
the |
grid.bounds |
the grid bounds that were used for profiling. |
call |
the matched call. |
Methods specific to objects of class "profileModel" are
print, see print.profileModel.
signedSquareRoots, see signedSquareRoots.
profConfint, see profConfint.
plot, see plot.profileModel.
pairs, see pairs.profileModel.
Ioannis Kosmidis <email: ioannis.kosmidis@warwick.ac.uk>
Lindsay, B. G. and Qu, A. (2003). Inference functions and quadratic score tests. Statistical Science 18, 394–410.
Chambers, J. M. and Hastie, T. J. (1992) Statistical Models in S. Chapman \& Hall/CRC.
## Begin Example 1
library(MASS)
m1 <- glm(Claims ~ District + Group + Age + offset(log(Holders)),
data = Insurance, family = poisson)
# profile deviance +-5 estimated standard errors from the estimate
prof0 <- profileModel(m1, objective = "ordinaryDeviance")
# profile deviance over a grid of values
gridd <- rep(c(-1,1), length(coef(m1)))
prof1 <- profileModel(m1, grid.bounds = gridd,
objective = "ordinaryDeviance")
# profile deviance until the profile reaches qchisq(0.95,1)
prof2 <- profileModel(m1, quantile = qchisq(0.95,1) ,
objective = "ordinaryDeviance")
# plot the profiles of the deviance
plot(prof2)
# quite quadratic in shape. Just to make sure:
plot(prof2, signed = TRUE)
# Ok straight lines. So we expect first order asymptotics to work well;
## Not run:
# plot the profiles of the Rao score statistic
# profile Rao's score statistic
prof3 <- update(prof2, objective = "RaoScoreStatistic",
X = model.matrix(m1))
plot(prof3)
# The 95% confidence intervals based on prof2 and prof3 and the simple Wald
# confidence intervals:
profConfint(prof2)
profConfint(prof3)
stdErrors <- coef(summary(m1))[,2]
coef(m1)+ qnorm(0.975) * cbind(-stdErrors,stdErrors)
# They are all quite similar in value. The result of a quadratic likelihood.
## End Example
## End(Not run)
## Begin Example 2: Monotone likelihood; data separation;
library(MASS)
y <- c(0, 0, 1, 0)
tots <- c(2, 2, 5, 2)
x1 <- c(1, 0, 1, 0)
x2 <- c(1, 1, 0, 0)
m2 <- glm(y/tots ~ x1 + x2, weights = tots,
family = binomial)
prof <- profileModel(m2, quantile=qchisq(0.95,1),
objective = "ordinaryDeviance")
plot(prof)
profConfint(prof)
# profile.glm fails to detect the finite endpoints
confint(m2)
## End Example
## Not run:
## Begin Example 3: polr
library(MASS)
options(contrasts = c("contr.treatment", "contr.poly"))
house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
prof.plr0 <- profileModel(house.plr, objective = function(fm) fm$deviance)
plot(prof.plr0)
# do it with a quantile
prof.plr1 <- update(prof.plr0, quantile = qchisq(0.95, 1))
plot(prof.plr1)
## End Example
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.