Skip to content

Instantly share code, notes, and snippets.

@eric-czech
Last active February 22, 2021 18:46
Show Gist options
  • Save eric-czech/ec803e34d85377950ee9b2eba3d3d7da to your computer and use it in GitHub Desktop.
Save eric-czech/ec803e34d85377950ee9b2eba3d3d7da to your computer and use it in GitHub Desktop.
Code for figure 2 from Bellégo and Papeis 2019
library(ggplot2)
# For figure 2 of paper https://www.semanticscholar.org/paper/Dealing-with-Logs-and-Zeros-in-Regression-Models-Bellégo-Pape/546a6f45433b413721f9a60f0be8e3e2b69fe103
set.seed(1)
beta <- 1
x_max <- 1
x <- x_max * runif(10000)
y <- sapply(x, function(x) rpois(1, lambda=exp(beta * x + 1)))
df <- do.call(rbind, lapply(seq(.01, 1, by=.01), function(delta){
beta_hat <- unname(lm(log(y + delta) ~ x)$coefficients['x'])
bias <- 100 * abs(beta_hat - beta) / beta
data.frame(bias=bias, delta=delta)
}))
ggplot(df, aes(x=delta, y=bias)) + geom_line()
# Adaptation w/ realistic parameterization + log base
set.seed(1)
beta <- 1
x_max <- 7
x <- x_max * runif(10000)
y <- sapply(x, function(x) rpois(1, lambda=2**(beta * x + 1)))
df <- do.call(rbind, lapply(seq(.01, 1, by=.01), function(delta){
beta_hat <- unname(lm(log2(y + delta) ~ x)$coefficients['x'])
bias <- 100 * abs(beta_hat - beta) / beta
data.frame(bias=bias, delta=delta)
}))
ggplot(df, aes(x=delta, y=bias)) + geom_line()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment