# Analysis

    library(tidyverse)
    library(tidybayes)
    library(targets)
    library(marginaleffects)
    library(modelsummary)
    library(parameters)
    library(patchwork)

    tar_load(penguins)
    tar_load(c(model1_bayes, model2_bayes))

# Look at data

    ggplot(penguins, aes(x = bill_len, y = body_mass, color = species)) +
      geom_point()

# Look at models

Here is model 1

    model_parameters(model1_bayes, verbose = FALSE)

    # Fixed Effects

    Parameter   | Median |            95% CI |     pd |  Rhat |  ESS
    ----------------------------------------------------------------
    (Intercept) | 394.74 | [-203.64, 962.24] | 90.33% | 1.000 | 3888
    bill_len    |  86.72 | [  73.86, 100.09] |   100% | 1.000 | 3907

    # sigma Parameters

    Parameter | Median |           95% CI |   pd |  Rhat |  ESS
    -----------------------------------------------------------
    sigma     | 653.05 | [604.48, 708.34] | 100% | 1.000 | 3630

Table 1 shows both at the same time

    modelsummary(
      list(model1_bayes, model2_bayes),
      statistic = "[{conf.low}, {conf.high}]",
      ci_method = "hdi",
      metrics = c("R2"),
      fmt = 1
    )

Table 1: Both models at the same time

| (1) | (2)
--- | --- | ---
b_Intercept | 394.7 | 201.7
| [-209.7, 948.8] | [-332.1, 775.9]
b_bill_len | 86.7 | 90.3
| [74.1, 100.2] | [76.0, 104.4]
sigma | 653.1 | 375.5
| [601.1, 703.9] | [347.0, 404.5]
b_speciesChinstrap |  | -876.2
|  | [-1054.5, -700.1]
b_speciesGentoo |  | 599.3
|  | [443.2, 753.3]
Num.Obs. | 333 | 333
R2 | 0.347 | 0.784

Check out Figure 1 for predicted values

    p1 <- model1_bayes |>
      plot_predictions(condition = c("bill_len")) +
      coord_cartesian(ylim = c(2000, 6500))

    p2 <- model2_bayes |>
      plot_predictions(condition = c("bill_len", "species")) +
      coord_cartesian(ylim = c(2000, 6500))

    p1 | p2

Figure 1: Predicted weight across bill length

And look at the `bill_len` coefficient in both of the models in Figure 2.

    bill_len_model_1 <- model1_bayes |>
      gather_draws(b_bill_len) |>
      mutate(model = "Model 1")

    bill_len_model_2 <- model2_bayes |>
      gather_draws(b_bill_len) |>
      mutate(model = "Model 2")

    bind_rows(bill_len_model_1, bill_len_model_2) |>
      ggplot(aes(x = .value, y = model, fill = model)) +
      ggdist::stat_halfeye() +
      guides(fill = "none") +
      labs(x = "bill_len")

Figure 2: Coefficient for bill length across the two models