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

## ----setup--------------------------------------------------------------------
library(csranks)
data(parent_child_income)
head(parent_child_income)

## ----lmranks------------------------------------------------------------------
lmr_model <- lmranks(r(c_faminc) ~ r(p_faminc), data = parent_child_income)
summary(lmr_model)

## ----lm-----------------------------------------------------------------------
c_faminc_rank <- frank(parent_child_income$c_faminc, omega = 1, increasing = TRUE)
p_faminc_rank <- frank(parent_child_income$p_faminc, omega = 1, increasing = TRUE)
lm_model <- lm(c_faminc_rank ~ p_faminc_rank)
summary(lm_model)

## ----lmrankscov---------------------------------------------------------------
lmr_model_cov <- lmranks(r(c_faminc) ~ r(p_faminc) + gender + race, data = parent_child_income)
summary(lmr_model_cov)

## ----grouped_lmranks_simple---------------------------------------------------
grouped_lmr_model_simple <- lmranks(r(c_faminc) ~ r(p_faminc_rank):gender,
  data = parent_child_income
)
summary(grouped_lmr_model_simple)

## ----grouped_lm_simple--------------------------------------------------------
grouped_lm_model_simple <- lm(c_faminc_rank ~ p_faminc_rank:gender + gender - 1, # group-wise intercept
  data = parent_child_income
)
summary(grouped_lm_model_simple)

## ----grouped_lmranksgran------------------------------------------------------
parent_child_income$subgroup <- interaction(parent_child_income$gender, parent_child_income$race)
grouped_lmr_model <- lmranks(r(c_faminc) ~ r(p_faminc_rank):subgroup,
  data = parent_child_income
)
summary(grouped_lmr_model)

## ----grouped_lm---------------------------------------------------------------
grouped_lm_model <- lm(c_faminc_rank ~ p_faminc_rank:subgroup + subgroup - 1, # group-wise intercept
  data = parent_child_income
)
summary(grouped_lm_model)

## ----plot_CIs, message = FALSE, out.width = "90%", fig.width=6, fig.height=4----
library(ggplot2)
theme_set(theme_minimal())
ci_data <- data.frame(
  estimate = coef(lmr_model),
  parameter = c("Intercept", "slope"),
  group = "Whole sample",
  method = "csranks",
  lower = confint(lmr_model)[, 1],
  upper = confint(lmr_model)[, 2]
)

ci_data <- rbind(ci_data, data.frame(
  estimate = coef(grouped_lmr_model),
  parameter = rep(c("Intercept", "slope"), each = 6),
  group = rep(c(
    "Hispanic female", "Hispanic male", "Black female", "Black male",
    "Other female", "Other male"
  ), times = 2),
  method = "csranks",
  lower = confint(grouped_lmr_model)[, 1],
  upper = confint(grouped_lmr_model)[, 2]
))

ci_data <- rbind(ci_data, data.frame(
  estimate = coef(lm_model),
  parameter = c("Intercept", "slope"),
  group = "Whole sample",
  method = "naive",
  lower = confint(lm_model)[, 1],
  upper = confint(lm_model)[, 2]
))

ci_data <- rbind(ci_data, data.frame(
  estimate = coef(grouped_lm_model),
  parameter = rep(c("Intercept", "slope"), each = 6),
  group = rep(c(
    "Hispanic female", "Hispanic male", "Black female", "Black male",
    "Other female", "Other male"
  ), times = 2),
  method = "naive",
  lower = confint(grouped_lm_model)[, 1],
  upper = confint(grouped_lm_model)[, 2]
))

ggplot(ci_data, aes(y = estimate, x = group, ymin = lower, ymax = upper, col = method, fill = method)) +
  geom_point(position = position_dodge2(width = 0.9)) +
  geom_errorbar(position = position_dodge2(width = 0.9)) +
  geom_hline(aes(yintercept = estimate),
    data = subset(ci_data, group == "Whole sample"),
    linetype = "dashed",
    col = "gray"
  ) +
  coord_flip() +
  labs(title = "95% confidence intervals of intercept and slope\nin rank-rank regression") +
  facet_wrap(~parameter)