In mneunhoe/simloglm: Simulate Correct Quantities of Interest from Linear Models with Logged Dependent Variables

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(simloglm)

ff <- "y~x1*x2"
n <- 5000

set.seed(20220309)
pop_df <-
  simloglm:::setup_pop(
    n,
    ff,
    zero_centered = F,
    type = c("num", "bin"),
    coefs = c(5, 0.5,-5.5, 0.4),
    noise_sd = 1.5
  )

reg_ia <- lm(ff, data = pop_df)
x1_seq <- seq(min(pop_df$x1), max(pop_df$x1), length.out = 100)

# Set seeds to make sure simulations are based on same draws of coefficients
set.seed(20220310)
scenario_1 <- simloglm(reg_ia, scenario = list(x1 = x1_seq, x2 = 1))
set.seed(20220310)
scenario_0 <- simloglm(reg_ia, scenario = list(x1 = x1_seq, x2 = 0))

Make plots

m <-matrix(c(1, 2, 3, 4), ncol = 2, byrow = T)
layout(m,
       widths = c(1,1))


col_vec <- viridis::viridis(2)
# Plot 1: Regressions on log scale
par(mar = c(5,     # bottom (5)
            6.5,   # left   (4)
            2,     # top    (4)
            0.5    # right  (2)
),
xpd = T)




    plot(
      pop_df$x1,
      pop_df$y,
      type = "p",
      pch = 16,
      col = adjustcolor(col_vec[pop_df$x2+1], 0.5),
      main = "",
      xlab = "",
      ylab = "",
      yaxt = "n",
      xaxt = "n",
      bty = "n",
      ylim = c(-1, 16)
    )
    box(lty = "solid", col = "grey50")


    axis(
      2,
      at = pretty(pop_df$y)[2:(length(pretty(pop_df$y))-1)],
      las = 2,
      col.axis = "grey50",
      col = "grey50"
    )


    axis(1,
         at = pretty(pop_df$x1)[2:(length(pretty(pop_df$x1))-1)],
         col.axis = "grey50",
         col = "grey50")
    mtext(2,
          text = "log(Y)",
          line = 4.1)
    mtext(3,
          text = expression(bold("Log-Scale")),
          line = 0.4)
    mtext(1,
          text = "X",
          line = 3.5)

    text(9, 3, labels = "D = 1", col = col_vec[2], font = 2)
    text(4, 12, labels = "D = 0", col = col_vec[1], font = 2)

    # Median on log scale

    #quantiles_log_1 <- apply(log(scenario_1$geometric_mean),2, quantile, c(0.025, 0.975))
    scenario_1_log_scale <- lapply(get_summary(scenario_1), log)
    quantiles_log_1 <- scenario_1_log_scale$quantiles
    polygon(
      c(rev(x1_seq), x1_seq),
      c(rev(quantiles_log_1[1, ]), quantiles_log_1[2, ]),
      col = adjustcolor(col_vec[2], 0.6),
      border = NA
    )

    lines(x1_seq, scenario_1_log_scale$point_estimate)

    scenario_0_log_scale <- lapply(get_summary(scenario_0), log)
    quantiles_log_0 <- scenario_0_log_scale$quantiles

    polygon(
      c(rev(x1_seq), x1_seq),
      c(rev(quantiles_log_0[1, ]), quantiles_log_0[2, ]),
      col = adjustcolor(col_vec[1], 0.6),
      border = NA
    )

    lines(x1_seq, scenario_0_log_scale$point_estimate)

    par(mar = c(5,     # bottom (5)
            6.5,   # left   (4)
            2,     # top    (4)
            0.5    # right  (2)
),
xpd = T)


    plot(
      pop_df$x1,
      exp(pop_df$y),
      type = "p",
      pch = 16,
      col = adjustcolor(col_vec[pop_df$x2+1], 0.5),
      main = "",
      xlab = "",
      ylab = "",
      yaxt = "n",
      xaxt = "n",
      bty = "n"
    )
    box(lty = "solid", col = "grey50")


    axis(
      2,
      at = pretty(exp(pop_df$y))[2:(length(pretty(pop_df$y))-1)],
      las = 2,
      col.axis = "grey50",
      col = "grey50"
    )



    axis(1,
         at = pretty(pop_df$x1)[2:(length(pretty(pop_df$x1))-1)],
         col.axis = "grey50",
         col = "grey50")
    mtext(2,
          text = "Y",
          line = 4.1)
    mtext(3,
          text = expression(bold("Original-Scale")),
          line = 0.4)
    mtext(1,
          text = "X",
          line = 3.5)


    # Conditional Mean

    scenario_1_original_scale <- get_summary(scenario_1, which_qoi = "mean")
    quantiles_original_1 <- scenario_1_original_scale$quantiles

    polygon(
      c(rev(x1_seq), x1_seq),
      c(rev(quantiles_original_1[1, ]), quantiles_original_1[2, ]),
      col = adjustcolor(col_vec[2], 0.6),
      border = NA
    )

    lines(x1_seq, scenario_1_original_scale$point_estimate)

    scenario_0_original_scale <- get_summary(scenario_0, which_qoi = "mean")
    quantiles_original_0 <- scenario_0_original_scale$quantiles

    polygon(
      c(rev(x1_seq), x1_seq),
      c(rev(quantiles_original_0[1, ]), quantiles_original_0[2, ]),
      col = adjustcolor(col_vec[1], 0.6),
      border = NA
    )

    lines(x1_seq, scenario_0_original_scale$point_estimate)


    # Plot 3: FD on log scale
par(mar = c(5,     # bottom (5)
            6.5,   # left   (4)
            2,     # top    (4)
            0.5    # right  (2)
),
xpd = T)


#plot(x1_seq, apply((scenario_0$arithmetic_mean),2, function(x) exp(mean(log(x)))) - apply((scenario_1$arithmetic_mean),2, function(x) exp(mean(log(x)))), type = "l")

fd_quantiles <- apply(log(scenario_0$median) - log(scenario_1$median), 2, quantile, c(0.025, 0.975))

fd <- log(scenario_0$median_point_estimate) - log(scenario_1$median_point_estimate)

    plot(
    c(rev(x1_seq), x1_seq),
      c(rev(fd_quantiles[1, ]), fd_quantiles[2, ]),
      type = "n",
      pch = 16,
      col = adjustcolor(col_vec[pop_df$x2+1], 0.5),
      main = "",
      xlab = "",
      ylab = "",
      yaxt = "n",
      xaxt = "n",
      bty = "n"
    )
    box(lty = "solid", col = "grey50")


    axis(
      2,
      at = pretty(fd_quantiles)[2:(length(pretty(fd_quantiles))-1)],
      las = 2,
      col.axis = "grey50",
      col = "grey50"
    )


    axis(1,
         at = pretty(x1_seq)[2:(length(pretty(x1_seq))-1)],
         col.axis = "grey50",
         col = "grey50")
    mtext(2,
          text = "First Difference on Log-Scale",
          line = 4.1)
    mtext(3,
          text = expression(bold("Log-Scale")),
          line = 0.4)
    mtext(1,
          text = "X",
          line = 3.5)


    polygon(
      c(rev(x1_seq), x1_seq),
      c(rev(fd_quantiles[1, ]), fd_quantiles[2, ]),
      col = adjustcolor(col_vec[1], 0.6),
      border = NA
    )

    lines(x1_seq, fd)

       # Plot 4: FD on original scale
par(mar = c(5,     # bottom (5)
            6.5,   # left   (4)
            2,     # top    (4)
            0.5    # right  (2)
),
xpd = T)


#plot(x1_seq, apply((scenario_0$arithmetic_mean),2, function(x) exp(mean(log(x)))) - apply((scenario_1$arithmetic_mean),2, function(x) exp(mean(log(x)))), type = "l")

fd_quantiles <- apply((scenario_0$mean) - (scenario_1$mean), 2, quantile, c(0.025, 0.975))

fd <- (scenario_0$mean_point_estimate) - (scenario_1$mean_point_estimate)
    plot(
    c(rev(x1_seq), x1_seq),
      c(rev(fd_quantiles[1, ]), fd_quantiles[2, ]),
      type = "n",
      pch = 16,
      col = adjustcolor(col_vec[pop_df$x2+1], 0.5),
      main = "",
      xlab = "",
      ylab = "",
      yaxt = "n",
      xaxt = "n",
      bty = "n"
    )
    box(lty = "solid", col = "grey50")


    axis(
      2,
      at = pretty(fd_quantiles)[2:(length(pretty(fd_quantiles))-1)],
      las = 2,
      col.axis = "grey50",
      col = "grey50"
    )


    axis(1,
         at = pretty(x1_seq)[2:(length(pretty(x1_seq))-1)],
         col.axis = "grey50",
         col = "grey50")
    mtext(2,
          text = "First Difference on Original-Scale",
          line = 4.1)
    mtext(3,
          text = expression(bold("Original-Scale")),
          line = 0.4)
    mtext(1,
          text = "X",
          line = 3.5)


    polygon(
      c(rev(x1_seq), x1_seq),
      c(rev(fd_quantiles[1, ]), fd_quantiles[2, ]),
      col = adjustcolor(col_vec[1], 0.6),
      border = NA
    )

    lines(x1_seq, fd)

mneunhoe/simloglm documentation built on June 14, 2022, 5:18 p.m.