gllvm_posthoc_reprex.R

#gllvm reprex on posthocs
library(palmerpenguins)
penguins <- na.omit(penguins)

Y <- with(penguins, data.frame(bill_length_mm, bill_depth_mm, 
                          flipper_length_mm, body_mass_g))

X <- with(penguins, data.frame(species, sex))

mod_gllvm <- gllvm(y = Y, X = X,
                   formula = ~species*sex,
                   family = gaussian(link = "identity"),
                   num.lv=2,
                   starting.val = 'zero')

# Get the difference between sexes for each species
pred_frame <- expand.grid(species = levels(X$species),
                          sex = levels(X$sex))

#get predictions - but, no SE on predict
pred_frame <- cbind(pred_frame,
                    predict(mod_gllvm, newX = pred_frame,
                      level = 0) )

# no way around using dplyr and tidyr here!
library(dplyr)
library(tidyr)

diff_frame <- pred_frame %>%
    pivot_longer(bill_length_mm:body_mass_g) %>%
    pivot_wider(names_from = sex, values_from = value) %>%
    mutate(diff = female - male)

library(ggplot2)
ggplot(diff_frame, 
       aes(x = species, y = diff)) +
    facet_wrap(vars(name), scale = "free_x") +
    geom_point() +
    geom_hline(yintercept = 0, lty = 2) +
    coord_flip() +
    labs(x = "Difference between Female and Male")