tweedie.profile(formula, p.vec=NULL, xi.vec=NULL, link.power=0,
data, weights, offset, fit.glm=FALSE,
do.smooth=TRUE, do.plot=FALSE, do.ci=do.smooth,
eps=1/6,
control=list( epsilon=1e-09, maxit=glm.control()$maxit, trace=glm.control()$trace ),
do.points=do.plot, method="inversion", conf.level=0.95,
phi.method=ifelse(method == "saddlepoint", "saddlepoint", "mle"),
verbose=FALSE, add0=FALSE)response ~ predictors.
For details,
see the documentation for lm,
p values for consideration.
The values must all be larger than one
(if the response variable has exact zeros,
the values must all be between one and two).
If NULL (the default),
p.vec is sep.vec;
some authors use the $p$ notation for the index parameter,
and some use $\xi$;
this function detects which is used and then uses that notation throughoutlink.power=0 (the default)
refers to the logarithm link function.
See the documentation fas.data.frame to a data frame)
containing the variables in the model.
If not found in data,
the variables are taken from environment(formuNULL or a numeric vector.NULL or a numeric vector of length either one or
equal to the number of cases.
One or moTRUE,
the Tweedie generalized linear model is fitted using the value of $p$
found by the profiling function.
If FALSE (the default),
no model is fitted.TRUE (the default),
a spline is fitted to the data to smooth the profile likelihood plot.
If FALSE,
no smoothing is used
(and the function is quicker).
Note that p.vec must TRUE,
a plot of the profile likelihood is produce.
If FALSE (the default),
no plot is produced.TRUE,
the nominal 100*conf.level
is computed.
If FALSE,
no confidence interval is computed.
By default,
do.ci is the same value as do.smooth,
since aeps=1/6
(as suggested by Nelder and Pregibon, 1987).
Note eps is ignored unless the
method="saddlepoint"
as it makes no sense otherwise.glm.control and glm.
The default is to use the maximum number of iterations maxitp;
defaults to the same value as do.plot"series",
"inversion" (the default),
"interpolation"
or
"saddlepoint".
If there are any troubles using this function,
often conf.level=0.95.phi,
one of
"saddlepoint"
or
"mle".
A maximum likelihood estimate is used unless
method="saddlepoint",
when the saddlepoint approximation method is used.
0 or FALSE means minimal feedback (the default),
1 or TRUE means some feedback,
or 2 means to show all feedback.
Since the function can be slowTRUE, the value p=0 is used in forming the profile lohg-likelihood
(corresponding to the normal distribution);
the default value is add0=FALSEp.max.
Optionally (if do.plot=TRUE),
a plot is produced that shows the profile log-likelihood
computed at each value in p.vec
(smoothed if do.smooth=TRUE).
This function can be tempermental
(for theoretical reasons involved in numerically computing the density),
and this plot shows the values of $p$ requested on the
horizontal axis (using rug);
there may be fewer points on the plot,
since the likelihood some values of $p$ requested
may have returned NaN, Inf or NA.
A list containing the components:
y and x
(such that plot(x,y) (partially)
recreates the profile likelihood plot);
ht (the height of the nominal confidence interval);
L (the estimate of the (log-) likelihood at each given value of p);
p (the p-values used);
phi (the computed values of phi at the values in p);
p.max (the estimate of the mle of p);
L.max (the estimate of the (log-) likelihood at p.max);
phi.max (the estimate of phi at p.max);
ci (the lower and upper limits of the confidence interval for p);
method (the method used for estimation: series, inversion,
interpolation or saddlepoint);
phi.method (the method used for estimation of phi:
saddlepoint or phi).
If glm.fit is TRUE,
the list also contains a component glm.obj,
a glm object for the fitted Tweedie generalized linear model.p.vec,
the function computes an estimate of phi
and then computes the value of the log-likelihood for these parameters.
The plot of the log-likelihood against p.vec
allows the maximum likelihood value of p
to be found.
Once the value of p is found,
the distribution within the class of Tweedie distribution is identified.dtweedie,
dtweedie.saddle,
tweedielibrary(statmod) # Needed to use tweedie.profile
# Generate some fictitious data
test.data <- rgamma(n=200, scale=1, shape=1)
# The gamma is a Tweedie distribution with power=2;
# let's see if p=2 is suggested by tweedie.profile:
out <- tweedie.profile( test.data ~ 1,
p.vec=seq(1.5, 2.5, by=0.2) )
out$p.max
out$ciRun the code above in your browser using DataCamp Workspace