data |> 
  select(INCTOT,EMPSTAT) |> 
  mutate(incstat=as.integer(INCTOT),
         empstat=as.factor(EMPSTAT),
         incstat=na_if(incstat,9999999),
         empstat=case_when(empstat==0 ~ "NA",
                           empstat==1 ~ "Employed",
                           empstat==2 ~ "Unemployed",
                           empstat==3 ~ "Not in Labour Force"),
         incstat=incstat/100000) |> 
  ggplot(aes(x=incstat,y=..scaled..,fill=empstat,color=empstat))+
  geom_density(alpha=0.3,na.rm=TRUE)+
  xlim(0,8)+
  scale_y_continuous(labels = scales::label_number())+
  scale_fill_discrete(name="Employment Status")+
  scale_color_discrete(name="Employment Status")+
  theme_economist(dkpanel=TRUE)+
  labs(title="2019 Income Distribution by Employment Status",
       x="Income: Values Scaled by 1/100000",
       caption="Source: IPUMS,2019")