Adjusted mortality rates

Question

I have some data similar to this

died  age  hospital
0     75     AA
0     88     AA
1     81     AA
0     77     AA
1     65     AA
0     41     AA
0     66     BA
1     81     BA
0     82     BA
1     64     BA
0     65     BA
1     52     BA

I was asked to calculate "age adjusted mortality rates" for each hospital. There are around 150 hospitals and approx 1000 patients (observations) per hospital. Each row in the data concerns a particular patient.

I was told how this could be done in Stata:

• Perform logistic regression of died on age.
• Use the predict function to get patient-level probabilities of death.
• Summarise the patient-level probabilities by hospital to get the mortality rates for each hospital.

However, I am using R.

Is this the correct approach ? Are the alternatives ? Can I do the same thing in R with glm and predict ?

Edit:

I should perhaps add that there are several other variables that are going to be adjusted for in the model. I have shown age above, just for simplicity.

Answer 1

Seems right. Here are the codes for replication.

died     <- c(0,0,1,0,1,0,0,1,0,1,0,1)
age      <- c(75,88,81,77,65,41,66,81,82,64,65,52)
hospital <- c(rep("AA",6),rep("BA",6))
mydata   <- data.frame(died, age, hospital)
m01 <- glm(died~age, family=binomial)

Now, for prediction:

> predict(m01)
           1            2            3            4            5            6 
-0.398530343 -0.547998721 -0.467515748 -0.421525478 -0.283554668 -0.007613048 
           7            8            9           10           11           12 
-0.295052236 -0.467515748 -0.479013316 -0.272057101 -0.283554668 -0.134086291 

which gives predicted ln(odds). And

> predict(m01, type="response")
        1         2         3         4         5         6         7         8 
0.4016655 0.3663288 0.3852044 0.3961518 0.4295825 0.4980967 0.4267674 0.3852044 
        9        10        11        12 
0.3824851 0.4324021 0.4295825 0.4665286 

gives predicted probability.

Stata says the same thing (see instruction I used and image I got)

But the more important comment is that in epidemiology, "age-adjustment" can also mean rates that are standardized by a certain reference population make up. You said "you were asked" and "you were told to do that in Stata." If these are told by the same person then I think you're fine. If they are two different people you should go back and clarify if (s)he meant epidemiological adjustment or statistical adjustment.

Also, analytically... I feel that if hospital care is a main factor, then we shouldn't go into the picture assume that the mortality risk is even across hospitals. I think it may make more sense to do the prediction by hospital, and then apply an imaginary patient body across each of the hospitals, and the come up with their relative death rates. This approach seems more "age-adjusted." Alternatively, your other predictors should better include hospital characteristic or even the hospital variable itself.

Answer 2

Well, I've worked it out for myself, though I'm still a little unsure if it is the correct approach.

fit1 <- glm(died~age,family=binomial(link=logit),data=dt)
pr <- predict(fit1, type="response")
dt <- cbind(dt,pr)

hosp.mortality.rates <- by(dt$pr, dt$hospital, mean)