psu betaregression brms.R

library(pacman)
p_load(tidyverse, magrittr, pxweb, janitor, glmmTMB, brms, emmeans)

options(scipen=999)

psu <- get_pxweb_data(url = "http://api.scb.se/OV0104/v1/doris/sv/ssd/ME/ME0201/ME0201B/Partisympati17",
               dims = list(Kon = c('1', '2'),
                           UtbNivaSUN2000 = c('F', '3', 'K', 'L'),
                           Partisympati = c('m', 'c', 'l', 'kd', 'mp', 's', 'v', 'SD', 'övr'),
                           ContentsCode = c('*'),
                           Tid = c('2018M11')),
               clean = TRUE) %>% 
  clean_names() %>% as_tibble %>% spread(contents_code, values) %>% 
  clean_names %>% 
  mutate(utbildningsniva_sun_2000 = utbildningsniva_sun_2000 %>% 
           ordered(levels=c("förgymnasial utbildning", 
                            "gymnasial utbildning", 
                            "eftergymnasial utbildning mindre än 3 år", 
                            "eftergymnasial utbildning 3 år eller mer"))) %>% 
  mutate(mu = svarsfordelning_procent/100,
               variance = ((felmarginal/100)/1.96)^2,
               phi = (mu*(1-mu))/variance-1)

model_sd <- brm(bf(mu~kon*utbildningsniva_sun_2000, 
                   phi ~ offset(phi)-1), 
                data=psu %>% filter(partisympati=="SD"), 
                family=Beta(link_phi="identity"))

model_sd %>% 
  emmeans(~kon*utbildningsniva_sun_2000, transform="response") %>% 
  contrast(method="pairwise", by="utbildningsniva_sun_2000") %>% 
  plot+geom_vline(xintercept = 0)

model_c <- update(model_sd, newdata=psu %>% filter(partisympati=="C"))
model_kd <- update(model_sd, newdata=psu %>% filter(partisympati=="KD"))
model_l <- update(model_sd, newdata=psu %>% filter(partisympati=="L"))
model_m <- update(model_sd, newdata=psu %>% filter(partisympati=="M"))
model_mp <- update(model_sd, newdata=psu %>% filter(partisympati=="MP"))
model_s <- update(model_sd, newdata=psu %>% filter(partisympati=="S"))
model_v <- update(model_sd, newdata=psu %>% filter(partisympati=="V"))
model_ovriga <- update(model_sd, newdata=psu %>% filter(partisympati=="övriga"))

model_l %>% marginal_effects(effects="kon:utbildningsniva_sun_2000", robust=FALSE, method="predict") %>% 
  plot %>% .[[1]]+scale_y_continuous(labels=scales::percent, limits=c(0,NA))