In yasui-salmon/CQR:

knitr::opts_chunk$set(echo = TRUE)

read libraries

library(dplyr)
library(tidyr)
library(CQR)
library(quantreg)

read dataset

data <- read.csv(file = "C:/Users/A13375/Documents/GitHub/CQR/EXAMPLE/data/extramartial_affair.csv")

quantile regression

qr_result <- rq(data = data,
                tau = seq(0.5, 0.95, 0.05),
                method = "fn",
                formula = affairs ~ rate_marriage + age + 
                  yrs_married + children + religious + 
                  educ + occupation + occupation_husb)

plot(qr_result)

1st step

taus <- seq(0.5, 0.95, 0.05)
q1 <- 0.1
cut_value <- 0
q2 <- 0.09

YV <- "affairs"
XV <- c("rate_marriage", "age",
        "yrs_married", "children", "religious",
        "educ, occupation", "occupation_husb")

cqr_data <- data %>% select(affairs, rate_marriage, age,
                            yrs_married, children, religious, 
                            educ, occupation, occupation_husb)



sq_extra <- cqr_data %>% select_(paste("-",YV, sep = "")) %>% mutate_at(vars(-matches(YV)), pt, n = 2)
colnames(sq_extra) <- paste("sq_", colnames(sq_extra), sep = "")

cub_extra <- cqr_data %>% select_(paste("-",YV, sep = "")) %>% mutate_at(vars(-matches(YV)), pt, n = 3)
colnames(cub_extra) <- paste("cub_", colnames(cub_extra), sep = "")

inter <- inter_party(cqr_data %>% select_(paste("-",YV, sep = "")))
lr_extra <- cbind(cqr_data, sq_extra, cub_extra, inter) %>%
  mutate(binom_var = ifelse(affairs > 0, 1, 0))

first_logit <- glm(lr_extra, formula = binom_var ~ . - affairs, family = binomial())
step_logit <- step(first_logit, trace = 0)

• You need to build your own step 1, since this part is needed to be modify with your data.
• Basically you need to predict the probability of not censoring with your treatment and control variables.
• Logit model might be enough here but some other model is also your available option.
• in this example, X, X^2, X^3 and interaction term within X are used in the logit model then used step AIC variable selection to discard.(This process is not in the original paper, also there was not enough information about this process in there.)
• the output of 1st step must be the model object.

2nd and 3rd step

three_step_result <- two_three_step(first_step_model = step_logit, 
                                    taus = taus, 
                                    q1 = q1, 
                                    q2 = q2, 
                                    YV = YV, 
                                    XV = XV, 
                                    cqr_data = cqr_data)

• this part is the 2nd and 3rd step estimation.
• taus could be something between 0 and 1.
• q1 is c in the Step 1 of the original paper.
• 0.1 works well accoring to their simulation result.
• q2 is $\delta$ in the Step2 of the original paper.
• original paper did not explain about the optimal value in here.
• but stata implementation called CQIV hired 0.03 for this value, so it may work.
• this part produce a list object
• [[1]] contains the list of three-step estimator by tau
  • if you run a model with large dataset, this might be trouble some since model object contin dataset itself.
• [[2]] contains the robustness statistics which recommended to check.
  • J1_rate_total is N(J1)/N(total)
  • J1_rate_J0 is number of sample in J1 also included in J0 divided in N(J1)
  • J1_notin_J0 is number of sample in J1 not included in J0
• One of the advantages in this method is the standard error and confidence intervals produced in the last step Quantile Regression are all valid without any modification.

result plot

sum_coef <- lapply(three_step_result[[1]],"[[", 3)
taus <- sapply(three_step_result[[1]],"[[", 6)
coefficients <- lapply(sum_coef, function(x){x[,1] %>% t() %>% as.data.frame()}) %>% bind_rows
stder <- lapply(sum_coef, function(x){x[,2] %>% t() %>% as.data.frame()}) %>% bind_rows

coef_df <- data.frame(coefficients, tau = taus) %>% gather(vid, value, -tau)
stder_df <- data.frame(stder, tau = taus) %>% gather(vid, stder, -tau)
plot_df <- coef_df %>% inner_join(stder_df, by = c("vid", "tau"))

plot_df %>%
  ggplot(aes(y = value, x = tau)) +
  facet_wrap(~vid, scales = "free") +
  geom_point() + 
  geom_line() +
  geom_errorbar(plot_df, 
                mapping = aes(x = tau, ymin = value - (stder*2), ymax = value + (stder*2)),
                width = 0.0005, size = 1, color="blue")

• The Above plot is the reproduction of Figure 4. in the original paper.
• Since there was no explanation about delta in the section 3, delta might be differ from the original environment and which may cause the small gap in this result from the original result.
• Also, original paper reported the size of dataset as N = 6388, but the dataset used in here is N = 6366, so that this difference also might the reason of the gap in these results.
• Though, the tendency in the changing of the parameters and their values are similar with the original.

yasui-salmon/CQR documentation built on May 20, 2019, 12:32 p.m.