---
title: "SampcompR for Response Bias Estimation"
author:
  - name: "Rohr Björn"
image: images/logo_sampcompR.png 
date: "`r format(Sys.Date(), '%B %d, %Y')`"
format:
  html:     
    toc: true
    toc-location: left
    toc-depth: 2
    embed-resources: true
bibliography: references.bib
biblio-style: apa
editor: 
  citation: "Rohr, Björn (2025). *SampcompR for Response Bias Estimation*. GitHub Repository. URL:  https://github.com/BjoernRohr/KODAQS-Toolbox-SampcompR"
  markdown: 
    wrap: 72
prefer-html: true

---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  out.width = "100%"
)
```


# At a glance

<!-- Summary with main points of the tool in 2-3 bullet points -->

### Table of Content

[Introduction](#introduction)

[Setup](#setup)

[Tool application](#tool-application)

[Conclusion and recommendations](#conclusion-and-recommendations)

# Introduction {#introduction}

# Set-up

## Data for representation bias analysis

For this tutorial, we use data from an old American survey, called card (n=3010). Although this dataset is rather dated, it is easily available and also suitable do demonstrate the use of the sampcompR package, without much preperation. it consists of several demographic variables as well as some variables to generate groupings, as well as some additional variables.

::: callout-note
Although available in the wooldridge package [@SheaWooldridge] the Card data, was also added to sampcompR and originates from Professor Card (1995):

D. Card (1995), Using Geographic Variation in College Proximity to Estimate the Return to Schooling, in Aspects of Labour Market Behavior: Essays in Honour of John Vanderkamp. Ed. L.N. Christophides, E.K. Grant, and R. Swidinsky, 201-222. Toronto: University of Toronto Press
:::

Some interesting items, as described in the card sample, are:

| Items        | Institution                   |
|--------------|-------------------------------|
| **age**      | in years                      |
| **educ**     | years of schooling, 1976      |
| **fatheduc** | father's schooling            |
| **motheduc** | mother's schooling            |
| **wage**     | hourly wage in cents, 1976    |
| **IQ**       | IQ score                      |
| **married**  | 1 if married, 1976            |
| **KWW**      | knowledge world of work score |
| **south**    | 1 if in south, 1976           |
| **black**    | 1 if black                    |

## Getting started

To use sampcompR [@rohrSampcompRComparingVisualizing2024], we first need to install the package from the repository of CRAN, the Comprehensive R Archive Network. For installation, we can use the following commands:

```{r}
# Installing SampcompR
# install.packages("sampcompR")

# Loading SampcompR
library(sampcompR)


```

In addition, to the the package, we also need data for comparison. Typically a data frame is required, containing the variables for which we want to estimate bias, alongside benchmark information. Benchmark information is best also provided in the form of a data frame, although, for example, if the goal is to estimate bias on an univariate level, a named vector of means (e.g., from a population census) is sufficient. For the purpose of this application example, we perform the bias estimation based on a data frame named `card` in the `wooldridge`[@SheaWooldridge] package.

```{r}
#| warning: false
# Get Data for comparison
data("card")
whole_card<-card

# we recode the married variable to a dummy, idicating of the respondents are married or not, which will be needed later. 
whole_card$married[whole_card$married!=1]<-0

# Split the data frame to generate example data
north<-whole_card[whole_card$south==0,]
white<-whole_card[whole_card$black==0,]

# load some additional packages needed in this application

library(tidyverse)
library(knitr)
library(kableExtra)


```

As described above, for this application example we assume that (1) we have a sample that only contains respondents living in the `North` and want to estimate the bias if this should be generalized to the whole population, and further (2) we have a sample containing only `White` respondents that should be generalized. Although those scenarios are only for the purpose of demonstrating the tool, there can be rather similar scenarios in your own work. For example if you have conducted a mixed mode survey (e.g., online and offline mode) and want to know for future surveys, if it is okay, to only recruit respondents in one of the modes (e.g., online mode).

# Tool application {#tool-application}

## Univariate Comparison

Now that we have data frames we want to compare that are gathered similarly, we will start with a univariate comparison.

```{r univariate_comparison, fig.width= 6, fig.hight=12}
#| warning: false
library(sampcompR)

univar_data<-sampcompR::uni_compare(dfs = c("north","white"),
                                    benchmarks = c("whole_card","whole_card"),
                                    variables=c("age","educ","fatheduc",
                                                "motheduc","wage","IQ"),
                                    funct = "rel_mean",
                                    nboots=0,
                                    summetric="avg2",
                                    data=T,type = "comparison")

sampcompR::plot_uni_compare(univar_data,
                            name_dfs=c("North","White"),
                            name_benchmarks=c("Whole Card","Whole Card"))

```

Here we can see a plot showing the difference in mean for several variables (indicated in the `variables` parameter) in the data frames. The first data frame, named as a character in `dfs`, is compared to the first data frame in `benchmarks`, while the second is compared to the second, and so forth. Another important parameter is `nboots`. It indicates the number of bootstraps to calculate the confidence intervals. We recommend using at least 2000 to 10000 bootstraps, although this may take a while. It is also possible to set `nboots==0`, as we did in this example. In this case, the confidence intervals will be calculated analytically. The upper right corner of the plot shows the average relative difference in mean, which is chosen by using `summetric=="avg2"`. Last, the parameter data indicates if the function should return a `uni_compare_object` which can be used in other package functions (e.g., `plot_uni_compare`).

::: {.callout-note title="Be careful:"}
High or low amounts of bias in some estimates does not necessarily mean that the survey has similar amounts of bias in other variables. This would only be the case, if those variables are highly correlated to the estimates for which bias was estimated. For variables, where no benchmark is available careful consideration is needed. 
:::

Also, the difference can be shown in a table that can be helpful to document the results. We can easily get a good result table by adding the resulting table to other functions (e.g., the kable() function in RMarkdown).

```{r output table}
#| warning: false
uni_output_table<-sampcompR::uni_compare_table(univar_data)


```

**Table 1: Difference in Relative Means off different Survey Groups**

```{r echo=FALSE}
### use the kable function to get a nicer table

library(kableExtra)

table1 <- knitr::kable(
  uni_output_table,
  format = "html",
  booktabs = TRUE,
  col.names = c("Variables", "North", "White"),
  align = c("l", "c", "c")) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(row = 0, align = "c") %>%
  column_spec(1:ncol(uni_output_table), width = "4em") #%>% 
   # footnote(
   #  general_title = "Note.",
   #  general = "Differences between north and white respondents of the card survey, to the respondents of a survey containing the whole population; N is based on the sample size of the north and white sample respectively and may differ due to missing values",
   # footnote_as_chunk = TRUE
   #  )
 
 table1
```


The table displays the difference between dfs and benchmarks for all variables, the confidence intervals, and the summary metric specified when creating the `uni_compare_object` for every comparison. As we can see in the table, as well as in the plot before, surveying only respondents living in the `North` would have let to an overestimation of education, wage and IQ, however the age would have been estimated rather similar. We find similar results for the second comparison. Of cause the data frame for this example is rather dated and a newer data frame might have led to different results.

## Bivariate Comparison

As with the univariate comparison, we start with a general function to calculate the bivariate difference between the data frames and plot the results.

```{r bivariate_comparison, fig.width= 7, fig.hight=12}
#| warning: false
biv_data<-sampcompR::biv_compare(dfs = c("north","white"),
                                 benchmarks = c("whole_card","whole_card"),
                                 variables= c("age","educ","fatheduc",
                                             "motheduc","wage","IQ"),
                                 data=T, corrtype = "rho",
                                 weight = "weight",
                                 id="id")

sampcompR::plot_biv_compare(biv_data,plots_label=c("North","White"))
```

In this comparison, the Pearson's r correlations of all variables indicated in the `variables` parameter are compared between the data frames and benchmarks. The plot shows the difference in three categories. A correlation is green (Same) for one of two reasons. On the one hand, it can be that the Pearson's r correlation of this variable pair is not significant in the data frame and in the benchmark. On the other hand, it can be that both z standardized Pearson's r values are not significantly different. If neither of these conditions is true, the plot's correlations are marked yellow (Small Diff) or red (Large Diff). To be red, it must also fulfill one of two additional conditions. The Pearson's r correlations could be of different directions (e.g., one is positive, while the other is negative). Second, it could be that one of them is double the size of the other.

In our example, we can see on the left that the a survey containing only respondents living in the `North` is very different from the complete survey regarding the investigated bivariate correlations. Only 46.7% of all correlations are similar between both groups (`Same`) and would lead to similar interpretations measured with either group of the survey. 40.0% of correlations show minor differences (`Small Diff`), meaning that similar interpretations may not be guaranteed, while 13.3% of the correlations are very different (`Large Diff`) between the groups. The right half of the plot shows the comparison between `White` respondents and the whole survey, where similar results are found.

As before, we can also take a closer look at the exact values by putting the `biv_compare_object` into the `biv_compare_table` function of our package.

```{r difference matrix}
# Difference correlation matrix of the first comparison
table_biv1<-sampcompR::biv_compare_table(biv_data,type = "diff",comparison_number=1) 

```

**Table 2: Difference in Pearson's r for the North/South Sample**

```{r echo=FALSE}
### Use the kable function to get a nicer table

table2 <- knitr::kable(
  table_biv1,
  format = "html",
  booktabs = TRUE,
  align = c(rep("c",(ncol(table_biv1))))) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(row = 0, align = "c") %>%
  column_spec(1:(ncol(table_biv1)+1), width = "4em") #%>% 
   # footnote(
   #  general_title = "Note.",
   #  general = "Difference in Pearson's r for the respondents in the card Sample, living in the north, from the whole survey; significance tests are based on the difference of z-standardized Person's r values.",
   # footnote_as_chunk = TRUE
   #  )

table2

```

Using `type = diff` gives us a matrix for the difference in Pearson's r between the surveys, while `comparison_number = 1` indicates that the table should be for the first comparison (between the`north` and the whole sample). Here we can again see which correlations significantly differ between the surveys and to what extent. However, to know why the colors are as they are in the plot, we must also look at the individual correlation matrices for both surveys. Here we only look at the tables for comparing `north` versus the whole sample respondents as an example.

```{r correlation matrices of comparison 1}
# North correlation matrix of the first comparison
table_biv2<-sampcompR::biv_compare_table(biv_data,type = "dfs",comparison_number=1) 

# South correlation matrix of the first comparison
table_biv3<-sampcompR::biv_compare_table(biv_data,type = "benchmarks",comparison_number=1) 

```

**Table 3: Pearson's r correlation matrix for the North Sample**

```{r echo=FALSE,results = "asis"}
### use the kable function to get a nicer table

table3 <- knitr::kable(
  table_biv2,
  format = "html",
  booktabs = TRUE,
  align = c( rep("c",(ncol(table_biv2))))) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(row = 0, align = "c") %>%
  column_spec(1:(ncol(table_biv2)+1), width = "4em") #%>% 
   # footnote(
   #  general_title = "Note.",
   #  general = "Pearson's r correlation matrix for the respondents in the card sample, living in the north.",
   # footnote_as_chunk = TRUE,
   # threeparttable= TRUE
   #  )

table3
```

**Table 4: Pearson's r correlation matrix for the Whole Crad Survey**

```{r echo=FALSE}
### use the kable function to get a nicer table

table4 <- knitr::kable(
  table_biv3,
  format = "html",
  booktabs = TRUE,
  align = c( rep("c",(ncol(table_biv3))))) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(row = 0, align = "c") %>%
  column_spec(1:(ncol(table_biv3)+1), width = "4em") #%>% 
   # footnote(
   #  general_title = "Note.",
   #  general = "Pearson's r correlation matrix for the respondents in the card sample, whole card survey.",
   # footnote_as_chunk = TRUE
   # )

table4

```


This shows us the whole picture. Looking at the correlation between wage and fathers education, for example, we can see that the correlations in both surveys are different. While, the correlation is positive and small in the north group, stronger for the in the whole survey. Here, in addition to the conditions for a slight difference, one condition for a `Large Diff` are true (i.e., one Pearson's r value is at least double the size if the other, 0.09 vs. 0.19). Therefore, we could come to very different conclusions measuring those correlations for either group.

## Multivariate Comparison

When we want to know how different certain data frames or sub-data frames are, it may also be interesting to see if those differences exist in multivariate comparisons. For this, the first step is to choose the multivariate models we want to compare. This function `(multi_compare())` is restricted to any model commutable with `glm()` (without weighting) or `svyglm()` (when weights are provided). In this example, we want to see if there are differences between the groups when looking at the regression of `age`, `father's education`, `mother's education`, and `IQ` on wage and education.

After choosing independent and dependent variables, all information can be given into `multi_compare()`. For this comparison, three similar regression models are performed. One model with either survey and one model with the combined surveys, including an interaction variable that indicates what survey the respondents are from. After this, for every coefficient, a similar classification is used, as in the bivariate comparison:

```{r ols_models, fig.width= 8, fig.hight= 7}
#| results: hide
# designate dependent variables 
dependent_ols<-c("wage","educ","KWW")

# designate independent variables
independent<-c("age","fatheduc","motheduc","IQ")


# compare the north and south data frames
multi_data1_ols<-sampcompR::multi_compare(df=north, 
                                     bench=whole_card,
                                     independent = independent,
                                     dependent = dependent_ols,
                                     family = "ols")  

# compare the black and white data frames
set.seed(123)
multi_data2_ols<-sampcompR::multi_compare(df=white, 
                                     bench=whole_card,
                                     independent = independent,
                                     dependent = dependent_ols,
                                     family = "ols")
```

```{r ols_models_plot, fig.width= 8, fig.hight= 7}
# plot the results
sampcompR::plot_multi_compare(c("multi_data1_ols","multi_data2_ols"),
                             plots_label=c("North","White"))

```

The plot looks very similar to that of the bivariate comparison, except that every row shows a coefficient while every column shows a model. Here the green color (Same) indicates that a coefficient is not significantly different in any of the models performed with the individual surveys or that the interaction effect in the model on the pooled data frames is not significant. Similarly, yellow (Small Diff) or red (Large Diff) indicates that the coefficient is significant in one of the base models, and the interaction effect is significant in the pooled models. If those conditions are true, it is red (Large Diff) if the coefficients differ in direction or one is the size of the other and yellow (Small Diff) otherwise.

As we can see here in those models, other than before in the bivariate comparison, there are no differences, at least for those models compared. Also, we can see that the differences on the wage variable we found in previous comparisons are not present in the multivariate model. However, even in multivariate regression, education seems more prone to group differences.

In addition to the models, we can answer with an ols regression, it would also be interesting to add a model estimating the effect of the coefficients on the probability of marriage. However, for this, we need a logit regression, which now follows.

```{r adding_logit_models, fig.width= 7, fig.hight= 7}
#| results: hide
# logit dependent variable


dependent_log<-c("married")


# compare the north and south data frames
multi_data1_log<-sampcompR::multi_compare(df=north, 
                                     bench=whole_card,
                                     independent = independent,
                                     dependent = dependent_log,
                                     family = "logit")  

# compare the black and white data frames
multi_data2_log<-sampcompR::multi_compare(df=white, 
                                     bench=whole_card,
                                     independent = independent,
                                     dependent = dependent_log,
                                     family = "logit")


# merge those  multi_compare_objects to the ols based objects

final_multi1<-sampcompR::multi_compare_merge(multi_data1_ols,multi_data1_log)

final_multi2<-sampcompR::multi_compare_merge(multi_data2_ols,multi_data2_log)

```

```{r adding_logit_models_plot, fig.width= 8, fig.hight= 7}
# plot the models together
sampcompR::plot_multi_compare(c("final_multi1","final_multi2"),
                             plots_label=c("North","White"))
```

After performing the same comparison with the `method = logit`, and with the new dependent variable, it would still be interesting to plot all models together. This can be done with the function `multi_compare_merge`, which adds multi_compare_objects together and use the newly created objects for the plot. This shows us that in addition to the education model, in the marriage model, again no coefficient is different for white respondents from the same models conducted on the whole survey.

We also add the objects into the provided table function to look at the differences closer.

```{r multi_compare_table}
### difference table
multi_table1 <- sampcompR::multi_compare_table(c("final_multi1","final_multi2"),type="diff")
multi_table2 <-sampcompR::multi_compare_table(c("final_multi1","final_multi2"),type="dfs")
multi_table3 <-sampcompR::multi_compare_table(c("final_multi1","final_multi2"),type="benchmarks")

```

**Table 5: Multivariate Comparison of Subgroups of the Card Sample : Interaction Models**

```{r echo=FALSE}
### Use the kable function to get a nicer table

table5 <- knitr::kable(
  multi_table1,
  format = "html",
  booktabs = TRUE,
  align = c("l",rep("c",(ncol(multi_table1)-1)))) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(row = 0, align = "c") %>%
  column_spec(1:ncol(multi_table1), width = "4em") #%>%
   # footnote(
   #  general_title = "Note.",
   #  general = "Difference in regression coefficients for three multivariate regression models with the card sample, respondents living in the north were compared to the whole card survey; white respondents were compared against the whole card survey; results are based on an interaction model.",
   # footnote_as_chunk = TRUE
   #  )

table5

```

**Table 6: Multivariate Comparison of Subgroups of the Card Sample, North and White Sample Models**

```{r echo=FALSE}
### Use the kable function to get a nicer table

table6 <- knitr::kable(
  multi_table2,
  format = "html",
  booktabs = TRUE,
  align = c("l",rep("c",(ncol(multi_table2)-1)))) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(row = 0, align = "c") %>%
  column_spec(1:ncol(multi_table2), width = "4em") #%>%
   # footnote(
   #  general_title = "Note.",
   #  general = "Difference in regression coefficients for for three multivariate regression models with the card sample; results are for respondents living in the north, and white respondents.",
   # footnote_as_chunk = TRUE
   #  )

table6

```

**Table 7: Multivariate Comparison of Subgroups of the Card Sample, Whole Sample Models**

```{r, echo=FALSE}
### Use the kable function to get a nicer table
multi_table3<-multi_table3[1:8,]
table7<-knitr::kable(
  multi_table3,
  format = "html",
  booktabs = TRUE,
  align = c("l",rep("c",(ncol(multi_table3)-1)))) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(row = 0, align = "c") %>%
  column_spec(1:ncol(multi_table3), width = "4em") #%>%
   # footnote(
   #  general_title = "Note.",
   #  general = "Difference in regression coefficients for for three multivariate regression models with the card sample; results are for respondents living in the south, and white respondents.",
   # footnote_as_chunk = TRUE
   #  )

table7


```


Overall the comparison shows no differences between the sub-groups and the whole survey. Suppose those models play a role in the overall analyses of a research project, one would come to similar results using any survey. In other words, for those multivariate models, we would not find any bias, had we only sampled parts of the population. Nonetheless, there could be other models, where bias can be prevalent. <!--subheadings for tool application can be specified "freely" and depending on structure of your tool -->

# Conclusion and recommendations {#conclusion-and-recommendations}

**Univariate Bias**
In our exemplary comparison we found that all variables of comparison, except for age would be significantly biased, if one would generalize the estimates, of respondents living in the north to the whole population. Similar amounts of bias would be encountered, when generalizing the results of white respondents to the whole population.

**Bivariate Bias** Similarly, we also find bias in Pearson's r estimates, when generalizing the results to the whole population. However, we also see that many of the differences are rather small, and in most of the cases, one would come to similar conclusions using any survey. Nonetheless, as shown for the north survey, one often would underestimate some correlations, especially regarding wage, when one generalize this survey to the whole population.  

**Multivariate Bias**
Finally, for the multivariate models, compared one would get similar results using only the north or white survey, instead of the whole survey. The reason for this unbiased, even in the wage model might be that controlling for the demographic variables, such as age, IQ and parental education, the bias in the separate correlations is reduced. 

**Other use cases**
Of cause, the example cases in this tool-description, are only for demonstration purpose, however there are very similar examples, that one could encounter in the real world of survey science. Taking for example a situation where one has gathered a mixed mode panel survey, with an inexpensive online and an expensive offline mode. Then it might be interesting to know how much bias one would introduce into the survey, when, in further waves only continuing the online mode. Here a similar bias estimation could be worth it. Other applications would be the comparison between a full sample and the actually recruited respondents, to evaluate nonresponse bias, or comparing a freshly recruited panel against the same panel, after several years, to evaluate the effect attrition had on the panel. When external benchmark data is available, one could also compare several recruitment methods, such as a probability and a non-probability survey, to evaluate if both surveys lead to similar amounts of bias or not [i.e., @rohrComparingAccuracyUnivariate2024a]. 

**Recommendation:**
As acknowledged by the Total Survey Error Framework [@grovesTotalSurveyError2010], any survey estimate has potential for. While often this potential bias is only mentioned as a limitation of the study, often one could go a step further.

- Although it is seldom possible to measure the bias for all variables of interest, due to a lack of benchmarks, it might be possible to perform at least some bias analyses, e.g., for demographic variables. Reporting such findings as additional analyses in the Appendix would not only indicate transparancy in research, but might also take some weight of limitations. Although unbiased estimates for one variable do not guarantee for other variables to be unbiased, they might at least hint at it, especially when one can reasonably argue that the variables are correlated. 

- Although we have shown only a simple use case for this tool, SampcompR is also applicable to more complex data, as it can also make use of strata, cluster information, and weights. Further raking and post-stratification weights can be calculated during comparison, and bootstrap can be used to evaluate the difference between survey and benchmark. The later method is often recommended for complex surveys including nonprobability surveys [i.e., @mcpheeDataQualityMetrics2022]

# References

