Description Usage Arguments Details Value Note Author(s) References See Also Examples
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
prespecified 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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19  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 = 1e08, 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 = 1e08, ...)

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. 
fitted
has to be an object which supports the method
coef
and which has fitted$terms
with the same
meaning as, for example, in lm
and
glm
(see also terms
). coef(fitted)
has to be a vector of coefficients with each component
corresponding to a column of the model matrix returned by
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.
The fitting procedure that resulted fitted
has to support
offset
in formula
.
Also, fitted$call
has to be the call that generated fitted
.
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 userspecified 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_{i1},
D(F(b_{i1}))]. + 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 overconservative 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_{M1} with 2C instead of C in the step
calculation. If the three last slopes were less than 1e8
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 1e8 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
.
profileModel
has been tested and is known to work for fitted
objects resulting from lm
, glm
,
polr
, gee
, geeglm
, brglm
and BTm
.
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.
confintModel
, plot.profileModel
.
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  ## 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)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.