Biased Bivariate Correlations in Combined Survey Data Measured With Different Instruments

Mean Bias

A substantive problem caused by incomparable and insufficiently harmonized measurement units is mean bias (Kolen & Brennan, 2014; Singh, 2022). Mean bias is best understood when we consider the following thought experiment. Imagine applying two instruments X and Y to sufficiently large random samples A and B of the same population. If ex-post harmonization was successful, we would expect the mean responses to be approximately equal: $\stackrel{-}{X}\approx \stackrel{-}{Y}$. This is because the average true score should be the approximately the same in two random samples of the same population: ${\stackrel{-}{\tau }}_A\approx {\stackrel{-}{\tau }}_B$ (Kolen & Brennan, 2014; Singh, 2020). However, without adequate ex-post harmonization, we might find that the average response to two congeneric instruments for the same construct differs by some constant d: $\stackrel{-}{X}=\stackrel{-}{Y}+d$ (Price, 2017; Raykov & Marcoulides, 2011). In other words, combining data across the two instruments introduces a mean bias $d$. This can easily occur, if two instruments for the same construct have different item difficulties (Moosbrugger & Kelava, 2012). Respondents may find one instrument wording easier to agree to then the alternative and thus for the same population of respondents would choose a higher mean response on one instrument than the other.

In practical terms this means that even after (insufficient) harmonization, an average respondent for our measured population would be represented by different numerical scores in our combined data (Kolen & Brennan, 2014). Of course, mean bias can also occur if instruments X and Y are applied to different populations. However, in such cases we cannot easily isolate the bias for single-question instruments, because the difference between $\stackrel{-}{X}$ and $\stackrel{-}{Y}$ is then a composite of the true construct difference ${\stackrel{-}{\tau }}_A-{\stackrel{-}{\tau }}_B$ and the bias. Again, in practical terms, this means the mean differences between the populations are either over- or underestimated by an amount proportional to the mean bias (Kolen & Brennan, 2014).

For a concrete example, consider two very similar measures of political interest. In Germany, the International Social Survey Programme (ISSP) is fielded together with the German General Social Survey (ALLBUS). In 2014, both asked the respondents about their level of political interest. The ALLBUS asked “How strongly are you interested in politics” with a five-point scale (GESIS-Leibniz-Institut Für Sozialwissenschaften, 2018), and the ISSP asked “How interested would you say you are in politics?” with a four-point scale (ISSP Research Group, 2016). Due to the different number of response options alone, we would expect a different mean response. And indeed, the mean response differed by Cohen’s $\left|d\right|=0.69$ (Singh, 2022). To align these differences in scale points (and thus scale range), harmonization practitioners then often apply linear stretching (Cohen et al., 1999; de Jonge et al., 2017; Durand et al., 2021). This means that the scale ranges (i.e., maximum score minus minimum score) are aligned, by setting the minimum responses and the maximum responses as equal across instruments and then stretching all scores equidistantly in between. In our example, we would linearly stretch the scores 1, 2, 3, 4, of the four-point ISSP instrument to 1, 2.33, 3.67, and 5. However, after this transformation we still find that the average responses differed by $\left|d\right|=0.38$ between the two instruments (Singh, 2022). As it turns out, the two instruments differed in more than their scale range. Additionally, we find that both instruments have different item difficulties of $P=43$ for the ALLBUS instrument and $P=33$ for the ISSP instrument (Singh, 2022). In other words, the average respondent chose a score that was 43% along the range from 1 to 5 in the ALLBUS instrument but chose a score that was only 33% along the range from 1 to 4 in the ISSP instrument (Moosbrugger & Kelava, 2012; Singh, 2022). However, linear stretching only aligns the scale ranges but not the position of the average response within the scale range. Thus, differences in item difficulty between two instruments remain untouched and can cause a substantive remaining mean bias when data from the two instruments are combined.

Such differences in difficulty between single-question instruments can be mitigated if we apply more suitable harmonization method than linear stretching. One example is observed score equating in a random groups design: OSE-RG (Kolen & Brennan, 2014; Singh, 2022). Alternatively, such difficulty differences can also be mitigated via multiple imputation. However, in this paper, we want to explore what happens if we fail to mitigate mean bias. Or in terms of our example, what would have happened if we had only used linear stretching. After all, harmonization practitioners might be unaware of the limitations of linear stretching. Or they may find the more suitable harmonization techniques unfeasible in their projects. After all, both approaches have data requirements that are not always easy to meet. OSE-RG, for example, requires random groups data; that is samples for both instruments drawn randomly from the same population. Harmonizing two instruments for the same construct with multiple imputation, meanwhile, requires a calibration sample (Siddique et al., 2015): That is a sample in which each respondent answered both instruments, but in a way that does not lead to question order effects.

The first consequence of mean bias in composite data is straightforward: Under mean bias, scores derived from different instruments are biased by an additive constant. In data drawn from the same population, this introduces error variance. And even more worrisome, in data from different populations measured with different instruments, the mean bias is mingled with true population differences. Thus, we can no longer be certain that if two populations differ on average, that this is a true population difference. Instead, it could be a methodological artifact.

Biased Bivariate Correlations due to Mean Bias

In this paper, however, I will demonstrate with systematic simulations that bivariate correlations between two variables in a harmonized dataset can be biased as well, if one or both variables are subject to mean bias. This might seem surprising, because Pearson product-moment correlations are unaffected by linear transformations of variable scores. Specifically, an additive constant $d$ would not change the correlation coefficient $r$, because adding a constant to each score changes the arithmetic mean with the same constant (adapted from Gill, 2008):

$\frac{1}{n}{\sum }_{i=1}^n(x_i+d)=\stackrel{-}{x}+d$

Thus, the Pearson product moment correlation formula (Gill, 2008) remains unchanged by an additive constant applied to all values of x (or y):

$r_{xy}=\frac{{\sum }_{i=1}^n\mathrm{ }[\left(x_i+d\right)-\left(\stackrel{-}{x}+d\right)](y_i-\stackrel{-}{y})}{s_x{s}_y}=\frac{{\sum }_{i=1}^n(x_i-\stackrel{-}{x})(y_i-\stackrel{-}{y})}{s_x{s}_y}$

However, this intuition is misleading, because in harmonization we combine data from different sources. In our combined dataset, a vector of responses for a construct would thus be composed of scores derived from different instruments. Consequently, mean bias does not add a constant to the whole variable, as above. Instead, mean bias adds different additive biases to different segments of the combined variable. The argument above obviously no longer holds if we add a constant d to some $x_i$, but not all.

Imagine the following simplified combined dataset with data from surveys A and B. Both surveys measured the constructs P and Q. However, both surveys used different instruments for each construct, thus leading to a total of four instruments. If we combine these data, we arrive at the following combined dataset structure: a dataset with two variables, one for construct P and one for construct Q. Crucially, each construct variable (i.e., vector) is a composite of values from two surveys and thus two instruments. In summary, there are thus three vectors for construct P: $P_A$ from survey A, $P_B$ from survey B, and $P_C$ as the composite vector in the joint dataset with data from two different instruments combined. Analogously, for construct Q there are then the vectors $Q_A$, $Q_B$, and $Q_C$.

This composite structure of variables in the combined dataset is crucial for understanding the impact of mean bias. After all, mean bias introduced by instrument differences does not affect the whole composite variable. Instead, it would be as if we added a constant to only half of the composite variable. Figure 1 shows the data structure schematically on the left.

Click to enlarge

Figure 1

Schematic Overview of a Composite Dataset and the Resulting Simpson’s Paradox

In this example, we assume that the scale of instrument $P_B$ assigns scores that are on average $d_P=-2$ lower than the scores that instrument $P_A$ would assign for respondents with the same true score in construct P. Instrument $Q_B$, meanwhile, assigns scores that are on average $d_Q=2$ higher than the scores that instrument $Q_A$ would assign for respondents with the same true score in construct Q.

If we now plot this combined dataset but differentiate by source surveys (and thus source instruments) we observe a surprising pattern in Figure 1 on the right. Within the data from each survey, constructs P and Q are perfectly correlated with $r_A=r_B=1.0$ (pink and green trend lines). However, if we correlate $P_C$ and $Q_C$ across the combined data, the correlation drops to a mere $r_C=.11$ (blue trend line).

What has happened? The conditional correlations in two groups are identical, but the correlation across both groups is substantially different. Through the mean bias in P and in Q, we have introduced a version Simpson’s paradox into our combined data (Rücker & Schumacher, 2008). In general, Simpson’s paradox describes an empirical pattern where we observe the same relationship between two variables in two (or more) groups separately, but a substantially different relationship in analyses across the groups. Of course, Simpson’s paradox is not an actual paradox. Instead it is “a form of bias, resulting from heterogeneity in the data [that has not been] accounted for” (Rücker & Schumacher, 2008). If we decompose the problem, we see that the overall correlation $r_C$ is a composite of the individual source survey (and thus source instrument) correlations, $r_A$ and $r_B$, on the one hand, and of the spurious correlation introduced by the mean biases ($d_P$ and $d_Q$ in red). To visualize this mean biased induced competing correlation as a trendline through the two bivariate group means in surveys A and B (the plus signs on the trend lines of surveys A and B). In other words, if we do not account for mean bias, ideally by removing it with a suitable harmonization procedure, we bias correlations by $\mathrm{\Delta }r=r_C-r_u$.

This simple example showed that, in principle, mean bias between source instruments can result in biased correlations. However, under what circumstances does this correlation bias occur and with which intensity? In this paper, I set out to map the landscape of this bias with a series of systematic simulations. Specifically, the simulations will demonstrate that the correlation bias due to mean bias in composite data is determined by the interaction of three factors. First, by the mean bias in the composite variable $P_C$. Graphically, in Figure 1 above, we would shift the data of survey B left or right. Second, by the mean bias in the composite variable $Q_C$. Graphically, we would shift the data of survey B up or down. Third, by the strength of the unbiased correlation $r_u$ between the constructs P and Q. By unbiased correlation, I mean the correlation we would expect in the absence of mean biases: $r_u=r_C|d_P=d_Q=0$. Graphically, lower or higher values of $r_u$ would mean that the datapoints of A and B would fluctuate more or less around the diagonal trendlines for each survey.

By systematically varying all plausible combinations of those three factors, we can map out plausible boundaries of the mean bias induced correlation bias $\mathrm{\Delta }r$. With these simulations, I aim to provide practical insights into the following questions:

1. How large is the maximum potential bias for plausible mean biases $-1\le d\le 1$?

2. How do different combinations of mean biases $d_P$ and $d_Q$ impact the correlation bias?

3. How different unbiased correlations $r_u$ impact the correlation bias?

4. Under which conditions are absolute empirical effect size $\left|r_C\right|$ over- or underestimations of the unbiased absolute effect size $\left|r_u\right|$?

5. Can mean bias cause correlations to change direction?

The overarching goal of this paper is to clearly map out the extent and shape of a bias in survey data harmonization that practitioners might not have previously considered. Armed with this intuition, harmonization practitioners can better anticipate the risk of incurring a correlation bias in their analyses if substantive variables composed of different source variables are not adequately harmonized.

Software

All simulations, analyses, and plots were done in R (R Core Team, 2021) within the RStudio IDE (RStudio Team, 2022). The tidyverse package collection (Wickham et al., 2019) was used for data transformation, automation, and data visualization. The pairs of correlated variables were simulated using the faux package (DeBruine, 2021).

Simulation

To answer the research questions, a three-dimensional matrix of simulations was computed. The matrix thus contains simulated estimates of $\mathrm{\Delta }r$ for different bias configurations: i.e., combinations of plausible values of the mean biases $d_P$, $d_Q \mathrm{i}\mathrm{n}\mathrm{ }\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{ }\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{ }\mathrm{v}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{ }$and unbiased correlations $r_u$. This bias configuration matrix allows us to systematically map out the resulting correlation biases for all combinations of mean biases and unbiased bivariate correlations. This is crucial because, as we will see in the results, the correlation bias changes drastically depending on these three factors. In the following, I first describe the basic assumptions, the process of generating a single correlation simulation for a single bias configuration, and then the structure of the whole bias configuration matrix.

Simulation Parameters

To clearly isolate the effects of mean bias, the simulation uses continuous, standard normally distributed pairs of variables (i.e., $\stackrel{-}{x}=0, s=1$), and a predefined covariance and thus a predefined correlation $r_u$. These variables were randomly generated using the faux package (DeBruine, 2021). Each simulated vector had a length of 10,000 elements. Since each vector was then duplicated (see next section), each simulation created a combined dataset C with two variables and a sample size $N=\mathrm{20,000}$. Mean bias is introduced by adding a constant to every value of a variable. This creates pure mean bias in the sense that the mean is the only distribution parameter that changes. Furthermore, since all source variables have the same standard deviation of $s=1$, raw mean differences can be directly interpreted as Cohen’s $d$ values.

Simulated Mean and Correlation Bias

For this paper, I simulate many bivariate correlations with different parameters. Each correlation is determined by three parameters. First, the mean bias $d_P$ in the composite variable $P_C$, second, the mean bias $d_Q$ in the composite variable Q, and third the unbiased correlation coefficient between constructs P and Q, $r_u$. Based on those three parameters, a simple harmonized dataset is simulated analogous to the one presented in Figure 1.

The algorithm works like this. First, two vectors of simulated data are created. Both vectors contain standard normally distributed continuous values. The vectors exhibit a bivariate correlation that reflects the unbiased correlation $r_u$ that we aim for in this specific simulation. This setup represents data from survey A with variables $P_A$ and $Q_A$, measured by survey A’s instruments for the constructs P and Q. Second, we duplicate the simulated vectors for survey A and modify it with the mean bias parameters $d_P$ and $d_Q$. This creates a simulated survey B with variables $P_B$ and $Q_B$, that is identical to survey A except the different measurement unit: For example, a respondent who chose a response $x$ in $P_A$ would have chosen a response $x+d_P$ in $P_B$. This also ensures that the correlations in each survey are the same and equal to the unbiased correlation we aim for: $r_A=r_B=r_u$. Third, we combine the simulated data for survey A and survey B together to generate a simulated harmonized dataset C. This means we append $P_B$ to $P_A$ to form $P_C$ and we append $Q_B$ and $Q_A$ to form $Q_C$. Figure 2 below summarizes this process.

Click to enlarge

Figure 2

Schematic Overview of the Simulation Process for a Specific Combination of $d_P,\mathrm{ }{d}_Q,$ and $r_u$

Based on these combined vectors $P_C$ and $Q_C$, we can calculate the biased correlation $r_C$. The correlation bias $\mathrm{\Delta }r$ can then be calculated as the difference between the correlation biased by mean bias in P and Q and the unbiased correlation: $\mathrm{\Delta }r=r_C-r_u$. The correlation bias $\mathrm{\Delta }r$ can be interpreted as follows: A positive $\mathrm{\Delta }r$ values mean that the biased correlation is numerically higher than the unbiased correlation and negative $\mathrm{\Delta }r$ values mean that the biased correlation is numerically lower than the unbiased correlation. Please note that this is not the same as the absolute correlation effect size being stronger or weaker. For example, a positive $\mathrm{\Delta }r$ in the context of a strong negative unbiased correlation means that the correlation is weaker (e.g., $r_u=-1; r_C=-0.6;\mathrm{\Delta }r=0.4)$.

Correlation Bias as a Function of the Combined Bivariate Mean Bias Strength

Our first question was: How large is the maximum potential bias for plausible mean biases $-1\le d\le 1$? To answer this, we can calculate a measure of bivariate mean bias. Specifically, I calculated the Euclidian distance between two points defined by the means of the two variables in survey A and B. In Figure 2 above, this would be the distance between the green and pink plus-signs signifying the bivariate means in surveys A and B. Since the bivariate mean biases are two-dimensional, we can apply the Pythagorean theorem (Gill, 2008) to calculate the Euclidean distance as follows:

If both composite variables are subject to a mean bias of $d_P=d_Q=1$, then the distance would be $\sqrt{2}=1.41$, for example. Meanwhile, if one mean bias is zero, then the distance is equal to the other mean bias. And, of course, all other combinations of mean biases work as well. For example, if $d_P=0.5$ and $d_Q=1$ then the distance is 1.12. Figure 3 below now shows all simulations plotted by their bivariate mean distance (x-axis) and their correlation bias $\mathrm{\Delta }r$ (y-axis). For easier interpretation, I have added two x-axis scales. The bottom scale shows the raw distance measure. The top x-axis, meanwhile, gives an example for the special case where both mean biases have the same absolute value. So a distance of 0.71 can, for example, be the result of the mean biases $d_P=d_Q=0.5$.

Click to enlarge

Figure 3

Correlation Bias as a Function of Bivariate Mean Bias

The graph illustrates that the range of possible correlation biases increases as mean biases increase. Specifically, the range of possible correlation biases increases quadratically. To show this, I have selected only the highest correlation biases for each distance and fitted a linear model which regressed the correlation bias on distance and quadratic distance. The black trendline shows the result. The numbers represent the maximum positive correlation bias at selected distances. If both variables have a mean bias of $d=0.5$, then we would expect a bias range of $-.12\le \mathrm{\Delta }r\le .12$, or in other words a span of .22. If both mean biases are $d=1$, then this range increases to $-.4\le \mathrm{\Delta }r\le .4$, or a span of .8.

At the same time, the graph intuitively shows that the correlation bias depends on more than just the distance. In fact, even for the highest distance, there are still cases without any mean bias. This is because correlation bias depends on the direction of the mean bias in relation to the direction of the unbiased correlation $r_u$ as well as on the strength of the unbiased correlation ${|r}_u|$. In the following sections, we will unravel these interactions step by step.

Correlation Bias as a Function of the Two Mean Biases Separately

To get a feeling for the interaction between the mean biases $d_P$ and $d_Q$ in both variables as well as the unbiased correlation $r_u$, let us consider a series of mean bias grids in Figure 4. Each grid varies mean biases systematically along the x and y axes. The colors meanwhile indicate the correlation bias for each mean bias combination. The five different panels, meanwhile, show the correlation bias pattern for a different underlying unbiased correlation $r_u$. From left to right, $r_u=-1; -.5;0;.5;1$. We can thus interpret each panel as a two-dimensional cross-section of the three-dimensional simulation matrix for a given value of $r_u$.

Click to enlarge

Figure 4

Correlation Bias as a Function of $d_P$ and $d_Q$ for Selected Values of $r_u$

Let us first consider the panel in the middle. Here, the two constructs are uncorrelated with $r_u=0$. However, mean bias introduces spurious correlations. Since correlation bias is calculated as ${\mathrm{\Delta }r=r_C-r}_u$, we see that if the mean biases are both positive or both negative, then we create a spurious positive correlation. In the extreme cases (the upper-right and lower-left corners), we see that $d_P=d_Q=1$ and $d_P=d_Q=-1$ result in a spurious correlation of $r_C=.2$. In the opposite case, where $d_P=-1$ and $d_Q=1$ or vice versa, we see a negative spurious correlation of $r_C=-.2$. In other words, the positive and negative mean bias diagonals are especially prone to creating correlation biases. In contrast, if either $d_P=0$ or $d_Q=0$, then no amount of mean bias creates a correlation bias. However, as we will see, this only remains true if $r_u=0$.

Next, we consider edge cases where the constructs are perfectly correlated. First, the rightmost panel, where the unbiased correlation is perfectly positive with $r_u=1$. Here we observe the strongest bias if the mean biases align diametrical to the correlation direction. In other words, we observe the strongest bias if $d_P=1$ and $d_Q=-1$, or when $d_P=-1$ and $d_Q=1$. At the same time, we observe a positive diagonal corridor that is completely free of correlation bias. If the bivariate mean bias moves along the positive diagonal, it merely aligns with the perfectly positive correlation. Second, leftmost panel, we see the same pattern, only mirrored. Bivariate mean bias along the positive diagonal creates the strongest correlation bias, bivariate mean bias along the negative diagonal results in no bias at all. We also find the same pattern as in the previous section. In the case of perfect correlations, we can realize positive or negative correlation biases of $\left|\mathrm{\Delta }r\right|=.4$. In both cases, the bias means that the perfect correlation is reduced in strength.

Now consider the second and fourth bias grid, with more realistic unbiased correlations of ${|r}_u|=.5$. Here, the problem of mean bias creating a bivariate correlation bias reveals its full extent. The only mean bias configuration with a correlation bias of zero is the middle devoid of mean bias ($d_P=d_Q=0$). Mean biases that move in the opposite direction as the unbiased correlation still results in stronger bias than mean biases that move in the same direction as the unbiased correlation. In these cases, the unbiased correlation is severely underestimated ($|\mathrm{\Delta }$r|=.3). However, mean biases that move in the same direction as the unbiased correlation also led to bias; this time a bias that overestimates the correlation by $\left|\mathrm{\Delta }\mathrm{r}\right|=0.1.$

Lastly, we see that the correlation bias intensity seems to vary for different levels of unbiased correlations. However, this is not quite true. In fact, the overall correlation bias range remains a steady $\max \left(\mathrm{\Delta }r\right)-\min \left(\mathrm{\Delta }r\right)\approx .4$ along the whole range of unbiased correlations $-1\le r_u\le 1$. It is just that the direction of the correlation bias shifts. Negative perfect unbiased correlations $r_u=-1$ only allow for positive bias values, positive perfect unbiased correlations $r_u=1$ only allow for negative bias values, while an unbiased correlation of $r_u=0$ allows for symmetrical biases in both directions. All other unbiased correlation values in between follow a linear pattern defined by these three selected cases. In the Supplementary Materials, I have plotted this dynamic in greater detail in Supplementary Figure A (see Singh, 2024b).

Is the Correlation Under- or Overestimated?

Now that we have gained a better understanding of the correlation bias dynamic, we can address two issues of practical relevance. First, practitioners might wonder which mean bias configurations inflate or reduce absolute correlations. In other words, whether absolute empirical effect sizes overestimate ($\left|r_C\right|<\left|r_u\right|$) or underestimate ($\left|r_C\right|<\left|r_u\right|$) the absolute unbiased effect size. In Figure 5 below we see the familiar mean bias grids, but this time colors represent areas of reduced effect sizes in blue, inflated effect sizes in yellow, and areas with negligible bias ($\mathrm{\Delta }r<.01$) in green. The unbiased $r_u$ values range from zero to one, because negative values are just mirrored along the x-axis. We see that when $r_u=0$, then every corner inflates absolute correlations. However, as soon as $r_u$ has a direction, then we have inflation in sectors where the mean bias direction aligns with the correlation direction, and reduction in sectors where the mean bias direction is opposite to the correlation direction. There are two further patterns worth noting. First, as $r_u$ increases in strength, the inflated sector contracts, the reduction sector expands. Second, as $r_u$ increases, bias that reduces the effect size gains in intensity and bias that inflates the effect size loses intensity.

Click to enlarge

Figure 5

Reduced or Inflated Absolute Effect Sizes

Can Mean Bias Change the Correlation Direction?

A last aspect of practical importance is if mean bias can invert the direction of correlations. Some researchers are more interested in the direction of effects then the absolute effect size. In such cases, it would be fatal if mean bias changed the direction of an effect. And indeed, there are simulations where the direction of $r_u$ is inverted. The minimum absolute distance was 0.57. In other words, in the worst case, an inverted correlation can already occur if both variables are biased with $d=0.4$. However, such a change in direction can only occur for unbiased r values ${|r}_u|<.25$ as long as the mean biases remain between $-1\le d\le 1$. Supplementary Figure B shows the mean bias areas with inverted correlations in detail (see Singh, 2024b).

Conclusion

The paper has two main practical implications: (1) Wherever possible, mean bias should be mitigated by applying appropriate ex-post harmonization procedures, such as observed score equating (Singh, 2022) or multiple imputation using a calibration sample (Siddique et al., 2015). (2) Where mean biases cannot be mitigated, correlations based on multi-source data should be interpreted with caution: The direction of small correlations may have been inverted and comparisons of the relative correlation strength across different instruments might be misleading.