Using standardisation to adjust means

Question

I'm familiar with direct and indirect standardisation to calculate adjusted rates and standardised ratios (e.g. SMRs).

Is it valid to use similar techniques to adjust means?

I'm interested in the duration of hospital admissions for a specific condition and want to compare the mean duration of two groups with different age profiles.

The indirect approach would be:

+-----------+----------------------+----------------+---------------+----------------------+----------------+---------------+-------------------+---------------+
|           | REFERENCE POPULATION                                  |   STUDY POPULATION                                                                        |
+-----------+----------------------+----------------+---------------+----------------------+----------------+---------------+-------------------+---------------+
| Age group | Number of admissions | Total bed-days | Mean bed-days | Number of admissions | Total bed-days | Mean bed-days | Expected bed-days | Expected mean |
| 15-34     | 1000                 | 1500           | 1.5           | 300                  | 300            | 1.0           | 450               | 1.5           |
| 35-44     | 2000                 | 4000           | 2.0           | 350                  | 525            | 1.5           | 700               | 2.0           |
| 45-55     | 3000                 | 7500           | 2.5           | 400                  | 800            | 2.0           | 1000              | 2.5           |
| TOTAL     | 6000                 | 13000          | 2.2           | 1050                 | 1625           | 1.5           | 2150              | 2.0           |
+-----------+----------------------+----------------+---------------+----------------------+----------------+---------------+-------------------+---------------+

Which would allow us to calculate a 'standardised mean ratio' of 1.5 / 2.0 * 100 = 80.

The direct approach would be:

+-----------+----------------------+----------------------+-------------------------+
|           |   STUDY POPULATION   | REFERENCE POPULATION                           |
+-----------+----------------------+----------------------+-------------------------+
| Age group | Actual mean bed-days | Number of admissions | Expected total bed-days |
| 15-34     | 1                    | 1000                 | 1000                    |
| 35-44     | 1.5                  | 2000                 | 3000                    |
| 45-55     | 2                    | 3000                 | 6000                    |
| TOTAL     |                      | 6000                 | 10000                   |
+-----------+----------------------+----------------------+-------------------------+

Which would allow us to calculate an 'age adjusted mean' of 10000/6000 = 1.7 days. Which is higher than the actual mean of 1.5 days in the study population, because the study population is younger and younger patients tend to have shorter stays.

Is that valid?

And if it is valid, do you know any method of estimating confidence intervals around the 'age adjusted mean'? I have got the raw data.

Thanks

Answer 1

OK... I think I've worked it out... for anyone else who happens to be interested...

I believe the age-adjusted mean is a valid measure. The 'Concise Encyclopedia of Biostatistics for Medical Professionals' describes it on p5 in the same way I've described the direct approach above.

In R, the 'survey' package includes a function 'svystandardize' that I believe can be used to construct standard errors of age-adjusted means. Here's an example I've created:

# generate data in which older people and women have longer spell durations, but men are older

n <- 250
datm <- data.frame(agegrp = sample(c('15-34', '35-44', '45-55'), n, T, c(0.1, 0.3, 0.6)), sex = rep('m', 50))
datf <- data.frame(agegrp = sample(c('15-34', '35-44', '45-55'), n, T, c(0.6, 0.3, 0.1)), sex = rep('f', 50))
dat <- rbind(datm, datf)
dat$agegrp <- as.factor(dat$agegrp)
dat$spell <- runif(n * 2, 0, 5)
dat$spell <- dat$spell * ifelse(dat$agegrp == '35-44', 1.5, 1)
dat$spell <- dat$spell * ifelse(dat$agegrp == '45-55', 2.0, 1)
dat$spell <- dat$spell * ifelse(dat$sex == 'f', 1.3, 1)

# calculate raw means: men higher than women; women higher than men within strata

aggregate(dat$spell, by = list(dat$sex), FUN = 'mean')
x <- aggregate(dat$spell, by = list(dat$agegrp, dat$sex), FUN = 'mean')

# calculate age-adjusted means: women higher than men

popage <- as.numeric(table(dat$agegrp))
sum(popage * x[1:3,3]) / sum(popage) # male age-adjusted mean
sum(popage * x[4:6,3]) / sum(popage) # female age-adjusted mean

# standard errors of age-adjusted means

library(survey)

design <- svydesign(ids = ~1, strata = ~agegrp, data = dat)
stdes <- svystandardize(design, by = ~agegrp, over = ~sex, population = popage)
y <- svyby(~spell, ~sex, svymean, design = stdes)
y$lower <- y$spell - 1.96 * y$se
y$upper <- y$spell + 1.96 * y$se

y$spell[2] / y$spell[1]

The final line shows that the age-adjusted mean spell for women is approx 1.3 times that of men, as you would expect from the data. The object y includes the age-adjusted means and confidence intervals.