I have been attempting to replicate an experiment from Computer Age Statistical Inference by Bradley Efron and Trevor Hastie on page 411.

In this experiment 100 datasets are populated normal random values. These values are then then corrected by subtracting from Tweedie's estimates. In the book the following plot is produced.

enter image description here

I made an attempt to replicate this experiment, but I believe my Tweedie estimates are wrong.

I was hoping someone could point out my error.



make_z <- function(){
  # mu is e^{-mu} distributed where (mu > 0)
  u <- -rexp(n = 1000)
  z <- rnorm(u)
  tibble(z, u)

get_tweedie <- function(data){
  z <- data$z
  u <- data$u
  # Get Tweedie's estimate using a natrual spline model with 5 degrees of freedom
  fit <- glm(z ~ ns(u, 5), data = data, family=tweedie(var.power=0, link.power=1L))
  exp(predict(fit, type = "response"))

df_z <- tibble(dataset = seq(1:100)) %>% 
  rowwise() %>% 
  mutate(Z = list(make_z())) %>% 
  ungroup() %>% 
  mutate(tweedie = map(Z, get_tweedie)) %>% 
  unnest(c(Z, tweedie)) %>% 
  arrange(desc(u)) %>% 
  group_by(dataset) %>% 
  # For each dataset, the 20 largest z values and the corresponding tweedie
  # estimates and actual mu values were recorded
  top_n(20, z) %>% 
  # uncorrected and corrected differences were calculated. The hope being
  # that empirical Bayes shrinkage would correct the selection bais in the
  # z values.
  mutate(uncorr = z - u,
         corrected = tweedie - u)

# ggplot
df_z %>% 
  ungroup() %>% 
  select(uncorr, corrected) %>% 
  gather() %>% 
  ggplot(aes(x = value, fill = key)) +
  geom_histogram(color = "white", bins = 45) +

enter image description here


I am putting the experiment details here for convenience.

enter image description here

  • 2
    $\begingroup$ I am the author of the statmod::tweedie function used in your code. I think you might be confusing different statistical methods attributed to Tweedie. Tweedie's formula used in Efron & Hastie (2016) actually has nothing to do with Tweedie distributions or with Tweedie glms. Given the settings for var.power and link.power, your code line using tweedie is actually the same as fit <- glm(z ~ ns(u, 5), data = df_1, family=gaussian()). $\endgroup$ – Gordon Smyth Nov 24 '19 at 6:02
  • 1
    $\begingroup$ ... so you're actually just fitting an ordinary linear model by least squares. You may as well have used lm. You are not doing any bias correction as far as I can see. $\endgroup$ – Gordon Smyth Nov 24 '19 at 6:12
  • $\begingroup$ Thanks! I did not realize the distinction. Would you have any insight on how I can calculate Tweedie's formula in this situation? All descriptions I have found in the book are analytical formulas for $\hat{\mu_i} = z_i + \hat{l}'(z_i)$. I have not been able to find any code that estimates this value. $\endgroup$ – Alex Nov 24 '19 at 12:48
  • $\begingroup$ That's because there isn't any code. Efron & Hastie's book is about research methods that are not implemented in standard public code packages. BTW, another problem with your code is that you have not generated the u from an exponential distribution. Any reason why you wouldn't use rexp? $\endgroup$ – Gordon Smyth Nov 24 '19 at 22:23
  • $\begingroup$ I should have used rexp thanks for pointing that out. $\endgroup$ – Alex Nov 25 '19 at 1:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.