Difference in Outcome among treatment groups, with a binary control variable

[Sarah F]

How do I check for differences in outcome (ordinal) among 2 treatment groups, while controlling for gender? Also, each treatment group was replicated 3 times, how does that affect the analysis?

[Clyde Schechter]

Please clarify what you mean by "Also, each treatment group was replicated 3 times." Do you mean that you actually have 6 treatment groups, 3 receiving one intervention and 3 receiving another? If so, what distinguishes the treatment groups that received each intervention: time period, location, person or other entity administering the intervention, person or other entity evaluating the outcome, or anything else that could potentially influence the outcome?

If you meant the outcome was measured 3 times in each person/firm/other unit of analysis, are all three measurements carried out in the same way, or are there relevant factors that differ among the three replicate measurements? And also, was one of these measurements undertaken as a "baseline" before the interventions were actually begun?

[Sarah F]

Hello Clyde,

Apologies for the late reply, did not get a notification for your message! Yes, i have 6 treatment groups where 3 received the control and 3 received the treatment. All 6 treatments were administered at the same time, by the same people just in different chicken 'houses'.

[Clyde Schechter]

OK, that clarifies a lot. The remaining large question in my mind is whether there is reason to think that the effect of the treatment might differ among the three chicken houses. That, of course, is a substantive question of the science, not a statistical question, and depends on what the intervention and outcome are and how the treatment is thought to work. So I'll leave that to you, as I don't know much about chickens. The simpler case is if there is no reason to think the effect would differ among the houses. In that case, your basic analysis would be

ologit outcome i.treatment i.gender i.chicken_house
margins treatment
margins, dydx(treatment)

The -ologit- is a simple regression model for ordinal outcomes. As coded there it adjusts for gender differences in outcome and also for any chicken_house differences in outcome (but not for chicken_house differences in treatment effect). The first -margins- command will give you the predicted probabilities of each outcome level in each treatment group (on average), and the second will give you the expected effect of treatment on those probabilities.

If, however, there is reason to think the treatment effect is not necessarily the same in all three chicken houses, then it's a bit more complicated:

ologit outcome i.treatment##i.chicken_house i.gender
testparm i.treatment#i.chicken_house
margins treatment#chicken_house
margins, dydx(treatment) at(chicken_house = ( 1 2 3 ))
margins, dydx(treatment)

Here we have an interaction model that estimates a separate treatment effect in each chicken house. The first thing you should do as soon as you run the -ologit- is look at the coefficients of the interaction terms, to help decide whether the treatment effects actually do differ among the three houses. If the interaction coefficients are close to zero, then the effect of the treatment actually doesn't differ (much) across the chicken houses and you can consider reverting to the simpler model shown at the top. If you want a test of that, the -testparm- command shows you a p-value for the null hypothesis that the effect of treatment is the same in all three houses. (Interpret this hypothesis test with the same caution you would and should apply to any null hypothesis test--how large are the effects regardless of significance issues? Do you have adequate sample size to power this test? etc.)

Assuming that the data have supported the notion that the treatment effect differs among the houses, proceed with the -margins- commands for more detailed output. The first -margins- command will show you the predicted probabilities of the outcomes in each treatment in each house. (A lot of numbers: 3 outcome levels X 2 treatment groups X 3 houses). The second -margins- command will give you the estimated effect of treatment on those probabilities in each house. And the third command will give you a single estimate of the treatment effect averaged over all three houses. People disagree about the usefulness of that last one--personally I'm not a fan, but some people like it.

[Sarah F]

Thanks a lot Clyde!

There is actually no reason to believe that the effect of the treatment might differ among the chicken houses. Running the first regression model yields:


ologit lesion i.treatment i.gender i.house [fw=count], or

note: 6.house omitted because of collinearity
Iteration 0: log likelihood = -230.6573
Iteration 1: log likelihood = -213.22468
Iteration 2: log likelihood = -212.74664
Iteration 3: log likelihood = -212.74499
Iteration 4: log likelihood = -212.74499

Ordered logistic regression Number of obs = 300
LR chi2(6) = 35.82
Prob > chi2 = 0.0000
Log likelihood = -212.74499 Pseudo R2 = 0.0777

------------------------------------------------------------------------------
lesion | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
2.treatment | .2411479 .1065256 -3.22 0.001 .1014543 .5731873
2.gender | .7621192 .1935357 -1.07 0.285 .4633031 1.253662
|
house |
2 | .9208921 .3613151 -0.21 0.834 .4268123 1.986921
3 | .8305073 .3284217 -0.47 0.639 .3825941 1.802804
4 | 1.850332 .8362245 1.36 0.173 .7630696 4.486782
5 | .4867716 .269044 -1.30 0.193 .1647619 1.438115
6 | 1 (omitted)
-------------+----------------------------------------------------------------
/cut1 | -.297631 .3202051 -.9252215 .3299595
/cut2 | 2.673613 .4242624 1.842074 3.505152
------------------------------------------------------------------------------



which shows that only treatment is significant.

However, the margins command yields:

margins treatment [fw=count]

Predictive margins Number of obs = 300
Model VCE : OIM

Expression : Pr(lesion==0), predict()

------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
treatment |
1 | . (not estimable)
2 | . (not estimable)
------------------------------------------------------------------------------

Did i do anything wrong?

[Sarah F]

I tried omitting gender and house from the model since they were not significant, than ran margins again:

. ologit lesion i.treatment [fw=count], or

Iteration 0: log likelihood = -230.6573
Iteration 1: log likelihood = -217.12517
Iteration 2: log likelihood = -216.99049
Iteration 3: log likelihood = -216.99042
Iteration 4: log likelihood = -216.99042

Ordered logistic regression Number of obs = 300
LR chi2(1) = 27.33
Prob > chi2 = 0.0000
Log likelihood = -216.99042 Pseudo R2 = 0.0593

------------------------------------------------------------------------------
lesion | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
2.treatment | .2774711 .0702242 -5.07 0.000 .1689625 .4556648
-------------+----------------------------------------------------------------
/cut1 | -.0439427 .1631093 -.363631 .2757457
/cut2 | 2.903306 .3311728 2.254219 3.552392
------------------------------------------------------------------------------

. margins treatment

Adjusted predictions Number of obs = 300
Model VCE : OIM

Expression : Pr(lesion==0), predict()

------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
treatment |
1 | .4890161 .0407576 12.00 0.000 .4091326 .5688996
2 | .7752324 .0338878 22.88 0.000 .7088134 .8416513
------------------------------------------------------------------------------

. margins, dydx(treatment)

Conditional marginal effects Number of obs = 300
Model VCE : OIM

Expression : Pr(lesion==0), predict()
dy/dx w.r.t. : 2.treatment

------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
2.treatment | .2862163 .0528512 5.42 0.000 .1826298 .3898027
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.


How do i interpret the results now?

[Clyde Schechter]

OK. The inestimability you initially encountered usually results from empty cells. That is, there is probably some combination of treatment, gender, and house that does not occur in your data. In that case, the predictive margins are not fully identified and the results, were you to force Stata to give them, would be sensitive to the particular parameterization of the model (including things like which levels of the variables were used as the reference categories.)

I don't like the approach of removing covariates that are "not significant." Statistical significance is not really the point when adjusting for the effects of potential confounders. You may notice that chicken houses 4 and 5 both have odds ratios that are markedly different from 1. Associations that large cannot be casually dismissed just because they don't have a p-value less than 0.05. Even the gender odds ratio is pretty impressive. These variables are clearly moderately to strongly associated with the outcome. The non-significance reflects the smaller sample sizes when breaking things down by gender and chicken house--and p-values are irrelevant to the issue of confounding bias in any case. The real test, however, is whether the treatment odds ratio changes much when you remove those variables. Fortunately for you, it doesn't. It only went from 0.24 to 0.28. So there is no material issue with removing gender and chicken house from the model here. (In other words, although the p-value was the wrong question to ask, you got the right answer anyway.)

So, your headline result is that the odds ratio associated with treatment is 0.28 (95% CI 0.17-0.466). (You probably also want to report the adjusted odds ratio from the model that includes chicken house and gender. If treatment assignment was randomized, then I think it's not necessary, but if treatment was not randomized, then I think it's mandatory to show the adjusted odds ratio as well, as the crude odds ratio is not very credible in a non-randomized study.) Unfortunately, your -margins- output does not look like what it should. I'm guessing that you are using an older version of Stata (current version is 14). In current Stata, you would get predicted probablities for each level of your outcome variable--but you have them only for one level I don't have access to an older version of Stata at this point, so I'm working from memory here, but I think that you can control which outcome level you get probabilities for by specifying -predict(outcome(#))- [replacing # by the desired level] as an option in the -margins- command. But you should read -help margins- to figure this out. Anyway, it seems you will have to run a separate -margins- command for each outcome level. (If you had more than 3 levels I'd recommend doing that in a loop, but for three levels it's just as easy and more transparent to spell out three -margins- commands.) You will then have three outputs that look more or less like the one you have now, with different numbers.

As for interpreting those, the output you have states that it refers to the 0 level for outcome lesions. So Stata is saying that averaging over the entire sample, treatment level 1 is associated with a 0.48 probability (95% CI 0.41-0.57) of outcome lesion = 0, whereas treatment level 2 is associated with 0.78 probability of the same outcome (95% CI 0.71-0.84). The output of the -margins, dydx(treatment)- command reports that being in the treatment = 2 condition is associated with an increase of 0.29 (95% CI 0.18-0.39) in the probability that outcome lesion will take on the 0 value compared to the treatment = 1 condition.

When you run the -margins- commands for the other outcome levels you will get additional similar outcome that is interpreted in the same way.

By the way, in the future, please paste your Stata output into code blocks so that it will align nicely and be easier to read. If you don't know about code blocks, read FAQ #12.

[Sarah F]

Thank you Clyde, very grateful for your input and help. I am indeed using an older version, 12.1. Will read up on code blocks also...

[Sarah F]

. ologit lesion i.treatment i.gender [fw=count], or

Iteration 0:   log likelihood = -82.377278  
Iteration 1:   log likelihood = -73.296479  
Iteration 2:   log likelihood = -70.197648  
Iteration 3:   log likelihood = -70.140164  
Iteration 4:   log likelihood = -70.140092  
Iteration 5:   log likelihood = -70.140092  

Ordered logistic regression                       Number of obs   =        381
                                                  LR chi2(2)      =      24.47
                                                  Prob > chi2     =     0.0000
Log likelihood = -70.140092                       Pseudo R2       =     0.1486

------------------------------------------------------------------------------
      lesion | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
 2.treatment |   .1294484   .0832427    -3.18   0.001      .036705    .4565287
    2.gender |   .1887607   .0964155    -3.26   0.001     .0693639    .5136763
-------------+----------------------------------------------------------------
       /cut1 |   1.428942   .3171903                      .8072602    2.050623
       /cut2 |   4.557518   1.020984                      2.556426    6.558609
------------------------------------------------------------------------------

. margins treatment

Predictive margins                                Number of obs   =        381
Model VCE    : OIM

Expression   : Pr(lesion==1), predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   treatment |
          1  |   .9027994   .0217048    41.59   0.000     .8602588      .94534
          2  |   .9854575   .0082985   118.75   0.000     .9691928    1.001722
------------------------------------------------------------------------------

. margins, dydx(treatment)

Average marginal effects                          Number of obs   =        381
Model VCE    : OIM

Expression   : Pr(lesion==1), predict()
dy/dx w.r.t. : 2.treatment

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
 2.treatment |   .0826581   .0232437     3.56   0.000     .0371013     .128215
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.

.

For this model, both treatment and gender were significant. Can i run margins for treatment given gender is in the model? If yes, how do I interpret the results?

[Clyde Schechter]

Can i run margins for treatment given gender is in the model? If yes, how do I interpret the results?

Yes. Your -margins treatment- output indicates that if everyone in your sample were to receive treatment 1, the model-predicted probability of lesion = 1 would be about 0.903, whereas if everyone were to receive treatment 2, it would be about 0.985. In calculating these predicted probabilities, Stata assumes that everybody has the same value of gender that is actually in your sample data. So, this is a predicted prevalence (or incidence as the case may be) of lesion = 1 conditional on treatment adjusted to the observed gender distribution of your sample. The output of -margins, dydx(treatment)- indicates that the expected difference in probability of lesion = 1 between the treatment = 2 and the treatment != 2 conditions is about 0.083, again adjusted to the observed gender distribution in your sample.

I notice that your outcome variable is ordinal. If you want the predicted probabilities or marginal effects on those probabilities for the other levels of your outcome, you will have to run a separate command for each level and specify the level in each command. (In Stata version 14, the default behavior of -margins- is to give results for all of the outcome levels, but back in Stata 12, it was one at a time, with level 1 being the default if no level was specified.)

[Sarah F]

Thank you very much Clyde, always very helpful!