# Comparing multiple groups with shared data

I am trying to assess for differences between groups where a participant may appear in one or multiple groups.

The aim is to assess if an outcome is worse in people with an abdominal injury, limb or thoracic injury. However, participants may have a limb and a thoracic, or a abdominal and a limb etc. I have thought about placing in a ling format, whereby the pathology (abd, thor, limb) is under one column. This means that those with multiple will appear twice. I can say the mean/median for those who have an abdominal injury (for instance) as PART of their injury pattern, but am wondering how to compare these to those with a thoracic injury as PART of their injury pattern as there will be some who have both.

Many thanks

## 3 Answers

With overlapping groups, you can still use a t-test based on the following:

\begin{aligned} \text{Under }H_0:\space \bar{Y}_1 = \bar{Y}_2 &\\ \\ \frac{\hat{\bar{Y}}_1 - \hat{\bar{Y}}_2}{\sqrt{\operatorname{var}(\hat{\bar{Y}}_1 - \hat{\bar{Y}}_2)}} &\sim t(df) \end{aligned}

The tricky part is that in order to estimate the variance term in this test statistic, you have to make sure to take into account the correlation between $$\hat{\bar{Y}}_1$$ and $$\hat{\bar{Y}}_2$$ caused by the overlap. Also, it's not clear what you should use for the degrees of freedom, $$df$$. For a simple random sample, you could perhaps just use the sample size minus 1. For a complex sampling design, you could probably do well to just use the rule-of-thumb design degrees of freedom $$(\#\text{ of clusters} - \#\text{ of strata})$$.

Below I present a couple options for getting the estimated standard error of the difference using the survey package in R. Note though that you can use it even if your data aren't a survey.

## Option #1: Use replicate weights

You can do this easily using replication methods (jackknife, bootstrap, etc.) if you have sets of replicate weights: simply use each set of replicate weights to calculate the two estimates and their difference, and then calculate the variance of the estimated difference across the sets of replicate weights, using the appropriate method for how your replicate weights were formed. If you're working with the survey R package, you can use the withReplicates() function to do this.

library(survey)
library(srvyr)

# Load survey of schools ----
data(api, package = "survey")

# Create survey design object with replicate weights ----
dstrata <- apistrat %>%
as_survey(strata = stype, weights = pw) %>%
# Generate two overlapping groups
mutate(meals_0_to_60 = ifelse(meals >= 0 & meals <= 60,
1, 0),
meals_50_to_100 = ifelse(meals >= 50 & meals <= 100,
1, 0))

dstrata_rep <- as_survey_rep(dstrata)

# Use the withReplicates() function to estimate standard error of differnces ----
diff_in_means <- function(svy_weights, svy_vbls) {
ybar_1_num = sum(svy_weights * svy_vbls[['meals_50_to_100']] * svy_vbls[['api00']])
ybar_1_dem = sum(svy_weights * svy_vbls[['meals_50_to_100']])
ybar_1 <- ybar_1_num / ybar_1_dem

ybar_2_num = sum(svy_weights * svy_vbls[['meals_0_to_60']] * svy_vbls[['api00']])
ybar_2_dem = sum(svy_weights * svy_vbls[['meals_0_to_60']])
ybar_2 <- ybar_2_num / ybar_2_dem

ybar_1 - ybar_2
}

withReplicates(dstrata_rep, theta = diff_in_means)
#>        theta    SE
#> [1,] -159.18 13.11


## Option #2: Use the method of linearization by influence functions

You can also use linearization methods to get this variance estimate. Here's a blog post explaining how to do this with the survey R package and why it works.

https://www.practicalsignificance.com/posts/how-correlated-are-survey-estimates-from-overlapping-groups/

An example with the school survey dataset from the example code is as follows:

library(survey)
library(srvyr)

# Load survey of schools ----
data(api, package = "survey")

# Create survey design object ----
dstrata <- apistrat %>%
as_survey(strata = stype, weights = pw) %>%
# Generate two overlapping groups
mutate(meals_0_to_60 = ifelse(meals >= 0 & meals <= 60,
1, 0),
meals_50_to_100 = ifelse(meals >= 50 & meals <= 100,
1, 0))

# Calculate means for the two groups ----

ybar_1 <- svyratio(numerator = ~ I(meals_50_to_100 * api00),
denominator = ~ meals_50_to_100,
design = dstrata, influence = TRUE)

ybar_2 <- svyratio(numerator = ~ I(meals_0_to_60 * api00),
denominator = ~ meals_0_to_60,
design = dstrata, influence = TRUE)

diff_in_means <- ybar_1$$ratio - ybar_2$$ratio

# Obtain influence functions of the two statistics ----

dstrata$$variables[['ybar_1_infl']] <- attr(ybar_1, 'influence') * dstrata$$prob
dstrata$$variables[['ybar_2_infl']] <- attr(ybar_2, 'influence') * dstrata$$prob

# Use the influence functions to estimate variances and covariance ----

infl_totals <- svytotal(x = ~ ybar_1_infl + ybar_2_infl,
design = dstrata)

vcov_of_ybar_estimates <- vcov(infl_totals)

# Estimate standard error of difference ----

var_of_diff <- t(c(1,-1)) %*% vcov_of_ybar_estimates %*% c(1,-1)
std_err_of_diff <- sqrt(var_of_diff)

#>          [,1]
#> [1,] 13.02953


I'm not sure this is correct... I would fit an ordinary linear model (or a generalized linear model, if more appropriate) using outcome as response variable and patient + injury as explanatory variables. That is, your input data would have three columns (outcome, patient, and injury) with "patient" appearing more than once if more injuries are present. Then you can compare different injuries accounting for the fact that the same patient could have multiple injuries.

This would be the syntax in R using the emmeans package:

fit <- lm(data=dat, outcome ~ patient + injury)
emm <- emmeans(fit, spec=~injury)


Get the estimates of all pairwise contrasts, e.g. "abd vs thor", "abd vs limb", etc:

contrast(emm, method='pairwise')


With this setup a patient with multiple injuries would have the same "outcome" acorss the same injuries (if I understand the question correctly), e.g.:

patient injury outcome
A       abd    1.234
A       limb   1.234


and I'm not sure if this is right.

• Hi, thanks for your response. My outcomes are either binary or ordinal so would this method still be appropriate for a glm or polr? Oct 26, 2022 at 12:17
• Yes, I think so, in principle the method is valid (conditional to my uncertainty in the question!) Oct 26, 2022 at 12:29

One way, probably not the most sophisticated one but a way, would be to create one injury variable with 5 levels: "limb", "thoracic", "abdominal", "two of them", and "all of them". This way, each participant appears only once and interpretation is straightforward, though maybe not as nuanced as you need.

Another way might be having three binary variables, one for each injury type, each coded for "yes" for injury present, "no" for injury not present. Then, you could predict the outcome with these 3 variables and their interactions (if you have enough data for testing interactions).

If you are specifically interested in the role of thoracic injury, maybe a variable with three levels: "only thoracic injury", "thoracic injury plus other injury", "non-thoracic injury/injuries" would work best.

Edited to add: these types of variables can then be used as predictors in a logistic regression.