2
$\begingroup$

I'm looking at the implementation of the tweedie.profile function from the R package tweedie. I have a few questions. When trying to find the best power parameter $p$, by using the function tweedie.profile for example,

tweedie.profile(y ~ .,
                method = 'series',
                phi.method = 'mle'
                weights = w, 
                eps = 0,
                data = data)

Where w is a vector of weights used in the fitting the GLM. The code to obtain the MLE of $\phi$ is,

ans <- optimize(f = dtweedie.nlogl, maximum = FALSE, 
                interval = c(low.limit, 10 * phi.est), power = p, 
                mu = mu, y = ydata)

How come that the weights don't play any rule in the function? Shouldn't it be taken into account when computing $\hat{\phi}_{mle}$?

Here is a concrete example, it is well known that if the power parameter of tweedie $p = 0$ the distribution is Normal. Therefore, we can compare the appropriate GLM with the tweedie function.

 library(tweedie)
set.seed(1)
normal_y <- exp(rnorm(1000))
random_w <- runif(1000,0,10)

# No weights 
glm_model <- glm(normal_y ~ 1, family = gaussian("log"))
mle_full_tweedie <- tweedie.profile(
  normal_y ~ 1 ,
  p.vec = 0,
  offset = NULL,
  link.power = 0, # default is logarithmic link function
  do.smooth = FALSE,
  do.plot = FALSE,
  method = "series",
  verbose = 0,
  eps = 0
)


# With weights
glm_model_w <- glm(normal_y ~ 1, weights = random_w, family = gaussian("log"))
tweedie_model_w <- mle_full_tweedie <- tweedie.profile(
  normal_y ~ 1 ,
  weights = random_w,
  p.vec = 0,
  offset = NULL,
  link.power = 0, # default is logarithmic link function
  do.smooth = FALSE,
  do.plot = FALSE,
  method = "series",
  verbose = 0,
  eps = 0
)

Resulting with

> print(paste0('logLiklihood of Tweedie model (no weights):', mle_full_tweedie$L.max))
[1] "logLiklihood of Tweedie model (no weights):-2297.19437626455"
> print(paste0('logLiklihood of GLM model (no weights):', logLik(glm_model)))
[1] "logLiklihood of GLM model (no weights):-2297.09892399843"
> 
> print(paste0('logLiklihood of Tweedie model (weights):', tweedie_model_w$L.max))
[1] "logLiklihood of Tweedie model (weights):-2297.19437626455"
> print(paste0('logLiklihood of GLM model (weights):', logLik(glm_model_w)))
[1] "logLiklihood of GLM model (weights):-2439.90755930638".

Second, the function optimizes over the function dtweedie.nlogl which is a function that mixes the series and interpolation methods for estimating the likelihood. Why does it differ from the method set in the tweedie.profile function?

$\endgroup$
3
  • 1
    $\begingroup$ What weights are you referring to? $\endgroup$
    – jbowman
    Commented Jan 22, 2023 at 18:57
  • $\begingroup$ added an example $\endgroup$
    – Kozolovska
    Commented Jan 22, 2023 at 19:19
  • 2
    $\begingroup$ You might also want to add which package this is in; it's an easy internet search, but still... $\endgroup$
    – jbowman
    Commented Jan 22, 2023 at 19:24

1 Answer 1

3
$\begingroup$
  1. "How come that the weights don't play any rule in the function? Shouldn't it be taken into account when computing $\phi_{mle}$?".
  2. "... the function optimizes over the function dtweedie.nlogl which is a function that mixes the series and interpolation methods for estimating the likelihood. Why does it differ from the method set in the tweedie.profile function?"

To answer 1: Yes, the weights should come into play in the estimation of $\phi_{mle}$ (as you demonstrate). This is an error (and feel free to send through a fix :->).

To answer 2: The documentation specifies that phi.method is 'one of "saddlepoint" or "mle"'... the options do not include series or inversion as implied.

Estimating the mle of $\phi$ is quite tricky, and sometimes mle fails... (hence the 'fallback' option of using the saddlepoint (mean deviance) estimate if the mle fails).

In the case phi.method = "mle", the function tries its best to get an estimate... using any means it can find that works (series; inversion), and this cannot be user-specified (at present anyway).

Hope that helps.

P.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.