Clarification about counter value in generate_pseudo_pop()

[m-qin]

Hello! Falco and I are working on a project and we are trying to understand the value of generate_pseudo_pop(), specifically the item counter. Is counter, a vector of length equal to the length of the pseudopopulation, representing how many times each real observation (that is untrimmed by the matching algorithm) is matched in the pseudopopulation?

Thanks!

[m-qin]

Also, if I run sum(counter), should I get length(counter) * {number of GPS bins}? And is {number of GPS bins} = (max(gps) - min(gps)) / delta_n? Unfortunately, I'm not seeing that this is true in my code.

Notation: gps, like counter, is one of the items in the value of generate_pseudo_pop(). delta_n is the parameter I input into generate_pseudo_pop().

I'm trying to understand the GPS matching algorithm (described on page 10 of Wu, Xiao, et al. "Matching on generalized propensity scores with continuous exposures." arXiv preprint arXiv:1812.06575 (2018)).

[Naeemkh]

@m-qin, If you use optimized_compile=TRUE, the counter is a vector of the same size as the original dataset, representing the number of times that specific data sample is matched in the process. If you use optimized_compile=FALSE, in that case, it will provide replicated data that the size of the output is equivalent to the sum(counter) if you chose to use the previous approach.

[Naeemkh]

sum(counter) should be the length(counter) * number w bins. But sometimes, it is not the case, because, for some exposures the algorithm cannot find a good match in the provided caliper.

[m-qin]

That makes sense, thank you! No wonder sum(counter)/length(counter) wasn't an integer.
Also, the algorithm is probably using a thinner caliper than the one we provided. (We set delta_n = 0.2, though it looks like the estimated gps values are all between 0 and 0.3.)

[Naeemkh]

If the algorithm cannot find a match for a data sample it will report it in the .log file. Please let me know if you need to see more information in the .log file. Thanks.

[m-qin]

Thank you, Naeem!

[m-qin]

Maybe for future iterations of the CausalGPS documentation, this could be added?

[Naeemkh]

Sure. I will add it to the to-do list. But also feel free to open a new branch and add this documentation wherever fits the best.

Also setting the logger flag to debug can help you see more internal processes.

set_logger("DEBUG")

[m-qin]

Hi Naeem! I have a follow-up question. I'm trying to winsorize values in counter. That is, I'm trying to replace counter values that are larger than the 95th quantile with the 95th quantile of counter.

I'm wondering if it's enough to do something like

winsorize_counter_onesided <- function(counter, quantile){
  cutoff <- quantile(counter, quantile)
  return(ifelse(counter > cutoff, cutoff, counter))
}

matched_pop$pseudo_pop$counter <- winsorize_counter_onesided(matched_pop$pseudo_pop$counter, 0.95)

and then send matched_pop into functions like estimate_pmetric_erf and estimate_npmetric_erf?

I'm asking because I've done this and my results from estimate_pmetric_erf seem to depend on the quantile I winsorize at, as expected, while my results from estimate_npmetric_erf and the covariate balance plot seem to be the same regardless of the quantile I winsorize at, which is strange. For the covariate balance plot at least, I'm guessing I need to modify matched_pop$adjusted_corr_results in some way.

[Naeemkh]

Tagging @kevjosey and @macork.

[kjosey]

Hey @m-qin, is there a test case that you can write that demonstrates this phenomenon so we can try to debug it?

[macork]

Agree with @kevjosey a test case would be useful, also curious when you are running this winsorize function in your code pipeline and agree that the adjusted correlation results would need to be changed

[Naeemkh]

@m-qin, that is correct, and you are right about the covariate balance values. We compute the covariate balance and store it in the pseudo_pop object. When you change the values, it is not going to change those values automatically. The best way to do that is not directly changing the values yourself but rather changing them via a function; inside that function, the values of the covariate balance should be adjusted.
This falls within the trimming discussions, and we need to come up with a universal design to reduce the complexity for the users.
In the meantime, you can compute the covariate balance using check_covariate_balance function.

https://github.com/fasrc/CausalGPS/blob/master/R/check_covar_balance.R

Here is an example:

set.seed(422)
n <- 100
mydata <- generate_syn_data(sample_size=100)
year <- sample(x=c("2001","2002","2003","2004","2005"),size = n, replace = TRUE)
region <- sample(x=c("North", "South", "East", "West"),size = n, replace = TRUE)
mydata$year <- as.factor(year)
mydata$region <- as.factor(region)
mydata$cf5 <- as.factor(mydata$cf5)
pseudo_pop <- generate_pseudo_pop(mydata$Y,
                                  mydata$treat,
                                  mydata[c("cf1","cf2","cf3","cf4","cf5","cf6","year","region")],
                                  ci_appr = "matching",
                                  pred_model = "sl",
                                  gps_model = "non-parametric",
                                  trim_quantiles = c(0.01,0.99),
                                  optimized_compile = TRUE,
                                  sl_lib = c("m_xgboost"),
                                  covar_bl_method = "absolute",
                                  covar_bl_trs = 0.1,
                                  covar_bl_trs_type = "mean",
                                  max_attempt = 1,
                                  matching_fun = "matching_l1",
                                  delta_n = 1,
                                  scale = 0.5,
                                  nthread = 1)


# Modify the counter
# Your code should be here.
# 
# 


# Recompute the corvariate balance
adjusted_corr_obj <- check_covar_balance(pseudo_pop$pseudo_pop,
                                         ci_appr="matching",
                                         nthread=1,
                                         covar_bl_method = "absolute",
                                         covar_bl_trs = 0.1,
                                         covar_bl_trs_type = "mean",
                                         optimized_compile=FALSE)

# Reassign it to the pseudo_pop object
pseudo_pop$adjusted_corr_results <-  adjusted_corr_obj$corr_results

Technical details aside, a higher counter number does not mean an outlier, and I think you should be careful in trimming data at this level. I would dig deep into the data to see why some units are getting higher counter values. Please let me know if we can take a look at data and discuss it more.

Please make sure that you do not copy any L3 data here.

[m-qin]

Thank you all for your replies! @Naeemkh I'll try your suggestion; great point that we should be careful about how we trim/winsorize/stabilize weights!