Aggregating Standard Errors for Predicted Probability Estimates

I obtain predicted values from a logistic regression for a certain outcome (e.g., mortality) at the hospital level – the data is at the patient level – and need to compute the average across hospitals.

The model for the logistic regression is:

$$Y_{ih} = \beta_{o} + \beta_{1}X_{ih} + U_h +\epsilon$$

The dependent variable, $Y_{ih}$ is the outcome of interest for individual $i$ receiving services at hospital $h$. In the equation, $U_h$ are hospital fixed effects, and $X_{ih}$ includes patient characteristics (i.e., age, gender, race, etc.).

The predicted probabilities are estimates of $Y$ that hold independent variables at their means with the exception of which is a vector of dummies for each hospital in the sample. NOTE: I actually implement this in stata and get these predicted probabilities from the "margin hospital, atmeans" command and the standard error for these estimates come from the estimated variance matrix. Assume there are 30 hospitals in the sample. Thus, for any given hospital, the dummy variable for that hospital will have the value of 1 and all other hospital dummies will have the value of 0 when obtaining the predicted probability of outcome Y.

So, after running this model, I have “probability” of mortality at hospital h and a standard error for this estimate (coming from the variances of regression estimates), and the respective sample sizes for each hospital.

These 30 hospitals can be classified as belonging in one of two groups, Group A (20 hospitals) and Group B (10 hospitals). I want to obtain the average predicted probability and standard error at the group level.

The average is straightforward as it is a weighted average of hospital-level predicted probabilities. I am not sure I am computing the standard error of this average of estimates correctly. Currently, I

1. Obtain the hospital-level variance of the predicated probability by squaring the hospital-level standard error : $se_{h}= \sqrt{\hat\sigma^2}$
2. Then, I “weight” the estimated variance by each hospital’s sample size : $\frac{\hat\sigma^2}{n_{h}}$ for each hospital
3. Then, I sum the weighted variances for all hospitals in a group and take the square root of that sum to obtain the group-level standard error

Is this the correct way of determining the standard error?

Here are a few other points to consider:

• The hospital samples are independent
• I ran a variation of the model above that replaced the hospital fixed effects with group-level (A vs B) fixed effects. The predicted probability from this model nearly matched the weighted average of predicted probabilities above (about 0.03 percentage point difference for outcome X). The standard error for the group-level predicted probability from this model - for outcome X for example - was 0.3807. The standard error from the hospital level prediction (generated from the three step process above) was closer to 0.02753. I would not expect the standard errors to be the same but I would not expect them to be so different.
• The sample size for this patient level data is just under 1million obs for any given measure and sample sizes, though different, are fairly similar across hospitals.