0
$\begingroup$

I am attempting to write R code to generate bootstrap likelihood as described in section 3 of this paper https://arxiv.org/pdf/1510.07287.pdf.

I am confident that I performed the bootstraps correct, but I am new to KDE and I am confused on how to aggregate the KDE densities to bootstrap likelihoods. I have attached a completely reproducible example of my attempt at bootstrap likelihood.

The steps they use in section 3 of the paper are as follows.

  1. Generate K = 100 bootstraps of the data.

  2. Calculate a theta (parameter of interest) for each of the K bootstraps

  3. Generate L = 1000 second level bootstraps for each K boostraps

  4. Calculate theta for each of the L bootstraps

  5. Pass the theta calculated above to a kernel function to estimate the distribution

  6. Log the KDE distribution to get the log likelihood of the K thetas

Using the mtcars data and the model $mpg = \theta*cyl$. I wanted to build a bootstrap likelihood for $\theta$. The point estimate from lm(mpg ~ cyl, data = mtcars)is about -2.9.

My likelihood function seems reasonable and is shown below:

enter image description here

# install.packages("strapgod")
# install.packages("tidyverse")
# install.packages("KernSmooth")
library(strapgod)
library(tidyverse)
library(KernSmooth)
data(mtcars)

# check point estimate for cyl coeff
f <- lm(mpg ~ cyl, data = mtcars)
summary(f)

###########################
# Set up vars and functions
###########################
K = 100  # firest level boot
L = 1000 # second level boot

# model and pull out cyl coeff
lm_fun <- function(data){
  fit <- lm(mpg ~ cyl, data = data)
  summary(fit)$coefficients[2] 
}

# kernel density functions
kern_funy <- function(params){
  p <- params$model
  est <- bkde(p, bandwidth=0.25)
  est$y
}

kern_funx <- function(params){
  p <- params$model
  est <- bkde(p, bandwidth=0.25)
  est$x
}

# Execute the method
# ------------------
set.seed(42)

# K bootstrap samples
P <- bootstrapify(mtcars, K) %>% 
  collect() 

# generate theta_i estimates
thetai <- P %>% 
  group_by(.bootstrap) %>% 
  nest() %>% 
  mutate(theta = map_dbl(data, lm_fun)) 

# generate theta_ij estimates
thetaij <- P %>% 
  group_by(.bootstrap) %>% 
  bootstrapify(L) %>% 
  collect() %>% 
  group_by(.bootstrap, .bootstrap1) %>% 
  nest() %>% 
  mutate(model = map_dbl(data, lm_fun)) %>% 
  unnest(model)

# get kernel density from thetaij
kd <- thetaij %>% 
  select(.bootstrap, model) %>% 
  group_by(.bootstrap) %>% 
  nest() %>% 
  mutate(kernx = map(data, kern_funx)) %>% 
  mutate(kerny = map(data, kern_funy)) %>% 
  unnest(kernx, kerny) %>% 
  group_by(kernx) %>% 
  summarise(avg_p = mean(kerny))

# look up fun to find likelihood for a theta
get_p <- function(x){
  kd %>% 
    filter(abs(kernx - x) == min(abs(kernx - x))) %>% 
    pull(avg_p)
}

# join params to kern
ll <- thetai %>% 
  rowwise() %>% 
  mutate(kerny = get_p(theta)) %>% 
  mutate(ll = log(kerny))

# Plot
ll %>% 
  ggplot(aes(x = theta, y = ll)) +
  geom_point(shape = "o") +
  stat_smooth(method = "lm", formula = y ~ x + I(x^2), size = 1, se = F) +
  ggtitle("log bootstrap likelihood") +
  theme_classic()

$\endgroup$

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.